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:
- $m$-by-$m$ permutation matrix $P$;
- $m$-by-$n$ lower trapezoidal matrix $\tilde L$ with ones on the diagonal; and,
- $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.
