**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](https://math.stackexchange.com/a/4502026/998). 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](https://mathematica.stackexchange.com/questions/255970/cholesky-like-rank-revealing-decomposition-on-singular-matrices/256531#256531) 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](https://stats.stackexchange.com/questions/579315/statistical-interpretation-of-diagonal-of-cholesky-decomposition) 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.`