In[]:=
deploy
Thu 14 Jul 2022 15:53:30
Also see https://www.wolframcloud.com/obj/yaroslavvb/newton/whitening_inverse_transposed.nb
In[]:=
X={{10,0,4},{0,10,-2},{4,-2,15}};​​Print["X=",X//MatrixForm];​​(*forX,findlowertriangularCandsetsglobal"diag"suchthatX=C.diag.diag.C*)​​mchol[cov_]:=Module[{},​​chol=CholeskyDecomposition[cov];​​diag=DiagonalMatrix@Diagonal@chol;​​chol.Inverse[diag]​​];​​​​C0=mchol[X];​​D0=diag;​​Print["C0=",C0//MatrixForm]​​Print["D0=",D0*D0//MatrixForm]​​Print["decomp X=",MatrixForm/@{C0,D0.D0,C0}]​​Print["numcheck: ",C0.D0.D0.C0==X]​​​​Print["
-1
X
=",Inverse@X//MatrixForm];​​T=Inverse[C0];​​Print["decomp inv(X)=",MatrixForm/@{Inverse[X],T,Inverse[D0.D0],T}];​​Print["numcheck2: ",T.Inverse[D0.D0].T==Inverse@X]​​​​C1=mchol@Inverse@X;​​D1=diag;​​Print["T=",T//MatrixForm];​​Print["C1=",C1//MatrixForm];​​Print["D1=",D1*D1//MatrixForm];​​Print["decomp2 inv(X)=",MatrixForm/@{Inverse@X,C1,D1.D1,C1}]​​Print["numcheck3: ",C1.D1.D1.C1==Inverse@X]​​​​
X=
10
0
4
0
10
-2
4
-2
15
C0=
1
0
0
0
1
0
2
5
-
1
5
1
D0=
10
0
0
0
10
0
0
0
13
decomp X=
1
0
0
0
1
0
2
5
-
1
5
1
,
10
0
0
0
10
0
0
0
13
,
1
0
2
5
0
1
-
1
5
0
0
1

numcheck: True
-1
X
=
73
650
-
2
325
-
2
65
-
2
325
67
650
1
65
-
2
65
1
65
1
13
decomp inv(X)=
73
650
-
2
325
-
2
65
-
2
325
67
650
1
65
-
2
65
1
65
1
13
,
1
0
-
2
5
0
1
1
5
0
0
1
,
1
10
0
0
0
1
10
0
0
0
1
13
,
1
0
0
0
1
0
-
2
5
1
5
1

numcheck2: True
T=
1
0
0
0
1
0
-
2
5
1
5
1
C1=
1
0
0
-
4
73
1
0
-
20
73
2
15
1
D1=
73
650
0
0
0
15
146
0
0
0
1
15
decomp2 inv(X)=
73
650
-
2
325
-
2
65
-
2
325
67
650
1
65
-
2
65
1
65
1
13
,
1
0
0
-
4
73
1
0
-
20
73
2
15
1
,
73
650
0
0
0
15
146
0
0
0
1
15
,
1
-
4
73
-
20
73
0
1
2
15
0
0
1

numcheck3: True