On Invertibility of Large Binary Matrices

1. Introduction

The 1969 paper by Strassen [1] showed that matrix multiplication can be done with a sub-cubic complexity. More recent algorithms such as the Method of Four Russians and the Strassen-Winograd multiplication as well as various low-level strategies such as the bit-level parallelization for multiplying dense binary matrices can be found in [2]. The Schur complement is often a basis for blockwise recursive methods to invert large real-valued matrices [3]. For base cases of the recursive methods, say, the matrices of sizes five or less, it is more efficient to use the direct Gauss-Jordan elimination (GJE) [4]. The caveat is that the submatrices created during the recursion must themselves be invertible. Unfortunately, for the matrices defined over finite Galois fields, this cannot be guaranteed, and it is also the key reason why the recursion often fails in practice.

In some applications, it is useful to explore the high-throughput strategies for enumerating the invertible matrices of a given size. The invertible matrices of size $(n\times n)$ over a Galois field, q, form the general linear group $GL(n,q)$[5]. The GL group can be ordered using a Grey encoding, so that any two subsequent matrices are related by simple operations on rows [6]. The underlying group structure can be studied using the corresponding Cayley graphs.

The current literature contains mainly the results on optimizing the primitives in matrix multiplications, and also provides deep theoretical insights into the structure of the $GL(n,q)$ groups. However, to our best knowledge, there are presently no investigations how to efficiently invert very large matrices over finite fields, and the related problems. In particular, it is sometimes sufficient to only detect when such a matrix is singular. In other cases, it may be necessary to also find the closest matrix that is invertible. Moreover, any recursive algorithm for inverting the matrices over finite fields needs to address the issue of submatrices being singular.

This paper fills the gap in the literature by reporting a robust and fast algorithm for inverting very large binary matrices. It is based on the pivoted PLU factorization [4]. On the other hand, the efficient enumeration of all invertible matrices in the $GL(n,2)$ group is achieved by the Bruhat decomposition [7]. In addition, a method for finding an invertible matrix that is closest to a given binary singular matrix is devised. This method determines the positions of the minimum number of bit-flips to make the matrix non-singular. The developed algorithms are implemented in C++, and they are available from the Github repository [8]. All implemented algorithms were validated. Even though the focus in this paper is on binary matrices, the proposed methods are also effective for the matrices defined over other finite fields by updating the arithmetic operations to these fields.

This paper is organized as follows. The rest of this section summarizes the terminology and notations, which are important in understating the paper. Section 2 introduces a robust blockwise matrix inversion algorithm based on the PLU factorization. The algorithm for efficiently enumerating elements of general linear groups is presented in Section 3. It is based on the Bruhat group decomposition. An algorithm for repairing singular binary matrices to make them invertible by performing the minimum number of bit-flips is proposed in Section 4. The paper is concluded in Section 5. The examples illustrating the proposed algorithms are provided in Appendix A and B.

1.1. Preliminaries and Notations

The Galois field, $\mathrm{GF}\left(2\right)$, is a ring, $\{0,1\}$, with addition denoted by symbol, ⊕, and multiplication denoted by symbol, ∧. Additions over $\mathrm{GF}\left(2\right)$ can be effectively performed using logical exclusive-or (XOR), whereas multiplications are equivalent to logical AND. Note also that, in $\mathrm{GF}\left(2\right)$, the addition and subtraction are equivalent.

The matrices are denoted by bold capital letters. The symbol, $\mathit{A}(i,j)$, represents the i-th row and the j-th column element of the matrix, $\mathit{A}$. A bit-flip at position, $(i,j)$, is defined as, $\mathit{A}(i,j):=\mathit{A}(i,j)\oplus 1$. The $(n\times n)$ identity matrix is denoted by ${\mathit{I}}_n$.

The $(n\times n)$ permutation matrix, $\mathit{P}$, is created from the identity matrix, ${\mathit{I}}_n$, by permuting its row and/or columns. There are $n!$ (n factorial) unique permutation matrices. The permutation matrices are orthogonal, i.e., $\mathit{P}{\mathit{P}}^T={\mathit{P}}^T\mathit{P}=\mathit{I}$.

A $(n\times n)$ unit lower-triangular binary matrix, $\mathit{L}$, have all diagonal entries equal to 1, while all the entries above the main diagonal are 0. Similarly, a $(n\times n)$ upper-triangular binary matrix, $\mathit{U}$, have all diagonal entries equal to 1, and all the entries below the main diagonal are 0.

The rank of $(n\times n)$ matrix, $\mathit{A}$, is the maximum number of linearly independent rows (or, equivalently, columns). The matrix has full rank, if $rank\mathit{A}=n$. The square matrix is invertible, i.e., its inverse, ${\mathit{A}}^{-1}$, exists, if and only if, it has full rank; otherwise it is said to be singular. The singular matrix has the rank deficiency, $d=n-rank\mathit{A}$. The rank deficiency represents how many linearly independent rows or columns are missing to have a full span of the n-dimensional vector space.

A general linear group $GL(n,2)$ is the set of all invertible $(n\times n)$ binary matrices that is closed under matrix multiplication. The unit element of the group is the identity matrix, ${\mathit{I}}_n$. The group has order (i.e., the cardinality),

$$|GL(n,2)|=\prod _{k=0}^{n-1}(2^n-2^k).$$

(1)

The $(n\times n)$, n is even, matrix can be partitioned into four equal-size submatrices (i.e., blocks) as,

$$\mathit{A}=\left(\begin{array}{cc}{\mathit{A}}_{11}& {\mathit{A}}_{12}\\ {\mathit{A}}_{21}& {\mathit{A}}_{22}\end{array}\right).$$

(2)

Provided that ${\mathit{A}}_{11}$ is invertible, its Schur complement is defined as, $\mathit{S}={\mathit{A}}_{22}-{\mathit{A}}_{21}{\mathit{A}}_{11}^{-1}{\mathit{A}}_{12}$.

2. Accelerated Inversion of Binary Matrices

The objective is to devise a fast algorithm for the robust inversion of very large binary matrices. The proposed algorithm relies on the following crucial techniques. The pivoted PLU factorization enables the proper factorization of any invertible matrix. Importantly, it avoids the fundamental limitation of simpler blockwise inversion methods that often fails at singular submatrices for matrices defined over finite fields. The recursive blockwise inversion can be further optimized for speed by considering the sub-cubic Strassen’s algorithm for matrix multiplication in order to improve the asymptotic time complexity.

For smaller matrices, the recursion overhead exceeds the complexity of the direct GJE. The GJE assumes the augmented matrix, $\left[\mathit{A}\right|\mathit{I}]$, and performs the row operations until it arrives at the matrix, $\left[\mathit{I}\right|{\mathit{A}}^{-1}]$. In $\mathrm{GF}\left(2\right)$, it only involves the row swaps and XORing the rows. If the matrix is found to be singular, i.e., there is no pivot for a given matrix column, the algorithm correctly reports a failure. The cubic time complexity of the GJE is acceptable for smaller matrices.

For larger matrices, other algorithm strategies such as divide-and-conquer and dynamic programming are required. These strategies provide much better control over the time complexity as the problem size increases. In this paper, we first consider a standard blockwise recursive inversion with the Schur complement to show that it is fundamentally flawed when considering matrices over finite fields. This will motivate much more robust block-recursive inversion method that is based on the PLU factorization. Recall also that, in this paper, only the binary matrices are considered; generalizing the methods to other finite fields should be straightforward by modifying the underlying arithmetic operations.

The standard blockwise inversion of matrices assumes the Schur complement to reduce the $(n\times n)$ problem into the four sub-problems, each of size $(n/2)\times (n/2)$. In particular, assuming that both ${\mathit{A}}_{11}$ and its Schur complement $\mathit{S}$ are non-singular, the inverse ${\mathit{A}}^{-1}$ is computed recursively using its Schur complement. The Schur complement requires two matrix multiplications, and its inverse, ${\mathit{S}}^{-1}$, is again computed recursively. The inverse matrix is then constructed as,

$${\mathit{A}}^{-1}=\left(\begin{array}{cc}{\mathbf{B}}_{11}& {\mathbf{B}}_{12}\\ {\mathbf{B}}_{21}& {\mathbf{B}}_{22}\end{array}\right)$$

(3)

where

$${\mathit{B}}_{22}={\mathit{S}}^{-1},{\mathit{B}}_{12}={\mathit{A}}_{11}^{-1}{\mathit{A}}_{12}{\mathit{S}}^{-1},{\mathit{B}}_{21}={\mathit{S}}^{-1}{\mathit{A}}_{21}{\mathit{A}}_{11}^{-1},and,{\mathit{B}}_{11}={\mathit{A}}_{11}^{-1}\oplus {\mathit{B}}_{12}{\mathit{A}}_{21}{\mathit{A}}_{11}^{-1}.$$

(4)

Computing the matrices in (4) requires six matrix multiplications.

Furthermore, in order to improve the asymptotic time complexity of the standard matrix blockwise inversion, the six required matrix multiplications can be efficiently performed using the Strassen’s technique, which has a sub-cubic time complexity. Subsequently, only seven rather than eight matrix multiplications are required in computing each sub-problem in the divide-and-conquer matrix inversion method described above; the asymptotic time complexity is then reduced from $\mathcal{O}\left(n^3\right)$ to $\mathcal{O}\left(n^{{log}_{2}7}\right)$. However, and importantly, for $\mathrm{GF}\left(2\right)$ and other finite fields, the submatrices, ${\mathit{A}}_{11}$, and, $\mathit{S}$, frequently become singular, which terminates the algorithm, even though the original matrix, $\mathit{A}$, is non-singular. The solution is to combine the Schur complement with the PLU factorization in performing the blockwise matrix inversion as will be described next. We also prove the correctness of this algorithm, and analyze its time complexity.

2.1. Block-Wise Inversion with PLU Factorization

The PLU matrix factorization can be used to achieve a fast matrix inversion as originally proposed in [9]. Specifically, the PLU factorization theorem asserts that, for any non-singular matrix, $\mathit{A}$, there exists a permutation matrix $\mathit{P}$ such that $\mathit{P}\mathit{A}=\mathit{L}\mathit{U}$, where $\mathit{L}$ and $\mathit{U}$ are the unit lower-triangular and and upper-triangular matrices, respectively, representing the LU factorization [4]. The key idea is to first recursively compute the PLU factorization of the matrix, $\mathit{A}$, to facilitate its fast inversion, since inverting a triangular matrix amounts to solving a set of linear equations being already in the reduced row echelon form [4].

In particular, the first step is to recursively compute the PLU factorization of the submatrix, ${\mathit{A}}_{11}$. Thus, provided that, ${\mathit{P}}_1{\mathit{A}}_{11}={\mathit{L}}_1{\mathit{U}}_1$, then the inverse, ${\mathit{A}}^{-1}$, is straightforward to obtain. The Schur complement can be inverted similarly using its PLU factorization, ${\mathit{P}}_2\mathit{S}={\mathit{L}}_2{\mathit{U}}_2$. The inverse matrix, ${\mathit{A}}^{-1}$, is then obtain using the structure (3) and (4). Most importantly, the PLU matrix factorization is guaranteed to always yield the correct result for all non-singular binary matrices, $\mathit{A}$.

The matrix bisection that is performed repeatedly at each recursive step of the blockwise matrix inversion and the Strassen’s multiplication implicitly assume that the original matrix, $\mathit{A}$, has the size, $n=2^i$, for some integer, $i\gg 1$. Even though the power of two size cannot be guaranteed, the original matrix can be always padded with zeros to such a size, i.e.,

$$\tilde{\mathit{A}}=\left(\begin{array}{cc}\mathit{A}& \mathbf{0}\\ \mathbf{0}& {\mathit{I}}_{p-n}\end{array}\right)$$

(5)

where $p=2^{\lceil {log}_{2}n\rceil }$, and $\lceil ·\rceil$ is the ceiling function. The inverse of the $(p\times p)$ matrix, $\tilde{\mathit{A}}$, can be then written as,

$${\tilde{\mathit{A}}}^{-1}=\left(\begin{array}{cc}{\mathit{A}}^{-1}& \mathbf{0}\\ \mathbf{0}& {\mathit{I}}_{p-n}\end{array}\right).$$

(6)

Once the complete factorization, $\mathit{P}\mathit{A}=\mathit{L}\mathit{U}$, is obtained, the matrix inverse is computed as,

$${\mathit{A}}^{-1}={\mathit{U}}^{-1}{\mathit{L}}^{-1}\mathit{P}={\left(\mathit{L}\mathit{U}\right)}^{-1}\mathit{P}.$$

(7)

2.2. Complexity and Correctness of Matrix Inversion via PLU Factorization

Lemma 1 

(Correctness of the PLU-based matrix inversion). The non-singular matrices over finite Galois fields can be correctly inverted using the PLU factorization with blockwise matrix decomposition. Moreover, the singular matrices are always correctly detected, and the algorithm terminates.

(Proof outline) 

The GJE, which is adopted for inverting the base cases, is guaranteed to find the inverse, or to correctly detect that the submatrix is singular. The PLU factorization theorem guarantees that any non-singular matrix, $\mathit{A}$, can be factored as, $\mathit{P}\mathit{A}=\mathit{L}\mathit{U}$,4]. Importantly, the submatrices, ${\mathit{P}}_1{\mathit{A}}_{11}={\mathit{L}}_1{\mathit{U}}_1$, and the Schur complements, ${\mathit{P}}_2\mathit{S}={\mathit{L}}_2{\mathit{U}}_2$, at each step of the blockwise recursion are always invertible by the PLU factorization theorem. Hence, the PLU factorization, $\mathit{P}\mathit{A}=\mathit{L}\mathit{U}$, can be always correctly obtained by properly combining the submatrices, ${\{{\mathit{P}}_i,{\mathit{L}}_i,{\mathit{U}}_i\}}_{i=1,2}$. The unit triangular matrices are always invertible [4]. Consequently, the non-singular matrix, $\mathit{A}$, can be inverted as, ${\mathit{A}}^{-1}={\mathit{U}}^{-1}{\mathit{L}}^{-1}\mathit{P}$. Moreover, the singular matrix is always detected when inverting the base case submatrices as the GJE fails to find all pivots.    □

It should be noted that Lemma 1, and its proof are independent of the specific field. Thus, Lemma 1 is valid for the matrices defined over all finite fields, and even infinite fields including real numbers.

Theorem 1 

(Complexity of blockwise matrix inversion with PLU factorization). The blockwise recursive inversion of $(n\times n)$ matrices with the PLU factorization and Strassen’s matrix multiplication has the complexity, $\mathcal{O}\left(n^{{log}_{2}7}\right)$.

Proof. 

Let $I\left(n\right)$ be the time to invert a $(n\times n)$ matrix, and $M\left(n\right)$ be the time to multiply two $(n\times n)$ matrices. The complexity is dominated by the operations performed at the first blockwise decomposition, i.e., by the PLU factorization and the inversion of the $(n/2)\times (n/2)$ submatrix, ${\mathit{A}}_{11}$, and its Schur complement, $\mathit{S}$. The seven required multiplications of submatrices has the complexity, $\mathcal{O}\left(n^2\right)$. Consequently, the overall complexity is dominated by the two recursive calls, and the matrix multiplications, i.e.,

$$I\left(n\right)=2I(n/2)+c·M(n/2)+O\left(n^2\right)$$

(8)

where c is the total number of required submatrix multiplications. The Strassen’s algorithm for matrix multiplication has the complexity, $M\left(n\right)=\mathcal{O}\left(n^{{log}_{2}7}\right)$. The recursion (8) can be solved by the Master Theorem [10], which yields the total complexity, $I\left(n\right)=\mathcal{O}\left(n^{{log}_{2}7}\right)$. Moreover, since the base cases are assumed to be sufficiently small, they do not affect the overall asymptotic complexity.    □

2.3. Implementation

The blockwise matrix inversion with the PLU factorization and Strassen’s matrix multiplication was implemented in C++ programming language [8]. The implementation adopts several optimizations to improve the actual run-time of inverting large binary matrices. In particular, the low-level C++ data structure, std::bitset<N>, is used to store the rows of binary matrices. This encourages compiled code to exploit the single instruction, multiple data (SIMD) CPU instructions, which provides the bit-level parallelism involving the bit arrays. In turn, this offers a substantial speed-up of adding and multiplying the rows of binary matrices.

The size of submatrices that are treated as base cases, and for which the GJE instead of the blockwise recursion is used, is a hyper-parameter that can be optimized for faster inversion. The former has better data locality, and avoids the overhead of recursive function calls.

The singular matrices, for which the GJE fails to find a pivot when inverting one of the base case submatrices, are detected using the function signature with, std::optional, parameter. This has the benefit that it avoids the complexities of exception handling.

The matrix to be inverted is instantiated at a compile time using the C++ template, template<size_t N>. It provides the type safety, and also reduces the runtime overhead.

The empirical measurements of the runtime and memory usage for inverting the binary matrices of size, $2^n$, $n\le 11$, are shown in Figure 1 and Figure 2, respectively. Both metrics appear to have asymptotically linear growths on the log-log scale plots suggesting the asymptotic complexity, $\mathcal{O}\left(2^{nlog2}\right)$.

Even though the reference implementation in C++ considers specifically the matrices over $\mathrm{GF}\left(2\right)$, the blockwise recursive inversion with PLU factorization can be readily modified to other fields. However, only the computations over the finite Galois fields are exact, whereas the algorithms involving the floating point arithmetic suffer from numerical instabilities and numerically ill-conditioned problems. For example, the pivoting used in the GJE must be made more robust for inverting the real-valued matrices. Moreover, the bit-level parallelization can be much more easily exploited for the finite fields than for real numbers, which has a noticeable impact on the practical runtime and the memory usage of the implemented algorithms.

The matrix multiplications in Strassen’s algorithm can be parallelized to support the implementation of inverting very large matrices on multi-core computing platforms. For such matrices, the data structures could be further optimized to improve the data locality, and to exploit fast memory caching mechanisms. There can be multiple implementations for different values of the matrix sizes, $2^n$. Similar considerations can be used for the high-throughput inversions of very large matrices in hardware, for example, FPGAs.

3. Enumerating $GL(n,2)$ Groups

In this section, the objective is to exhaustively enumerate the unique elements of a linear group, $GL(n,2)$. A naive approach may again exploit the PLU factorization, which has been assumed in the previous section for inverting very large binary matrices. The caveat is how to avoid a considerable redundancy by removing the duplicated matrices. Ideally, every unique invertible matrix of a given size is generated exactly once. The proposed solution is to consider the Bruhat decomposition theorem to define unique factorization of every invertible matrix [7]. This theorem can be specifically restated for general linear groups as follows.

Theorem 2 

(Bruhat decomposition of general linear groups). The general linear group, $GL(n,q)$, is a disjoint union of double cosets of the subgroup of upper triangular invertible matrices, $\mathit{B}$, i.e.,

$$GL(n,q)=\bigcup _{\mathit{P}\in {\mathit{S}}_n}\mathit{B}\mathit{P}\mathit{B}$$

(9)

where ${\mathit{S}}_n$ is the symmetric group of $(n!)$ permutation matrices.

Theorem 2 allows partitioning the group, $GL(n,2)$, into the disjoint subgroups (called Bruhat cells), which are indexed by the permutation matrices. More importantly, every invertible matrix, $\mathit{A}$, has a unique factorization, $\mathit{A}=\mathit{L}\mathit{P}\mathit{U}$, where $\mathit{P}$ is a permutation matrix, and $\mathit{L}$ and $\mathit{U}$ are the lower and upper triangular matrices, respectively, having the non-zero elements that are strictly constrained by the structure of $\mathit{P}$. The uniqueness of this decomposition enables the effective non-redundant enumeration of all invertible matrices.

In particular, the matrix, $\mathit{L}$, is defined by so-called inversions of the permutation, $\pi$, of a given permutation matrix, $\mathit{P}$. The inversion is an index pair, $(i,j)$, such that, $i<j$, whereas $\pi \left(i\right)>\pi \left(j\right)$. The matrix, $\mathit{L}$, is only allowed to have non-zero entries at positions, $(i,j)$. The number of inversions of permutation, $\pi$, is denoted as, $k=\mathrm{inv}\left(\pi \right)$. Thus, there are k elements of matrix, $\mathit{L}$, that can be freely selected from a given field, whereas other elements must be zero.

There are $(n!)$ disjoint cells corresponding to one permutation matrix, $\mathit{P}$. The order of each cell is a function of the number of inversions in $\mathit{P}$. For the field of cardinality, q, the cell order is given by the formula,

$$|C_{\mathit{P}}|=\prod _{k=0}^{n-1}{(q-1)}^n{q}^{{\displaystyle \frac{n(n-1)}{2}}}{q}^k\stackrel{q=2}{=}\prod _{k=0}^{n-1}{2}^{n(n-1)/2+k}.$$

(10)

The term, $2^{n(n-1)/2}$, accounts for the number of unconstrained upper-triangular matrices, $\mathit{U}$, and the term, $2^k$, accounts for the number of constrained lower-triangular matrix, $\mathit{L}$, respectively.

The enumeration algorithm proceeds by systematically enumerating all permutation matrices, $\mathit{P}$. For each matrix, $\mathit{P}$, all invertible matrices in the disjoint cell, $C_{\mathit{P}}$, are generated. Assuming that the two matrix multiplications to construct a matrix, $\mathit{A}$, have the complexity, $\mathcal{O}\left(n^3\right)$, and noting that the double cosets are disjoint and cover the whole GL group, which has the order given by (1), the overall complexity is,

$$\mathcal{O}\left(n^3\right)\sum _{\mathit{P}\in {\mathit{S}}_n}|C_{\mathit{P}}|=\mathcal{O}\left(n^3\right)|GL(n,2)|=\mathcal{O}\left(n^3\prod _{k=0}^{n-1}(2^n-2^k)\right)=\mathcal{O}\left(n^3{2}^{n^2}\right).$$

(11)

For comparison, the naive enumeration of a $GL(2,n)$ group using the PLU factorization is performed by creating all combinations of products between $n!$ permutation matrices, and $2^{n(n-1)}$ lower and upper triangular matrices, respectively. The overall complexity is dominated by the number of these products, each having, $\mathcal{O}\left(n^3\right)$, complexity. In addition, each resulting matrix must be checked whether it has already been generated previously, which has the complexity, $\mathcal{O}\left(n^2\right)$. Consequently, the total complexity of the naive enumeration via the PLU factorization is, $\mathcal{O}\left(n!2^{n(n-1)/2}{n}^3\right)=\mathcal{O}\left(n!n^3{2}^{n^2}\right)$. Consequently, even though the complexity of the Bruhat-based enumeration is still exponential, it avoids the factorial scaling of the PLU based enumeration, and at the same time, it is guaranteed to systematically produce a sequence of the unique invertible matrices.

3.1. Implementation

The rows of binary matrices can be stored in std::vector<int>, or in more efficient std::bitset<N> data structure to enable the compiled C++ code to use the low-level SIMD instructions. The Bruhat-based enumeration of the general linear group, $GL(2,n)$, was implemented in C++ [8]. For complexity comparison, and despite its inefficiency, the PLU-based enumeration was also implemented, but validated only for $n\le 5$. For smaller values of n, the uniqueness of the generated matrices can be automatically enforced by exploiting the std::set data structure with hashing or tree-based indexing. For larger n, more advanced data structures are required such as the hash-tries. The complete sets of all lower and upper triangular matrices can be pre-generated, which results in the large memory footprint, but faster generation. Alternatively, these matrices can be generated sequentially.

It should be emphasized that the exhaustive enumeration of general linear groups is only feasible for invertible matrices over finite fields, otherwise the group order becomes infinite. Moreover, deciding whether the two real-valued matrices are identical is problematic, as it depends on the floating-point representation of their values. The PLU-based enumeration can be effective, provided that it is modified, so that only unique matrices are implicitly or explicitly generated. Moreover, the PLU-based enumeration is already “embarrassingly parallel”, allowing the distributed matrix generation.

4. Repairing Singular Matrices to Become Invertible

In many scenarios, large binary matrices may be nearly invertible, however, they are strictly singular by definition. Such a case can be practically resolved by finding a `closest” invertible matrix that is non-singular. The problem can be formulated as follows.

Problem 1. 

Given a $(n\times n)$ binary matrix, $\mathit{A}$, find the binary matrix, ${\mathit{A}}^{\prime }\in GL(n,2)$, having the minimal Hamming distance, ${\mathcal{D}}_{\mathrm{H}}\left(\mathit{A},{\mathit{A}}^{\prime }\right)$.

Solution 1. 

$${\mathit{A}}^{\prime }=\left\{\begin{array}{cc}\mathit{A},& \mathit{A}\in GL(n,2)\\ \underset{{\mathit{A}}^{\prime }\in {\mathcal{Z}}_2^n}{arg\; min}{\mathcal{D}}_{\mathrm{H}}\left(\mathit{A},{\mathit{A}}^{\prime }\right),& \mathrm{otherwise}.\end{array}\right.$$

(12)

Specifically, in the first step, the rank deficiency, $d=n-rank\mathit{A}$, is computed to decide whether $\mathit{A}$ is invertible, or not. If $\mathit{A}$ is already invertible, i.e., $d=0$, the algorithm terminates. Otherwise, the invertible matrix, ${\mathit{A}}^{\prime }$, having the Hamming distance precisely, ${\mathcal{D}}_{\mathrm{H}}\left(\mathit{A},{\mathit{A}}^{\prime }\right)=d$, is obtained by flipping exactly d bits in $\mathit{A}$ as explained below.

The rank (i.e., also the rank deficiency) can be computed by the standard GJE. The matrix rank is equal to the number of pivots [4]. Specifically, when transforming the matrix into a row-echelon form, the j-th column does not have a pivot, if there is no row, $i\ge j$, such that the element, $\mathit{A}(i,j)=1$. The columns without a pivot are referred to as being pivot-free. Equivalently, the row-echelon form contains exactly d all-zero rows. It is also important to record, which rows are swapped in the GJE process. The rows that were swapped and zeroed are referred as dependent rows, since they are linear combinations of the other rows. The identified pivot-free columns and dependent rows are then used to define the bit positions that need to be flipped in order to make the original matrix invertible.

Lemma 2 

(Rank change due to single bit-flip). A single bit-flip changes the rank of binary matrix at most by one.

Proof. 

Let ${\mathit{A}}^{\prime }$ be a square binary matrix obtained by flipping the $(i,j)$ entry of $\mathit{A}$, so that, ${\mathit{A}}^{\prime }=\mathit{A}\oplus {\mathbf{E}}_{ij}$, and the matrix, ${\mathbf{E}}_{ij}$, has exactly one non-zero element at position, $(i,j)$. The rank of sum, $rank\left(\mathit{A}\oplus {\mathbf{E}}_{ij}\right)=rank\left(\mathit{A}+{\mathbf{E}}_{ij}\right)\le rank\mathit{A}+rank{\mathbf{E}}_{ij}=1+rank\mathit{A}$, since $rank{\mathbf{E}}_{ij}=1$. Similarly, $rank\left({\mathit{A}}^{\prime }\oplus {\mathbf{E}}_{ij}\right)\le 1+rank{\mathit{A}}^{\prime }$, and thus, $rank{\mathit{A}}^{\prime }\ge rank\mathit{A}-1$. Consequently, $|rank\mathit{A}-rank{\mathit{A}}^{\prime }|\le 1$.    □

Assuming Lemma 2, the task now is to find the set of d bit-flips that reduce the rank deficiency of $\mathit{A}$ to zero, so that it becomes invertible. These bit-flips are defined by the following theorem.

Theorem 3 

(Minimal bit-flips to make binary matrix invertible). The binary square matrix, $\mathit{A}$, with rank deficiency, $d>0$, can be made invertible by performing no less than and exactly d bit-flips at the intersection of the dependent rows and the pivot-free columns.

Proof. 

By Lemma 2, every bit-flip can increase the matrix rank by at most one. Thus, the minimum number of bit-flips required to make the matrix invertible is precisely equal to its rank deficiency, d. The GJE produces the row-echelon matrix,

$$\mathrm{GJE}\left(\mathit{A}\right)=\left(\begin{array}{cc}\mathit{B}& \mathbf{C}\\ {\mathbf{0}}_1& {\mathbf{0}}_2\end{array}\right)$$

(13)

where $\mathit{B}$ is a $(r\times r)$ invertible submatrix, $r=rank\mathit{A}$, and ${\mathbf{0}}_2$ is a $(d\times d)$ all-zero matrix. The d all-zero rows are the swapped dependent rows, and the last d columns are pivot-free [4]. Since every bit-flip at the intersection of dependent row and pivot-free column is guaranteed to increase the rank exactly by one, it is necessary as well as sufficient to define the d unique bit-flips. This is achieved by replacing the matrix, ${\mathbf{0}}_2$, with a permutation matrix, ${\mathit{P}}_d$, of the same size, i.e., the modified matrix, ${\mathit{A}}^{\prime }$, becomes,

$${\mathit{A}}^{\prime }=\left(\begin{array}{cc}\mathit{B}& \mathbf{C}\\ {\mathbf{0}}_1& {\mathit{P}}_d\end{array}\right).$$

(14)

For instance, ${\mathit{P}}_d={\mathit{I}}_d$, is one possible solution. The rank of ${\mathit{A}}^{\prime }$ is the sum of ranks of its diagonal blocks, i.e., $rank{\mathit{A}}^{\prime }=rank\mathit{B}+rank{\mathit{P}}_d=r+d=n$, so it is invertible. Furthermore, the modified matrix, ${\mathit{A}}^{\prime }$, has the minimal Hamming distance, ${\mathcal{D}}_{\mathrm{H}}\left({\mathit{A}}^{\prime },\mathit{A}\right)=d$, which is equal to the rank deficiency of $\mathit{A}$.    □

4.1. Implementation

The C++ implementation again uses std::bitset<N> data structure to store the matrix rows to allow the compiler to utilize the SIMD instructions for bit-level parallelism on a CPU [8]. For instance, a binary matrix containing $3,{000}^2$ bits requires approximately $1.13$ MB of memory. The time complexity of making the matrix invertible is dominated by the GJE, which has the complexity, $\mathcal{O}\left(n^3\right)$. The GJE computations can be at least partially distributed, which is important in high-throughput applications requiring to invert large binary matrices such as designing forward error correction coding, cryptography, and quantized machine learning algorithms. Moreover, the GJE arithmetic operations for the matrices defined over $\mathrm{GF}\left(2\right)$ field are readily amenable to fast hardware implementations including, for example, FPGAs.

The empirical runtimes and the memory requirements of repairing randomly generated binary matrices to become invertible are shown in Figure 3 and Figure 4, respectively. Note that a binary random matrix is invertible with the probability, $|GL(n,2)|/2^{n^2}\doteq 28.9\%$, for $n\gg 1$.

It is unclear whether the required bit-flips to make the binary matrix invertible could be found more directly in order to reduce the cubic cost of the GJE. The complexity of the GJE itself can be reduced by adopting faster matrix multiplications, for example, using the Method of Four Russians, and other Strassen-like algorithms [2].

5. Conclusions

The paper investigated the problem of inverting very large binary matrices. Specifically, a robust blockwise matrix inversion was proposed that is based on the PLU matrix factorization. The time complexity of this algorithm was evaluated both theoretically and empirically. Next, it was shown that the enumeration of general linear groups can be done much more efficiently by assuming the Bruhat decomposition rather than the PLU decomposition in order to only generate the unique invertible matrices. Finally, it was shown how to find an invertible binary matrix that is closest to a given singular binary matrix by performing the minimum number of bit-flips at the determined positions. The number of these bit-flips is equal to the rank deficiency of the original matrix. All proposed algorithms were validated to show their correctness. The reference C++ implementation is available from the public Github repository.