Mathematica Numerics for developers

Cholesky decomposition

Stored Timings and Plots

Needs["Graphics`Graphics`"]

\(Ctimes\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n, True]; 
	
        t\  = \ time[cd\  = \(CholeskyDecomposition[A]\)[\([1]\)]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]\)

Ctimes  = 
       {{2, 0.00009704589843749999`, 8.881784197001252`*^-16}, {3, 
          0.00009887695312500006`, 1.7763568394002505`*^-15}, {4, 
          0.00010253906249999998`, 3.552713678800501`*^-15}, {6, 
          0.00011047363281250014`, 7.105427357601002`*^-15}, {8, 
          0.00013305664062499955`, 1.4210854715202004`*^-14}, {11, 
          0.0001660156250000008`, 2.842170943040401`*^-14}, {16, 
          0.00026611328125000083`, 5.684341886080802`*^-14}, {23, 
          0.0005029296875000006`, 2.2737367544323206`*^-13}, {32, 
          0.0010644531249999964`, 1.1368683772161603`*^-12}, {45, 
          0.002519531250000012`, 9.094947017729282`*^-13}, {64, 
          0.0064453124999999944`, 3.637978807091713`*^-12}, {91, 
          0.01781250000000001`, 1.4551915228366852`*^-11}, {128, 
          0.04874999999999985`, 3.2741809263825417`*^-11}, {181, 
          0.13125000000000053`, 4.729372449219227`*^-11}, {256, 
          0.3775000000000013`, 1.8917489796876907`*^-10}, {362, 
          1.0250000000000057`, 1.7462298274040222`*^-10}, {512, 
          3.0799999999999983`, 1.1059455573558807`*^-9}};

Cp  =  
      LogLogListPlot[#[[{1, 2}]]& /@ Ctimes, PlotJoined -> True, 
        DisplayFunction -> Identity];

times1  = 
       {{2, 0.00164062500000000666`, 0.`}, {3, 0.00318359374999999111`, 
          1.77635683940025046`*^-15}, {4, 0.00499999999999989341`, 
          3.55271367880050092`*^-15}, {6, 0.0109374999999998223`, 
          7.10542735760100185`*^-15}, {8, 0.0246875000000006394`, 
          1.42108547152020037`*^-14}, {11, 0.0390625`, 
          2.84217094304040074`*^-14}, {16, 0.0962499999999977262`, 
          5.68434188608080148`*^-14}, {23, 0.252499999999997726`, 
          1.13686837721616029`*^-13}, {32, 0.550000000000011368`, 
          2.27373675443232059`*^-13}, {45, 1.17999999999994997`, 
          4.54747350886464118`*^-13}, {64, 3.03999999999996362`, 
          1.81898940354585647`*^-12}, {91, 7.76999999999998181`, 
          1.81898940354585647`*^-12}, {128, 20.3100000000000591`, 
          7.2759576141834259`*^-12}, {181, 54.4400000000000527`, 
          1.09139364212751388`*^-11}, {256, 148.339999999999916`, 
          2.91038304567337036`*^-11}, {362, 419.850000000000012`, 
          4.36557456851005554`*^-11}, {512, 1166.99999999999981`, 
          1.16415321826934814`*^-10}};

p1  =  
      LogLogListPlot[#[[{1, 2}]]& /@ times1, 
        PlotStyle -> {RGBColor[0, 0, 0], PointSize[0.025]}, 
        DisplayFunction -> Identity];

times2  = 
       {{2, 0.00139648437499999978`, 0.`}, {3, 0.00242187499999999689`, 
          1.77635683940025046`*^-15}, {4, 0.00351562500000000133`, 
          3.55271367880050092`*^-15}, {6, 0.006249999999999992`, 
          7.10542735760100185`*^-15}, {8, 0.00882812500000001953`, 
          1.42108547152020037`*^-14}, {11, 0.0139843749999999932`, 
          2.84217094304040074`*^-14}, {16, 0.0254687499999999822`, 
          5.68434188608080148`*^-14}, {23, 0.0478124999999999289`, 
          2.27373675443232059`*^-13}, {32, 0.0912499999999998223`, 
          2.27373675443232059`*^-13}, {45, 0.188749999999999751`, 
          4.54747350886464118`*^-13}, {64, 0.447499999999999786`, 
          1.81898940354585647`*^-12}, {91, 1.10999999999999943`, 
          3.63797880709171295`*^-12}, {128, 2.89000000000000056`, 
          7.2759576141834259`*^-12}, {181, 8.46999999999999886`, 
          1.09139364212751388`*^-11}, {256, 23.270000000000004`, 
          2.91038304567337036`*^-11}, {362, 66.6600000000000125`, 
          4.36557456851005554`*^-11}, {512, 181.980000000000003`, 
          1.16415321826934814`*^-10}};

p2  =  
      LogLogListPlot[#[[{1, 2}]]& /@ times2, 
        PlotStyle -> {RGBColor[0, 1, 0], PointSize[0.025]}, 
        DisplayFunction -> Identity];

times3  = 
       {{2, 0.000280761718749911182`, 8.88178419700125232`*^-16}, {3, 
          0.000466308593750186517`, 1.77635683940025046`*^-15}, {4, 
          0.000742187499999769073`, 3.55271367880050092`*^-15}, {6, 
          0.00158203124999989341`, 7.10542735760100185`*^-15}, {8, 
          0.00300781249999992894`, 1.42108547152020037`*^-14}, {11, 
          0.00652343750000028421`, 2.84217094304040074`*^-14}, {16, 
          0.01953125`, 5.68434188608080148`*^-14}, {23, 
          0.0662500000000250111`, 1.13686837721616029`*^-13}, {32, 
          0.221249999999940882`, 2.27373675443232059`*^-13}, {45, 
          0.800000000000181898`, 4.54747350886464118`*^-13}, {64, 
          2.96999999999934516`, 9.09494701772928237`*^-13}, {91, 
          16.4200000000000727`, 1.81898940354585647`*^-12}, {128, 
          90.7499999999999928`, 3.63797880709171295`*^-12}};

p3  =  
      LogLogListPlot[#[[{1, 2}]]& /@ times3, 
        PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.025]}, 
        DisplayFunction -> Identity];

ptimes1  = 
       {{2, 0.0013183593750001332`, 8.881784197001252`*^-16}, {3, 
          0.002675781250000231`, 1.7763568394002505`*^-15}, {4, 
          0.004570312500000284`, 0.`}, {6, 0.010390625000001208`, 
          7.105427357601002`*^-15}, {8, 0.018906249999997016`, 
          1.4210854715202004`*^-14}, {11, 0.03812500000000085`, 
          2.842170943040401`*^-14}, {16, 0.08937500000000398`, 
          5.684341886080802`*^-14}, {23, 0.21250000000000568`, 
          1.1368683772161603`*^-13}, {32, 0.4824999999999591`, 
          2.2737367544323206`*^-13}, {45, 1.1600000000000819`, 
          4.547473508864641`*^-13}, {64, 2.9700000000000273`, 
          9.094947017729282`*^-13}, {91, 8.029999999999973`, 
          1.8189894035458565`*^-12}, {128, 21.04000000000019`, 
          3.637978807091713`*^-12}, {181, 56.649999999999864`, 
          7.275957614183426`*^-12}, {256, 154.84000000000015`, 
          1.4551915228366852`*^-11}, {362, 428.8299999999999`, 
          2.9103830456733704`*^-11}, {512, 1197.62`, 
          5.820766091346741`*^-11}};

pp1  =  
      LogLogListPlot[#[[{1, 2}]]& /@ ptimes1, 
        PlotStyle -> {RGBColor[ .25,  .25,  .25], PointSize[0.025]}, 
        DisplayFunction -> Identity];

ptimes2  = 
       {{2, 0.0013574218749998757`, 8.881784197001252`*^-16}, {3, 
          0.0022460937499997335`, 1.7763568394002505`*^-15}, {4, 
          0.0031249999999998224`, 3.552713678800501`*^-15}, {6, 
          0.0049218749999999645`, 7.105427357601002`*^-15}, {8, 
          0.006757812500000071`, 1.4210854715202004`*^-14}, {11, 
          0.009609375000000142`, 1.4210854715202004`*^-14}, {16, 
          0.014531250000000995`, 5.684341886080802`*^-14}, {23, 
          0.022187500000001137`, 1.1368683772161603`*^-13}, {32, 
          0.03531250000000341`, 2.2737367544323206`*^-13}, {45, 
          0.05187499999999545`, 4.547473508864641`*^-13}, {64, 
          0.08249999999999602`, 9.094947017729282`*^-13}, {91, 
          0.14124999999998522`, 1.8189894035458565`*^-12}, {128, 
          0.27000000000003865`, 3.637978807091713`*^-12}, {181, 
          0.5750000000000455`, 7.275957614183426`*^-12}, {256, 
          1.3900000000001`, 1.4551915228366852`*^-11}, {362, 
          3.709999999999809`, 2.9103830456733704`*^-11}, {512, 
          9.720000000000027`, 5.820766091346741`*^-11}}

pp2  =  
      LogLogListPlot[#[[{1, 2}]]& /@ ptimes2, 
        PlotStyle -> {RGBColor[0,  .5,  .5], PointSize[0.025]}, 
        DisplayFunction -> Identity];

ptimes3  = 
       {{2, 0.00013427734375001665`, 8.881784197001252`*^-16}, {3, 
          0.0002075195312499778`, 1.7763568394002505`*^-15}, {4, 
          0.0003027343750000022`, 3.552713678800501`*^-15}, {6, 
          0.00056640625000004`, 2.220446049250313`*^-16}, {8, 
          0.0009228515625000488`, 1.4210854715202004`*^-14}, {11, 
          0.001708984375`, 1.4210854715202004`*^-14}, {16, 
          0.0036328125000002487`, 5.684341886080802`*^-14}, {23, 
          0.007343750000000426`, 1.1368683772161603`*^-13}, {32, 
          0.015234375000000355`, 2.2737367544323206`*^-13}, {45, 
          0.03187499999999943`, 4.547473508864641`*^-13}, {64, 
          0.06937499999999375`, 9.094947017729282`*^-13}, {91, 
          0.16499999999999204`, 1.8189894035458565`*^-12}, {128, 
          0.37000000000000455`, 3.637978807091713`*^-12}, {181, 
          0.8699999999998909`, 7.275957614183426`*^-12}, {256, 
          2.4299999999998363`, 1.4551915228366852`*^-11}, {362, 
          6.230000000000018`, 2.9103830456733704`*^-11}, {512, 
          16.710000000000036`, 5.820766091346741`*^-11}}

pp3  =  
      LogLogListPlot[#[[{1, 2}]]& /@ ptimes3, 
        PlotStyle -> {RGBColor[ .5, 0, 1], PointSize[0.025]}, 
        DisplayFunction -> Identity];

ptimes4  = 
       {{2, 0.0007666015624999967`, 8.881784197001252`*^-16}, {3, 
          0.0013281249999999994`, 1.7763568394002505`*^-15}, {4, 
          0.0018457031250000144`, 3.552713678800501`*^-15}, {6, 
          0.002832031249999978`, 7.105427357601002`*^-15}, {8, 
          0.004140625000000009`, 1.4210854715202004`*^-14}, {11, 
          0.007656250000000031`, 1.4210854715202004`*^-14}, {16, 
          0.014140625000000018`, 5.684341886080802`*^-14}, {23, 
          0.020624999999999893`, 1.1368683772161603`*^-13}, {32, 
          0.024218750000000178`, 2.2737367544323206`*^-13}, {45, 
          0.04062499999999947`, 4.547473508864641`*^-13}, {64, 
          0.08000000000000007`, 9.094947017729282`*^-13}, {91, 
          0.19875000000000043`, 1.8189894035458565`*^-12}, {128, 
          0.6550000000000011`, 3.637978807091713`*^-12}, {181, 
          1.7700000000000102`, 7.275957614183426`*^-12}, {256, 
          4.949999999999989`, 1.4551915228366852`*^-11}, {362, 
          13.060000000000002`, 2.9103830456733704`*^-11}, {512, 
          36.79999999999998`, 5.820766091346741`*^-11}};

Deriving the decomposition

For a thorough explanation see "Matrix Computations", Gene Golub and Charles Van Loan, 3rd ed, Johns Hopkins, 1996

(A  =  
          Table[If[i  >=  j, (Subscript[a, i, j]), (Subscript[a, j, i])], {i, 1, 4}, {j, 1, 
              4}])  //  MatrixForm

(L  =   
          Table[If[i  >=  j, (Subscript[l, i, j]), 0], {i, 1, 4}, {j, 1, 4}])  // 
       MatrixForm

(prod  =  L . Transpose[L])  //  MatrixForm

\(rules\  = \ {l\_\(1, 1\) -> \@a\_\(1, 1\)}; 

    prod\  /. \ rules\  // MatrixForm\)

\(rules\  = \ 
      Join[rules, {l\_\(j_, 1\ \) :> a\_\(j, 1\)\/l\_\(1, 1\)}]; 
    
prod\  //. rules\  // \ MatrixForm\)

\(rules\  = \ 
      Join[rules, {
          l\_\(2, 2\ \) -> \@\(a\_\(2, 2\) - a\_\(2, 1\)\^2\/a\_\(1, 1\)\)}]; 
    
prod\  //. rules\  // \ MatrixForm\)

Implementation 1

To make coding more convenient, define a function that determines if the diagonal element is large enough to use.  A tolerance ϵ is used to avoid dividing by numbers that are too small.

\(setdiagonal[j_, ϵ_]\  := \ 
	
      If[ϵ\  === \ Symbolic\  || \ L[\([j, j]\)]\  > \ ϵ, 

        		L[\([j, j]\)]\  = \ \@L[\([j, j]\)], 
\ \ \ \  (*\ else\ *) 

        		Message[CholeskyDecomposition::posdef]; 
		Throw[$Failed], 

        		 (*\ inequality\ didn' t\ eval\ *) 
		
        Message[CholeskyDecomposition::neval, \ L[\([j, j]\)], ϵ]; 

        		Throw[$Failed]]\)

I'll explain the reason for this definition later

zero  =  0;

Here is a function that, given the matrix and a tolerance ϵ, will compute the Cholesky decomposition.  It makes a copy of the matrix and then does the computation in place in the matrix L.

\(CholeskyDecomposition1[a_, ϵ_]\  := \ 
	
      Block[{n\  = \ Length[a], L\  = \ a}, 
		setdiagonal[1, ϵ]; 
        
		Do[
			Do[L[\([j - 1, i]\)]\  = \ zero, {i, j, n}]; 
			
          Do[L[\([i, j - 1]\)]\  /= \ L[\([j - 1, j\  - 1]\)], {i, j, n}]; 

          			Do[
            L[\([i, j]\)]\  -= \ 
              Sum[L[\([i, k]\)]*L[\([j, k]\)], {k, 1, j - 1}], {i, j, n}]; 

          			setdiagonal[j, ϵ], 
			{j, 2, n}]; 
		L]\)

Verification

For exact input 0 tolerance is reasonable.

CholeskyDecomposition1[{{2, 1}, {1, 1}}, 0]

% . Transpose[%]

Because we have written this in Mathematica, the function is polymorphic for many types:  in particular, exact numeric data (as shown above), symbolic matrices, machine number matrices, and arbitrary precision matrices.  Some more examples

(A  =  
          Table[If[i  >=  j, (Subscript[a, i, j]), (Subscript[a, j, i])], {i, 1, 3}, {j, 1, 
              3}])  //  MatrixForm

In this case, the data is symbolic, and we will not be able to verify the positivity of the diagonal elements before taking the square root, so we give ϵ === Symbolic which causes setdiagonal to ignore this.

CholeskyDecomposition1[A, Symbolic]

Simplify[% . Transpose[%]]  //  MatrixForm

Here is a random symmetric matrix

\(A\  = \ Table[Random[Integer, {\(-9\), 9}], {4}, {4}]; 

    A\  += \ Transpose[A]; 
MatrixForm[A]\)

One (relatively fast) way to determine if the matrix is positive definite is to simply try to Cholesky decomposition.

CholeskyDecomposition1[A, 0]

This matrix will be diagonally dominant and thus positive definite

A  +=  100 IdentityMatrix[4]  //  MatrixForm

With an exact matrix, the decompostion will be exact

\(cd\  = \ CholeskyDecomposition1[A, 0]; 

    Simplify[A\  - \ cd . Transpose[cd]]\)

The problem with using exact quantities is that the computations will become increasingly complex for bigger matrices.  These problems can be reduced by using an L.D.Transpose[L] decomposition instead, but the complexity still becomes overwhelming for large matrices.  For most applications, numerical data is quite sufficient.  For example, with machine numbers, it is appropriate to use machine ϵ, and

\(cd\  = \ CholeskyDecomposition1[N[A], $MachineEpsilon]; 

    A\  - \ cd . Transpose[cd]\)

With arbitrary precision numbers, the tolerance should be associated with the precision of the numbers.  In analogy with machine ϵ, define

Nϵ[prec_]  :=  N[2^(1 - Floor[prec*Log[2, 10]])]

Nϵ[$MachinePrecision]  ===  $MachineEpsilon

Now we obtain the accuracy you would expect from higher precision input.

\(cd\  = \ CholeskyDecomposition1[N[A, 32], Nϵ[32]]; 

    A\  - \ cd . Transpose[cd]\)

So, with a relatively simple definition, we now have a function that computes a Cholesky decomposition for generic type.  The next question is how fast does it run?

Asymptotic timing

Theoretically, the complexity of the decomposition is \(TraditionalForm
      \`\(O(n\^3)\) . \ \ Lets\ do\ some\ timings\ and\ \ \(see . \)\)

\(makepdsmatrix[n_]\  := \ 
      Module[{A}, 
		A\  = \ Table[Random[], {n}, {n}]; 
		
        A\  += \ Transpose[A]; 
		A\  += \ n^2\ IdentityMatrix[n]; 
		A]
        \)

This makes sure we don't run into the granularity of the operating system's timing facilities

\(SetAttributes[time, \ HoldAll]; 

    time[expr_]\  := \ 
      Block[{Second\  = \ 1, t, \ tries\  = \ 1}, \ 
\ \ \ \ 
        t\  = \ \(Timing[expr; ]\)[\([1]\)]; \ \ \ \ 
\ \ \ \ 
        While[t\  < \ 1. , tries\  *= \ 2; 
          t\  = \ \(Timing[Do[expr, {tries}]; ]\)[\([1]\)]]; 
\ \ \ \ 
        t/tries]; \)

If you want to check asymptotic complexity, often a good way to do it is to try the function over a large range of magnitudes.  Here we use a loop which increases the system size exponentially.  This has the advantage that on a log-log plot, the data will be equally spaced.

\(times1\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition1[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]\)

These timings were done a a Pentium PRO 200, which is a tad faster than this machine -- for comparison,

\(timeshere\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition1[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 3,  .5}]\)

Needs["Graphics`Graphics`"]

p1  =  
      LogLogListPlot[#[[{1, 2}]]& /@ times1, 
        PlotStyle -> PointSize[0.025]];

\(Clear[n]; 

    Fit[Log[Take[\(#[\([{1, 2}]\)]&\)\ /@\ times1, \(-4\)]], {1, ξ}, 
      ξ]\)

The problem is that this is getting quite slow as n gets large.  Just for comparison, we compare timings with an implementation of the Cholesky decomposition written in C.

Show[p1, Cp, DisplayFunction -> $DisplayFunction, PlotRange -> All];

\(Clear[n]; 

    Fit[Log[\(#[\([{1, 2}]\)]&\)\ /@\ Take[Ctimes, \(-8\)]], {1, ξ}, 
      ξ]\)

The asymptotic complexity is the same O((n^3))  The difference in coefficient, however is huge:

times1[[-1, 2]]/Ctimes[[-1, 2]]

Implementation 2

Before we blame Mathematica as being "too slow", lets work on the program a bit.

Asymptotically, the majority of the time is going to be spent on the innermost loop, which occurs in the Sum command. An examination of this indicates that this command is really just implementing a dot product.  In fact, if we look to the next loop out, we can see that it is really implementing a (sub)matrix-vector product.  But since Mathematica has a command to do this which is optimized for machine numbers, it seems sensible to use it.

\(CholeskyDecomposition2[a_, \ ϵ_]\  := 
	
      Block[{L, \ \ n\  = \ Length[a], rc, rr}, 
		L\  = \ a; \ 
		
        setdiagonal[1, ϵ]; 
\ \ \ \ 
        Do[
			rc\  = \ Range[j, \ n]; \ 
			
          rr\  = \ Range[1, j - 1]; 
			
          L[\([rc, j - 1]\)]\  /= \ L[\([j - 1, j - 1]\)]; 
			
          L[\([rc, j]\)]\  -= \ L[\([rc, rr]\)] . L[\([j, rr]\)]; 
			
          setdiagonal[j, ϵ]; \ 
			L[\([rr, j]\)]\  = \ zero, \ 
          
			{j, \ 2, \ n}]; 
		\ L]\)

We do a similar set of timing for this routine...

\(\(times2\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition2[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]; \)\)

The plot shows that we have made a significant improvement.

Show[p1, p2, Cp, DisplayFunction -> $DisplayFunction, 
      PlotRange -> All];

Just to be sure, we can check the complexity

\(Clear[n]; 

    fit2\  = \ 
      Fit[Log[Take[\(#[\([{1, 2}]\)]&\)\ /@\ times2, \(-4\)]], {1, ξ}, 
        ξ]\)

It looks better, but we are still not even close!

{times1[[-1, 2]]/times2[[-1, 2]], 
      times2[[-1, 2]]/Ctimes[[-1, 2]]}

The Compiler

The next option we have is to use the Mathematica compiler.  Because Compile is a HoldAll function, we need to change out definition somewhat -- we have included the setdiagonal functionality inline, and have modified the If in that to return consistent type.

\(CholeskyDecomposition3\  = \ 
      Compile[{{a, \ _Real, \ 2}, {ϵ, \ _Real}}, 
		
        Block[{L\  = \ a, \ \ n\  = \ Length[a], rc\  = \ Range[1], 
            rr\  = \ Range[1], d}, 
			
          If[L[\([1, 1]\)] > \ ϵ, 
				
            \(d\  = \ \(L[\([1, 1]\)]\  = \ Sqrt[L[\([1, 1]\)]]\); \), 
            \(Return[{{0. }}]; \)]; 
			
          Do[	
				rr\  = \ Range[1, j - 1]; 	
				
            Do[	
					L[\([i, j - 1]\)]\  /= \ d; \ \ 
					
              L[\([i, j]\)]\  -= \ L[\([i, rr]\)] . L[\([j, rr]\)], 

              					{i, j, n}]; 
				
            If[L[\([j, j]\)] > \ ϵ, \ 
              \(d\  = \ \(L[\([j, j]\)]\  = \ Sqrt[L[\([j, j]\)]]\); \), 
              \(Return[{{0. }}]; \)]; 
				
            Do[L[\([i, j]\)]\  = \ 0. , \ {i, 1, j - 1}], 

            				{j, \ 2, \ n}]; 
			L]]\)

A similar timing ...

\(\(times3\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition3[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 7,  .5}]; \)\)

... shows that something is amiss!!

Show[p1, p2, p3, Cp, DisplayFunction -> $DisplayFunction, 
      PlotRange -> All];

A computation of the asymptotic complexity shows ...

\(Clear[n]; 

    fit3\  = \ 
      Fit[Log[Take[\(#[\([{1, 2}]\)]&\)\ /@\ times3, \(-4\)]], {1, ξ}, 
        ξ]\)

... that something is way wrong.  The reason why is subtle.

PackedArray objects.

A PackedArray object is a special internal implementation of a list essentially the same as what the Mathematica compiler uses.  Unless you really want to know, a PackedArray feels, smells, and tastes just like a list.  To the internals of Mathematica, however, they are quite different.  Just like in the Mathematica compiler, they must be rectangular tensors, and all of the same type amongst machine integers, real, or complex numbers.

Since the matrices we have been trying timings for satisfy this definition, this is a good candidate problem for PackedArrays.

First redefine the matrix generator we use

\(makepdsmatrix[n_, packed_:  False]\  := \ 
      Module[{A}, 
		A\  = \ Table[Random[], {n}, {n}]; 
		
        A\  += \ Transpose[A]; 
		A\  += \ n^2\ IdentityMatrix[n]; 
		
        If[packed, Developer`Pack[A], Developer`Unpack[A]]]\)

There is a function to check if a list is stored internally as a PackedArray or not:

Developer`PackedArrayQ[makepdsmatrix[3, False]]

Developer`PackedArrayQ[A  =  makepdsmatrix[3, True]]

It still looks just like a list

A

and we can give it to our command

\(setdiagonal[j_, ϵ_]\  := \ 
	
      If[ϵ\  === \ Symbolic\  || \ L[\([j, j]\)]\  > \ ϵ, 

        		L[\([j, j]\)]\  = \ Sqrt[L[\([j, j]\)]], 

        \ \ \ \  (*\ else\ *) 
		Message[CholeskyDecomposition::posdef]; 

        		Throw[$Failed], 
		 (*\ inequality\ didn' t\ eval\ *) 
		
        Message[CholeskyDecomposition::neval, \ L[\([j, j]\)], ϵ]; 

        		Throw[$Failed]]\)

zero  =  0;

\(CholeskyDecomposition1[a_, ϵ_]\  := \ 
	
      Block[{n\  = \ Length[a], L\  = \ a}, 
		setdiagonal[1, ϵ]; 
        
		Do[
			Do[L[\([j - 1, i]\)]\  = \ zero, {i, j, n}]; 
			
          Do[L[\([i, j - 1]\)]\  /= \ L[\([j - 1, j\  - 1]\)], {i, j, n}]; 

          			Do[
            L[\([i, j]\)]\  -= \ 
              Sum[L[\([i, k]\)]*L[\([j, k]\)], {k, 1, j - 1}], {i, j, n}]; 

          			setdiagonal[j, ϵ], 
			{j, 2, n}]; 
		L]\)

CholeskyDecomposition1[A, $MachineEpsilon]

Developer`PackedArrayQ[%]

To get this to work, we need to redefine

zero  =  0. ;

So the whole matrix consists of floating point numbers

CholeskyDecomposition1[A, $MachineEpsilon]

Developer`PackedArrayQ[%]

So what does this do to the timings?

\(\(ptimes1\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n, True]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition1[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]; \)\)

\(CholeskyDecomposition2[a_, \ ϵ_]\  := 
	
      Block[{L, \ \ n\  = \ Length[a], rc, rr}, 
		L\  = \ a; \ 
		
        setdiagonal[1, ϵ]; 
\ \ \ \ 
        Do[
			rc\  = \ Range[j, \ n]; \ 
			
          rr\  = \ Range[1, j - 1]; 
			
          L[\([rc, j - 1]\)]\  /= \ L[\([j - 1, j - 1]\)]; 
			
          L[\([rc, j]\)]\  -= \ L[\([rc, rr]\)] . L[\([j, rr]\)]; 
			
          setdiagonal[j, ϵ]; \ 
			L[\([rr, j]\)]\  = \ zero, \ 
          
			{j, \ 2, \ n}]; 
		\ L]\)

\(\(ptimes2\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n, True]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition2[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]; \)\)

\(\(CholeskyDecomposition3\  = \ 
      Compile[{{a, \ _Real, \ 2}, {ϵ, \ _Real}}, 
		
        Block[{L\  = \ a, \ \ n\  = \ Length[a], rc\  = \ Range[1], 
            rr\  = \ Range[1], d}, 
			
          If[L[\([1, 1]\)] > \ ϵ, 
				
            \(d\  = \ \(L[\([1, 1]\)]\  = \ Sqrt[L[\([1, 1]\)]]\);\), 
            \(Return[{{0. }}]; \)]; 
			
          Do[	
				rr\  = \ Range[1, j - 1]; 	
				
            Do[	
					L[\([i, j - 1]\)]\  /= \ d; \ \ 
					
              L[\([i, j]\)]\  -= \ L[\([i, rr]\)] . L[\([j, rr]\)], 

              					{i, j, n}]; 
				
            If[L[\([j, j]\)] > \ ϵ, \ 
              \(d\  = \ \(L[\([j, j]\)]\  = \ Sqrt[L[\([j, j]\)]]\); \), 
              \(Return[{{0. }}]; \)]; 
				
            Do[L[\([i, j]\)]\  = \ 0. , \ {i, 1, j - 1}], 

            				{j, \ 2, \ n}]; 
			L]]; \)\)

\(\(ptimes3\  = \ 
      Table[
	n\  = \ Round[2^j]; 
	A\  = \ makepdsmatrix[n, True]; 
	
        t\  = \ time[cd\  = CholeskyDecomposition3[A, $MachineEpsilon]]; 
	
        err\  = \ Max[Abs[A\  - \ cd . Transpose[cd]]]; 
	{n, t, err}, 

        	{j, 1, 9,  .5}]; \)\)

Show[p1, p2, p3, pp1, pp2, pp3, Cp, 
      DisplayFunction -> $DisplayFunction, PlotRange -> All];

The second version with packed arrays is almost as fast as compiled C code.

{times2[[-1, 2]]/ptimes2[[-1, 2]], 
      ptimes2[[-1, 2]]/Ctimes[[-1, 2]]}

Improving the user interface

There is really no reason that a user should have to worry about the value of the tolerance variable ϵ used to ensure there are no division problems.  Also, as you can see there are big advantages to having the matrix be represented as a PackedArray.  It is also a good idea to check that the matrix is symmetric, too.

CholeskyDecomposition::symm  =  
      "The matrix `1` at position `2` is not symmetric";

CholeskyDecomposition::neval  =  
      "The diagonal element `1` could not be compared with the tolerance \
`2`.  You may need to numericize the matrix.";

CholeskyDecomposition::posdef  =  
      "The matrix is not sufficiently positive definite for a numerically \
accurate computation of the Cholesky decomposition.";

MyCholeskyDecomposition[a_]  :=  
      With[{cd  =  Catch[MyCholesky[a]]}, cd  /;  cd  =!=  $Failed];

\(MyCholesky[a_]\  := \ 
      Module[{A\  = \ a, dims, \ prec, \ machine, \ symbolic, \ ϵ}, 
        
		dims\  = \ Dimensions[A]; 
		
        If[Length[dims]\  < \ 2\  || \ Not[MatrixQ[A]], 
			
          Message[CholeskyDecomposition::matrix, A, 1]; 
		
          Throw[$Failed]]; 
		
        If[dims[\([1]\)]\  != \ dims[\([2]\)], \ 
				
          Message[CholeskyDecomposition::matsq, A, 1]; 	Throw[$Failed]]; 

        		If[\(! SameQ[A, Transpose[A]]\), \ 
				
          Message[CholeskyDecomposition::symm, A, 1]; 	Throw[$Failed]]; 

        		prec\  = \ Infinity; 
		machine\  = \ False; 
		
        symbolic\  = \ False; 
		
        ftest\ \  = \ 
          Function[{x}, 
            If[MachineNumberQ[x], machine\  = \ True; \ True, 
				
              If[Not[NumericQ[x]], \ symbolic\  = \ True; \ False, 

                						
                If[x\  != \ 0, prec\  = \ Min[Precision[x], prec], 
                  prec\  = \ Min[Accuracy[x], prec]]; \ True]]]; 
		
        MatrixQ[A, ftest]; 
		
        If[machine, 
			A\  = \ N[A]; 
			
          ϵ\  = \ 
            If[MatrixQ[A, NumberQ], $MachineEpsilon, Symbolic], 

          			 (*\ else\ *) 
			
          If[symbolic, 
				ϵ\  = \ Symbolic, 

            				 (*\ else\ *) 
				A\  = \ N[A, prec]; 
            ϵ\  = \ 2^\((1\  - \ Floor[Log[2. , 10]*prec])\)]]; 

        		Block[{zero}, If[machine, zero\  = \ 0. , \ zero\  = \ 0]; 
          CholeskyDecomposition2[A, ϵ]]]\)

MyCholeskyDecomposition[{{1, 2}, {3, 4}}]

MyCholeskyDecomposition[{{1. , 1}, {1, 1}}]

MyCholeskyDecomposition[{{2, 1}, {1, 1}}]

MyCholeskyDecomposition[{{2. , 1}, {1, 1}}]

MyCholeskyDecomposition[{{2`20, 1}, {1, 1}}]

So with this definition, concerns about precision and numeric/symbolic data are handled automatically.