*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.