Cholesky-like decomposition that works on singular matrices? Is there a variant of Cholesky-like decomposition that works for singular matrices?
My problem is that all implementations of Cholesky I found fail when the matrix is singular in machine precision.
One idea for such decomposition is to consider duality between Cholesky and Gram-Schmidt. Once you encounter vectors that are linear combinations of previous vectors, skip them. This for instance produces the following decomposition for rank-3 matrix below.
$$\left(
\begin{array}{ccccc}
 9 & 1 & 0 & 0 & 0 \\
 1 & 1 & 2 & 2 & 2 \\
 0 & 2 & 5 & 5 & 5 \\
 0 & 2 & 5 & 5 & 5 \\
 0 & 2 & 5 & 5 & 5 \\
\end{array}
\right)=\left(
\begin{array}{ccc}
 1 & 0 & 0 \\
 \frac{1}{9} & 1 & 0 \\
 0 & \frac{9}{4} & 1 \\
 0 & \frac{9}{4} & 1 \\
 0 & \frac{9}{4} & 1 \\
\end{array}
\right)\left(
\begin{array}{ccc}
 9 & 0 & 0 \\
 0 & \frac{8}{9} & 0 \\
 0 & 0 & \frac{1}{2} \\
\end{array}
\right)\left(
\begin{array}{ccccc}
 1 & \frac{1}{9} & 0 & 0 & 0 \\
 0 & 1 & \frac{9}{4} & \frac{9}{4} & \frac{9}{4} \\
 0 & 0 & 1 & 1 & 1 \\
\end{array}
\right)$$
notebook

An approach using $QR$ decomposition: the following results in a decomposition of the form $A = MM^T$ with $M$ of full rank of size $n \times r$ (with $r$ equal to the rank of $A$), but $M$ is not necessarily in lower-triangular form. If it is desired, $M$ can be put into lower-triangular form with a further $QR$ decomposition $M^T = QR$.
Begin with a (symmetric) positive semidefinite $n \times n$ matrix $A$ of rank $r<n$. An exact $QR$ decomposition $A = QR$ must have the form
$$
Q = \pmatrix{Q_1 & Q_2}, \quad R = \pmatrix{R_1\\ 0}
$$
where $Q$ is square of size $n$, $Q_1$ has $r$ columns and $R_1$ has $R$ rows. Thus, the columns of $Q_1$ form an orthonormal basis for the column-space of $A$.
It follows that $Q_1^TAQ_1$ is strictly positive definite. Moreover, we have $Q_2^TA = (AQ_2)^T = 0$. Compute a Cholesky decomposition $LL^T = Q_1^TAQ_1$. From there, we may deduce that
$$
\pmatrix{L \\0} \pmatrix{L\\0}^T = \pmatrix{LL^T & 0\\0 & 0} = \pmatrix{Q_1^TAQ_1 & Q_1^TAQ_2\\ Q_2^TAQ_1 & Q_2^TAQ_2} = Q^TAQ.
$$
Conclude that
$$
A = Q\pmatrix{L\\0}\pmatrix{L \\ 0}^TQ^T = \left[Q\pmatrix{L \\0} \right]\left[ Q\pmatrix{L \\0}\right]^T
$$
