Cholesky Decomposition

There are a bunch of different Cholesky decomposition algorithms in the text.

RandomSym[n_]:=Module[
{Q,Evals},
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]];
Evals= RandomReal[{0.1,12},n];
Q.DiagonalMatrix[Evals].Transpose[Q]
]

n=25;
A=RandomSym[n];
MatrixPlot[A]

Inner

CholeskyInnerOne[A_]:=Module[
{n=Length[A],R=0*A,i,j, Tol =(10^(-15)),ShouldBePositive},
If[Norm[A-Transpose[A]] >Norm[A] Tol, Return[Print["CholeskyInnerOne: Not Symm!"]]];
Do[
ShouldBePositive=A[[i,i]]-Sum[((R[[k,i]])^2),{k,1,i-1}];
If[ShouldBePositive <0, Return[Print["CholeskyInnerOne: Not +Def!"]]];
R[[i,i]]=(Sqrt[ShouldBePositive]);
Do[
R[[i,j]]=(A[[i,j]]-Sum[R[[k,i]] R[[k,j]],{k,1,i-1}])/R[[i,i]],
{j,i+1,n}],
{i,1,n}];
R
]

Lets test the thing!

n=4;
A=RandomSym[n];
A[[3,3]]=-12.3;
L=CholeskyInnerOne[A]
A-Transpose[L].L

Data=Table[
{n,A=RandomSym[n];Timing[CholeskyInnerOne[A];][[1]]},
{n,20, 320,5}]

CubicDemo[n_]:=(10^(-6))(n^3);
ListLogLogPlot[
{Data,{{20,CubicDemo[20]}, {320,CubicDemo[320]}}},
PlotRange->All,
Frame->True,
GridLines->Automatic,
Joined->True,
Mesh->All
]

Lets see if we can make this faster.

CholeskyInnerTwo[A_]:=Module[
{n=Length[A],R=0*A,i,j,ri,rj,rii},
Do[
ri=R[[1;;i-1,i]];
rii=(Sqrt[A[[i,i]]-ri.ri]);
R[[i,i]]=rii;
Do[
rj=R[[1;;i-1,j]];
R[[i,j]]=(A[[i,j]]-ri.rj)/rii,
{j,i+1,n}],
{i,1,n}];
R
]

CholeskyInnerThree[A_]:=Module[
{n=Length[A],R=0*A,i,j,ri,rj,rii},
Do[
ri=R[[1;;i-1,i]];
rii=(Sqrt[A[[i,i]]-ri.ri]);
R[[i,i]]=rii;
rj=R[[1;;i-1,i+1;;n]];
R[[i,i+1;;n]]=(A[[i,i+1;;n]]-ri.rj)/rii,
{i,1,n}];
R
]

Like magic these go faster. Almost 9 times faster and then another factor of 7 or so.

n=123;A=RandomSym[n];
{Timing[R1=CholeskyInnerOne[A];][[1]],
Timing[R2=CholeskyInnerTwo[A];][[1]],
Timing[R3=CholeskyInnerThree[A];][[1]]}

n=5;A=RandomSym[n];
CholeskyInnerOne[A]
CholeskyInnerThree[A]

They are still cubic

DataInnerTwo=Table[
{n,A=RandomSym[n];Timing[CholeskyInnerTwo[A];][[1]]},
{n,20, 420,5}]
DataInnerThree=Table[
{n,A=RandomSym[n];Timing[CholeskyInnerThree[A];][[1]]},
{n,20, 420,5}]

CubicDemo[n_]:=(10^(-6))(n^3);
ListLogLogPlot[
{DataInnerTwo,DataInnerThree,
{{20,CubicDemo[20]}, {420,CubicDemo[420]}}
},
PlotRange->All,
Frame->True,
GridLines->Automatic,
Joined->True,
Mesh->All
]

Outer

On page 42 there is another version of the Cholesky decomposition.

CholeskyOuterOne[A_]:=Module[
{n=Length[A]},
AA=A;R=0*A;
Do[
R[[i,i]]=(Sqrt[AA[[i,i]]]);
s=AA[[i,i+1;;-1]]/R[[i,i]];
R[[i,i+1;;-1]]=s;
AA[[i+1;;-1,i+1;;-1]]=AA[[i+1;;-1,i+1;;-1]]-Outer[Times, s,s],
{i,1,n}];
R
]

CholeskyOuterTwo[A_]:=Module[
{n=Length[A]},
AA=A;R=0*A;
Do[
R[[i,i]]=(Sqrt[AA[[i,i]]]);
s=AA[[i,i+1;;-1]]/R[[i,i]];
R[[i,i+1;;-1]]=s;
AA[[i+1;;-1,i+1;;-1]]-=Outer[Times, s,s],
{i,1,n}];
R
]

CholeskyOuterThree[A_]:=Module[
{n=Length[A],r11=(Sqrt[A[[1,1]]]),s, AHat},
s=A[[1,2;;-1]]/r11;
AHat=A[[2;;-1,2;;-1]]-Outer[Times,s,s];
Flatten[
If[Length[A]>1,
Join[Join[{r11},s],CholeskyOuterThree[AHat]],
{r11}
]
]
]

Note that the outout form of the last one is different.  All the numbers match but the output is different.

n=5;A=RandomSym[n];
CholeskyInnerOne[A]
CholeskyOuterOne[A]
CholeskyOuterTwo[A]
CholeskyOuterThree[A]

Lets compare our  versions

$RecursionLimit=(10^5);
n=250;
A=RandomSym[n];
{Timing[CholeskyInnerOne[A];][[1]],
Timing[CholeskyInnerTwo[A];][[1]],
Timing[CholeskyInnerThree[A];][[1]],
Timing[CholeskyOuterOne[A];][[1]],
Timing[CholeskyOuterTwo[A];][[1]],
Timing[CholeskyOuterThree[A];][[1]]
}

So the fastest is over 30 times faster than the slowest (inner one)!

DataOuterOne=Table[
{n,A=RandomSym[n];Timing[CholeskyOuterOne[A];][[1]]},
{n,100, 420,5}];
DataOuterTwo=Table[
{n,A=RandomSym[n];Timing[CholeskyOuterTwo[A];][[1]]},
{n,100, 420,5}];
DataOuterThree=Table[
{n,A=RandomSym[n];Timing[CholeskyOuterThree[A];][[1]]},
{n,100, 420,5}];

There is now a funny pattern.

CubicDemo[n_]:=(10^(-8))(n^3);
ListLogLogPlot[
{DataOuterOne,DataOuterTwo,DataOuterThree,
{{100,CubicDemo[100]}, {420,CubicDemo[420]}}
},
PlotRange->All,
Frame->True,
GridLines->Automatic,
Joined->True,
Mesh->All,
PlotStyle->{Red, Green, Blue, Cyan}
]

n=560;
DataOuterThree=Table[
A=RandomSym[n];Timing[CholeskyOuterTwo[A];][[1]],
{120}];
ListPlot[DataOuterThree,
PlotRange->All,
Joined->True,
Mesh->All]