# Cholesky/LDLt that works for singular matrices?

TLDR When given $$A=X'X$$ decomposition, user298737 works by relying on QRDecomposition on $$X$$. If we didn't have this decomposition, one could use QR to get full-rank part of $$A$$ and use Cholesky on that, and then project back to original space like in this answer. However, this requires full QR decomposition whereas Mathematica only provides reduced QR decomposition.

Cholesky decomposition fails for matrices that are singular in machine precision. Earlier question suggested to use LDTLt in ResourceFunction["RationalCholeskyDecomposition"], which is more robust than CholeskyDecomposition, but still fails if the matrix is too close to singular.

Any tips what I can do in Mathematica to deal with effectively singular matrices?

LDLt decomposition of $$(X^T X)^{-1}$$ has a statistical interpretation if we treat $$X$$ as data matrix of rows $$x$$ :

• $$i$$th row of $$L$$ give coefficients of optimal predictive model that predicts $$i$$'th entry of $$x$$ from previous entries
• entries of $$D$$ gives residual variances of models
• For singular $$X^T X$$, bottom rows of $$L$$ have a lot of freedom, problem is ill-posed

Example below computes total residual variance using Cholesky/LDLt, but this fails for numerically singular matrices.

SeedRandom;
dataSize = 30;
dimSize = 30;
sigma = DiagonalMatrix@Table[Exp[-i^2], {i, 1, dimSize}];
X = Normalize /@
RandomVariate[MultinormalDistribution[sigma], dataSize];
gram = X\[Transpose] . X;
chol2 = ResourceFunction["RationalCholeskyDecomposition"];

(* Total residual variance of best model predicting x[i] from x[j], j<i *)
Total@Last@chol2@gram


Fails with Divide::indet: Indeterminate expression 0./0. encountered.

• Is efficiency important? Jul 27, 2022 at 20:44
• Ideally it would be faster that reducing it to solving $n$ least-squares problems (I'll post that solution in a bit) Jul 27, 2022 at 21:30
• Singular Value Decomposition (SVD) is exactly for isolating the (numerically) singular part of a matrix from the non-singular part. Jul 28, 2022 at 5:24
• Perhaps Householdertransformation QRDecomposition  helps? Jul 28, 2022 at 6:27

If I understand correctly, given a real square matrix $$X$$ you are asking for an upper triangular real matrix $$R$$ such that $$X^T X = R^T R$$. One could use a QR decomposition $$X = QR$$, since then $$X^TX = (QR)^T QR = R^T (Q^TQ) R = R^T R$$ If we also want the diagonal entries to be nonnegative, then we need a little adjustment:

CholeskyDecompositionOfXTX[X_]:=With[{R=Last[QRDecomposition[X]]},
Map[If[Negative[#],-1,1]&,Diagonal[R]]*R];


Example.

(* example of an almost singular matrix *)
X = With[{n=40}, Array[N[Range[0,1,1/n]]^#&,n]];

(* calculate Cholesky and check *)
R = CholeskyDecompositionOfXTX[X];
UpperTriangularMatrixQ[R] (* true *)
Norm[Transpose[R].R-Transpose[X].X] (* small *)

• Nice! That works. Also Ben Grossman posted a recipe for how to do this with QR when you don't have $X^T X$ decomposition provided, although my first attempt to turn this into Mathematica didn't work (ie, QRDecomposition gives rectangular matrices, his recipe expected square) math.stackexchange.com/a/4502026/998 Jul 30, 2022 at 1:41

You can find a rational approximation to a specified accuracy.

With[{ϵ = 10^-12},
gramrat = Rationalize[gram, ϵ];
cholrat =
CholeskyDecomposition[
gramrat + ϵ IdentityMatrix[Dimensions[gram]]];
];

gram - Transpose[cholrat] . cholrat // Abs // Max
(* 1.4726*10^-12 *)


I'm not sure it's a good idea, but it might be worth investigating.

In NCAlgebra we have a code for calculating the LDL decomposition that will handle singular matrices using a pivoting strategy. See documentation for details. Don't expect it to be efficient though. Also you can use LDLDecomposition instead of NCLDLDecomposition if you only have numerical entries.

Below is an approach (hackedLdl) which takes $$X$$ and produces LDLt decomposition of $$X^T X$$. Gives same results as ResourceFunction["RationalCholeskyDecomposition"][Transpose[X].X] most of the time, but also seems to work for nearly singular matrices. Unclear how to extend this to the case when $$A=X^T X$$ decomposition is not available.

X1 = {{-2, 0, 0}, {-1, -1, -2}, {-2, 0, 1}};
X2 = {{-2, 0, 0, 0}, {-1, -1, -2, -2}, {-2, 0, 1, 1}};

On[Assert];
(* Returns coefficients needed to predict i'th coordinate from \
previous coords *)
onlyFirstKColumns[data_, k_] := Module[{n},
n = Last@Dimensions@data;
Drop[data\[Transpose], -(n - k)]\[Transpose]
];
extractKColumn[data_, k_] := data[[All, {k}]];

phiL2R[i_, data_] := Module[{},
coefs = Switch[i,
1,
{},  (* empty list,
since no "previous coordinates" for the first coord *)
_,
(* Drop[data\[Transpose],-(n-i+1)]\[Transpose];  (* predictors,
drop columns > i *) *)
Y1 = onlyFirstKColumns[data, i - 1];
Y2 = extractKColumn[data, i];(* target variable, column i *)
p = LeastSquares[Y1\[Transpose] . Y1, (Y1\[Transpose] . Y2)];
residuals = Y2 - Y1 . p;
Flatten@p
];
(* pad column with zeros *)
];

hackedLdl[X_] := Module[{},
n = Length@First@X;
coefsLeftToRight = Table[phiL2R[i, X], {i, n}]\[Transpose];
l = Inverse[IdentityMatrix[n] - coefsLeftToRight];
R2 = X - X . coefsLeftToRight;
d = Total[#*#] & /@ Transpose[R2];
{l\[Transpose], DiagonalMatrix@d}
];
{l, d} = hackedLdl[X1];
Print["L=", l // MatrixForm];
Print["d=", d // MatrixForm];
Print["Equality test1: ",
l . d . l\[Transpose] == X1\[Transpose] . X1];

{l, d} = hackedLdl[X2];
Print["L=", l // MatrixForm];
Print["d=", d // MatrixForm];
Print["Equality test2: ",
l . d . l\[Transpose] == X2\[Transpose] . X2];

SeedRandom;
dataSize = 30;
dimSize = 30;
sigma = DiagonalMatrix@Table[Exp[-i^2], {i, 1, dimSize}];
X3 = Normalize /@
RandomVariate[MultinormalDistribution[sigma], dataSize];
gram = X3\[Transpose] . X3;
{l, d} = hackedLdl[X3];
Print["Equality test3: ",
Max[Abs[Flatten[l . d . l\[Transpose] - X3\[Transpose] . X3]]] <
10^-10];