How can one construct a sparse null space basis using recursive LU decomposition?

On the LUQ decomposition

The algorithm implemented in luq (see reference given below) computes bases for the left/right null spaces of a sparse matrix $A$. Unfortunately, as far as I can tell, there seems to be no thorough discussion of this particular algorithm in the literature. In place of a reference, let us clarify how/why it works and test it a bit.

The luq routine inputs an $m$-by-$n$ matrix $A$ and outputs an $m$-by-$m$ invertible matrix $L$, an $n$-by-$n$ invertible matrix $Q$ , and an $m$-by-$n$ upper trapezoidal matrix $U$ such that: (i) $A=LUQ$ and (ii) the pivot-less columns/rows of $U$ are zero vectors. For example, $$ \underbrace{\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}}_A = \underbrace{\begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}}_L \underbrace{\begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}}_U \underbrace{\begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}}_Q $$

Point (ii) allows one to construct bases for the left/right null spaces of $A$.

Bases for Left/Right Null Spaces of $A$

Let $r = \operatorname{Rank}(A)$. Suppose we can compute the exact $LUQ$ decomposition of $A$ as described above. Then,

• The $n-r$ columns of $Q^{-1}$ corresponding to the pivotless columns of $U$ are a basis for the null space of $A$. This follows from the fact that $\operatorname{null}(A) = \operatorname{null}(A Q^{-1}) = \operatorname{null}(L U)$ and that the pivotless columns of $U$ are zero vectors by construction.
• The $m-r$ rows of $L^{-1}$ corresponding to the pivotless rows of $U$ are a basis for the left null space of $A$. This follows from the fact that $\operatorname{null}(A^T) = \operatorname{null}((L^{-1} A)^T) = \operatorname{null}( (U Q)^T)$ and that the pivotless rows of $U$ are zero vectors by construction.

LUQ Algorithm

Assume that $m \ge n$. (If $m < n$, then the lu command mentioned below outputs a slightly different $PA=LU$ factorization. Otherwise the LUQ decomposition is almost the same, and so, we omit this case.)

Given an $m$-by-$n$ matrix $A$, the LUQ decomposition calls MATLAB command lu with partial (i.e., just row) pivoting. lu implements a variant of the LU decomposition that inputs $A$ and outputs:

1. $m$-by-$m$ permutation matrix $P$;
2. $m$-by-$n$ lower trapezoidal matrix $\tilde L$ with ones on the diagonal; and,
3. $n$-by-$n$ upper triangular matrix $\tilde U$

such that $PA = \tilde L \tilde U$. Write: $$ \tilde U = \begin{bmatrix} \tilde U_{11} & \tilde U_{12} \\ 0 & \tilde U_{22} \end{bmatrix} $$ where $\tilde U_{11}$ has nonzero diagonal entries, and hence, is invertible. Also, let $e_i$ denote unit $m$-vectors equal to $1$ in the $i$th component and zero otherwise. The algorithm then builds: $$ L = P^T \begin{bmatrix} \tilde L & e_{n+1} & \cdots & e_m \end{bmatrix} $$ which is an $m \times m$ invertible matrix, and $$ U = \begin{bmatrix} \tilde U_{11} & 0 \\ 0 & \tilde U_{22} \\ 0 & 0 \end{bmatrix} $$ which is upper trapezoidal, and $$ Q = \begin{bmatrix} I & \tilde U_{11}^{-1} \tilde U_{12} \\ 0 & I \end{bmatrix} $$ which is an $n$-by-$n$ invertible matrix. To summarize, we obtain: $$ A = L \begin{bmatrix} \tilde U_{11} & 0 \\ 0 & \tilde U_{22} \\ 0 & 0 \end{bmatrix} Q $$ For the most part, that is all the algorithm does. However, if there are any nonzero entries in $\tilde U_{22}$, then the algorithm will call luq again with input matrix containing all of the nonzero entries of $\tilde U_{22}$. This last step introduces more zeros into $U$ and modifies the invertible matrices $L$ and $Q$.

To understand this last step, it helps to consider a simple input to luq like $$ A = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} $$ The first call to luq with this input trivially gives $U=A$ with $L$ and $Q$ being the $3$-by-$3$ identity matrices. Since $U$ has nonzero entries, a second call is made to luq with input $1$, which outputs $L=U=Q=1$. This second decomposition is incorporated into the first one by making the second column of $L$ the first one and moving all the other columns to the right of it, and similarly, moving the third row of $Q$ to the first row and moving all the other rows below it. This yields, $$ A = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{pmatrix} $$

To be sure, consider another simple example $$ A = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & a & 0 & b \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & c \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$ where $a,b,c$ are nonzero reals. In the first pass through luq the algorithm again sets $U=A$ and $L$, $Q$ equal to the $5$-by-$5$ identity matrices. Since $U=\tilde U_{22}$ has nonzero elements, luq is called again with input matrix $$ B = \begin{pmatrix} a & b \\ 0 & c \end{pmatrix} $$ This is incorporated into the first decomposition by permuting $L$ and $Q$ as shown: $$ A = \begin{pmatrix} 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} a & b & 0 & 0 & 0 \\ 0 & c & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix} $$ In general, the columns of $L$ and the rows of $Q$ are permuted so that the the zero columns/rows of $\tilde U_{22}$ are moved to the end of the matrix. An LUQ decomposition is then performed on this nonzero sub-block.

A full explanation would be notation heavy (requiring index sets for the zero/nonzero elements) and not much easier to understand than the code itself.

Simple Test

In reality, the algorithm computes an approximate LUQ decomposition and approximate bases, i.e., with rounding errors. These rounding errors might be significant if some of the nonzero singular values of $A$ are too small for the algorithm to detect.

Here is a MATLAB script file that tests the luq code. The script is a slight modification of the demo file that the software comes with. I modified the original file so that it inputs a sparse, random, rectangular, rank deficient matrix and outputs bases for the left/right null spaces of this input matrix.

Here is a sample output from this demo file.

 elapsed time = 0.011993 seconds
Input matrix:
 size = 10000x500
 true right null space dimension = 23
 true left null space dimension = 9523
Output:
 estimated right null space dimension = 23
 estimated left null space dimension = 9523
 error in basis for right null space = 0
 error in basis for left null space = 2.2737e-13

"Extreme" Test

This example is adapted from Gotsman and Toledo [2008]. Consider the $(n+1)$-by-$n$ matrix: $$ A_1 = \begin{pmatrix} 1 & & & & \\ -1 & 1 & & & \\ \vdots & -1 & \ddots & & \\ \vdots & & \ddots & 1 & \\ -1 & -1 & \cdots & -1 & 1 \\ 0.5 & 0.5 & \cdots & 0.5 & 0.5 \end{pmatrix} $$ and in terms of this matrix, define the block diagonal matrix: $$ A = \begin{bmatrix} A_1 & 0 \\ 0 & A_2 \end{bmatrix} $$ where $A_2$ is an $n$-by-$n$ random symmetric positive definite matrix whose eigenvalues are all equal to one except $3$ are zero and one is $10^{-8}$. With this input matrix and $n=1000$, we obtain the following sample output.

 elapsed time = 1.1092
the matrix:
 size of A = 2001x2000
 true rank of A = 1997
 true right null space dimension = 3
 true left null space dimension = 4
results:
 estimated right null space dimension = 3
 estimated left null space dimension = 4
 error in basis for right null space = 9.2526e-13
 error in basis for left null space = 5.9577e-14

Remark

There is an option in the luq code to use LU factorization with complete (i.e., row and column) pivoting $PAQ=LU$. The resulting $U$ matrix in the $LUQ$ factorization may better reflect the rank of $A$ in more ill-conditioned problems, but there is an added cost to doing column pivoting.

Reference

Kowal, P. [2006]. "Null space of a sparse matrix."
https://www.mathworks.com/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix

Gotsman, C., and S. Toledo [2008]. "On the computation of null spaces of sparse rectangular matrices." SIAM Journal on Matrix Analysis and Applications, (30)2, 445-463.

Not the same algorithm, but here is an alternative that I have used in a recent paper; you can put it together quickly with standard Matlab functions if $A$ is a $m\times n$ matrix with $m\ll n$ (like in your example).

Compute a rank-revealing QR decomposition with [Q,R,E] = qr(A); now $E$ is a permutation matrix, and $A=QRE^T$. Set $$ R = \begin{bmatrix} R_1 & R_2\\ 0 & 0 \end{bmatrix}, $$ with $R_1$ square invertible and triangular. Sounds crazy when you think about it, but a matrix can be square and triangular at the same time. :)

Now $E \begin{bmatrix}-R_1^{-1}R_2\\ I \end{bmatrix}$ is a basis for the nullspace of $A$ with only $O(mn)$ nonzeros (and you can compute it in time $O(nm^2)$.

If $A$ is full rank, the second block row of zeros is going to be empty. You might need to set a threshold to decide which rows are zero and which are not in that decomposition; but if you want to compute a basis for the nullspace you have to take a numerical decision on its rank whatever you do.

At a first glance, it seems that the algorithm you linked to does something similar, but with a pivoted LU decomposition instead of a QR (which is indeed a better idea if $A$ is sparse in addition to being short and fat).

---

As for the name of your decomposition, if you don't impose any particular structure on $L,Q$ it could be everything from an SVD to a rank-revealing QR or LU. Probably, though, it is some variant of "rank-revealing LU decomposition"; this term is a good literature pointer to search for alternatives.