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

Given an $m$ by $n$ matrix $A$ I'm familiar with the standard method to compute a basis for the null space of $A$ by computing a QR factorization of $A^T$. If $A$ is large and sparse, we can use sparse-QR techniques with pivoting, but the resulting $Q$ matrix (and portion used for the null space basis) may be quite dense.

For simple running example, consider a $n \times 1$ row of ones $A = [1 ... 1]$. This matrix has all zero-average vectors in its null space. If I use (MATLAB's) QR factorization ([Q,R,E] = qr(A');) the resulting $Q$ matrix is a dense $n \times n$ matrix. In this case, we know a sparse basis for the null space exists:

$N = \left[\begin{array}{c} \begin{array} --1 & -1& \dots & -1\\ \end{array} \\ I \end{array}\right]$, where $I$ is the $n-1 \times n-1$ (sparse) identity matrix.

A preliminary question is, 1) how can one construct a(/the most) sparse null space basis given $A$?

But, actually I have found an algorithm's implementation online by Pawel Kowal* that computes sparse null spaces very well. Trying to trace through the code it seems to work by recursively applying an LU-decomposition with pivoting. However, I can't completely understand what it's doing and certainly don't understand why it's working. The comments and function names say it is computing an "LUQ" decomposition:

function [L,U,Q] = luq(A,do_pivot,tol)
%  PURPOSE: calculates the following decomposition
%
%       A = L |Ubar  0 | Q
%             |0     0 |
%
%       where Ubar is a square invertible matrix
%       and matrices L, Q are invertible.
%

Is this decomposition well known? Does it go by another name?

This algorithm seems to work very well. In the example above, this decomposition produces the "ideal" sparse basis for the null space $N$.

So my current question, 2) is there a corresponding academic paper describing this method for computing a sparse basis for the null space of a matrix via recursive LU decomposition?

*I've had no luck trying to contact Pawel Kowal for more information.

• I almost found a solution in "Strong rank revealing LU factorizations" [Miranian & Gu 2003]. The first equation supposes we have permutations $P_1$ and $P_2$ such that: $P_1A _2P_2 = \left[ \begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \end{array} \right]$ with $A_{11}$ an invertible $k \times k$ matrix where $k$ is the rank of $A$. Then the paper claims that $N = \left[\begin{array}{c}-A_{11}^{-1}A_{12}\\ I\end{array}\right]$ forms the right null space of $A$, i.e., $AN=0$. But by my reckoning, $AN= \left[\begin{array}{c}0\\ -A_{21}^{-1}A_{11}A_{12}+A_{22}\end{array} \right] \neq 0$. – Alec Jacobson Nov 6 '16 at 2:16
• Gu and Miranian's decomposition chooses permutations and block sizes to ensure that $S(A_1) = A_{22}-A_{21}A_{11}^{-1}A_{12}$ (which is the content of the second block -- there is a typo in your formula) has norm under a certain small threshold, so $AN$ is numerically zero. – Federico Poloni Nov 6 '16 at 9:33
• Re: typo, yes "$A_{21}^{-1}A_{11}$" should be $A_{21}A_{11}^{-1}$. Ah, I didn't read on enough to find that they're making the Schur complement vanish (they don't really mention that near the introduction of the null space basis, so it seemed like the claim was that this always applied). – Alec Jacobson Nov 6 '16 at 15:56

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 . 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. . "Null space of a sparse matrix."
https://www.mathworks.com/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix

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

• I guess the part that I'm most interested in is indeed the part left out in: "For the most part, that is all the algorithm does. However, if there are any nonzero rows or columns in Ũ22, then the algorithm will call luq again with input matrix given by these nonzero columns/rows of Ũ22. This last step introduces more zeros in U and modifies the invertible matrices L and Q". Does calling luq on U22 guarantee somehow that zeros move to the bottom corner? How is the luq decomposition of U22 incorporated? – Alec Jacobson Nov 7 '16 at 14:43
• Based on your feedback, I added more detail. – Nawaf Bou-Rabee Nov 7 '16 at 19:12

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.

• For the m<<n, this seems like a great alternative. Could you add the name of your recent paper or where I can read more about this? I'd still like to find the reference/understand this LU-based method since it's working very well for m<<n and even for m=n. – Alec Jacobson Nov 5 '16 at 21:16
• @AlecJacobson My paper is dx.doi.org/10.1002/nla.2018, but there are a few differences from this setup; in particular it is combined with another nonstandard of decomposition. I would suggest you to search for papers on rank-revealing LU decompositions; that is probably the more appropriate keyword here. – Federico Poloni Nov 5 '16 at 21:26
• Thanks. Using your qr decomposition of A (rather than A^T) is it also possible to get a feasible solution to A x = b, given some b? This way, one can generate an affine basis for the null space of Ax=b. – Alec Jacobson Nov 5 '16 at 21:45
• @AlecJacobson Sure, $QRE^Tx=b$ is equivalent to $Ry=Q^Tb$, $x=Ey$. If you set $Q^Tb=\begin{bmatrix}b_1 \\ b_2\end{bmatrix}$, then the first system is solvable iff $b_2=0$, and a special solution is $y=\begin{bmatrix}R_1^{-1}b_1\\0\end{bmatrix}$. – Federico Poloni Nov 5 '16 at 21:53
• @FedericoPoloni, poor engineer here, it is an interesting discussion, I came across. I did not truly understand how you shaped R into the format written above. More specifically, how can I compute null space of a sparse matrix with the following entries (1,2)>1 (2,3)>1 (3,4)>1 (1,7)>-1 (2,8)>-1 (3,9)>-1? – Umut Tabak Jun 8 '17 at 17:10