Skip to main content
added 87 characters in body
Source Link
Yaroslav Bulatov
  • 6.7k
  • 1
  • 21
  • 47

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.

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

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.

Source Link
Yaroslav Bulatov
  • 6.7k
  • 1
  • 21
  • 47

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