GF2::LUDecompose — Solutions Access

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

std::optional<Vector> operator()(const Vector &b) const;  (1)
std::optional<Matrix> operator()(const Matrix &B) const;  (2)
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 as a separate right hand side and the corresponding column of \(X\) is the solution vector.

Once you have the LU decomposition of \(A\), solving systems like these is relatively painless. 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 is_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 of course instead just return a std::nullopt but a dimension mismatch is almost surely an indication of a coding error somewhere.
Example
#include <GF2/GF2.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 = GF2::Matrix<>::random(N);
        auto b = GF2::Vector<>::random(N);

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

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

    // Stats ...
    auto p = GF2::Matrix<>::probability_singular(N);        (1)
    std::cout << "\n"
              << "Singularity stats ...\n";
    std::cout << "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 probability_singular
Output for a consistent system (details depends 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 ...
Matrix size: 16 x 16
P[singular]: 71.1207%
Trials:      32
Singular:    21 times
Expected:    22 times
See Also

singular