bit::lu — Solutions

We can use the LU decomposition of \(A\) to solve the system \(A \cdot x = b\):

1std::optional<bit::vector> operator()(const bit::vector &b) const;
2std::optional<bit::matrix> operator()(const bit::matrix &B) const;
1
If \(A\) is non-singular, this solves the system \(A \cdot x = b\).
2
If \(A\) is non-singular, this solves the systems \(A \cdot X = B\).

In the second case, each column of the bit-matrix B is considered a separate right-hand side, and the corresponding column of \(X\) is the solution vector.

Once you have the LU decomposition of \(A\), it is easy to solve systems like these. If \(A\) is \(n \times n\) each system solution takes just \(\mathcal{O}(n^2)\) operations.

These methods return std::nullopt if the underlying bit-matrix \(A\) is singular. You can avoid that by first calling the lu::singular method.

Both methods throw an exception if the number of elements in \(b\) or rows in \(B\) does not match the number of rows in \(A\). They could instead return a std::nullopt, but a dimension mismatch is likely an indication of a coding error somewhere.

Example

#include <bit/bit.h>
int
main()
{
    // Number of trials
    std::size_t trials = 32;

    // Each trial will run on a bit-matrix of this size
    std::size_t N = 16;

    // Number of non-singular matrices
    std::size_t singular = 0;

    // Start the trials
    for (std::size_t n = 0; n < trials; ++n) {

        // Create a random matrix & vector
        auto A = bit::matrix<>::random(N);
        auto b = bit::vector<>::random(N);

        // LU decompose the matrix & solve A.x = b
        auto LU = bit::lu(A);
        if (auto x = LU(b); x) {
            auto Ax = bit::dot(A, *x);
            std::cout << "x: "   << x->to_string() << "; ";
            std::cout << "A.x: " << Ax.to_string() << "; ";
            std::cout << "b: "   << b.to_string()  << "; ";
            std::cout << "A.x == b? " << (Ax == b ? "YES" : "NO") << "\n";
        }

        // Count the number of singular matrices we come across
        if (LU.singular()) singular++;
    }

    // Stats
1    auto p = bit::matrix<>::probability_singular(N);
    std::cout << "\n"
              << "Singularity stats ...\n";
    std::cout << "bit::matrix size: " << N << " x " << N << "\n"
              << "P[singular]: " << 100 * p << "%\n"
              << "Trials:      " << trials << "\n"
              << "Singular:    " << singular << " times\n"
              << "Expected:    " << int(p * double(trials)) << " times\n";
    return 0;
}
1
See matrix::probability_singular

Output for a consistent system (details depend on the values of the random inputs)

x: 0100010101110000; A.x: 0101011111010111; b: 0101011111010111; A.x == b? YES
x: 0110111000000101; A.x: 0001100110100101; b: 0001100110100101; A.x == b? YES
x: 1001000000111000; A.x: 0111110110111101; b: 0111110110111101; A.x == b? YES
x: 1011010000110100; A.x: 0100001001010100; b: 0100001001010100; A.x == b? YES
x: 0110100100110100; A.x: 1001111001100001; b: 1001111001100001; A.x == b? YES
x: 0101000101111100; A.x: 1001100000011101; b: 1001100000011101; A.x == b? YES
x: 0110000100100100; A.x: 0010100110010110; b: 0010100110010110; A.x == b? YES
x: 1011001101010000; A.x: 0010011101110000; b: 0010011101110000; A.x == b? YES
x: 1101101110001111; A.x: 0011010110110010; b: 0011010110110010; A.x == b? YES
x: 0110101001101110; A.x: 1011010001011010; b: 1011010001011010; A.x == b? YES
x: 1000011100010001; A.x: 0100111110001101; b: 0100111110001101; A.x == b? YES

Singularity stats
bit::matrix size: 16 x 16
P[singular]: 71.1207%
Trials:      32
Singular:    21 times
Expected:    22 times

See Also

lu::singular

Back to top