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](http://www.mathworks.com/matlabcentral/fileexchange/11120-null-space-of-a-sparse-matrix) 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.



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](http://www.crab.rutgers.edu/~nb361/mysoftware/demo2.m) 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.


