Cholesky/LDLt that works for singular matrices?

Question

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[1];
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.

Answer 1

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 *)
And@@Thread[Diagonal[R]>=0] (* true *)
Norm[Transpose[R].R-Transpose[X].X] (* small *)

Answer 2

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.

Answer 3

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.

Answer 4

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 *)
   PadRight[coefs, n]
   ];

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[1];
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];