GF2::LUDecompose — Bit-matrix Inverse

We can use the LU decomposition of \(A\) to find \(A^{-1}\)

std::optional<Matrix> invert() const;   (1)
1 Returns std::nullopt if \(A\) is singular otherwise returns \(A^{-1}\).
Example
#include <GF2/GF2.h>
int
main()
{
    // Number of trials
    std::size_t trials = 100;

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

    // 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 & decompose it
        auto A = GF2::Matrix<>::random(N);
        auto LU = GF2::lu_decompose(A);

        // See if we can invert the matrix and if so check A.A_inv == I
        if (auto A_inv = LU.invert(); A_inv) {
            auto I = GF2::dot(A, *A_inv);
            std::cout << "A.Inverse[A] == I? " << (I.is_identity() ? "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)
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES
A.Inverse[A] == I? YES

Singularity stats ...
Matrix size: 30 x 30
P[singular]: 71.1212%
Trials:      100
Singular:    68 times
Expected:    71 times
See Also

singular