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