What's the most "functional" way to do Cholesky decomposition?

Question

I can do Cholesky in a procedural style, such as:

ProceduralCholesky[matrin_List?PositiveDefiniteMatrixQ] := 
 Module[{dimens,
   ll},
  dimens = Length@matrin;
  ll = ConstantArray[0, {dimens, dimens}];
  Do[
   Do[
     If[i == j,
       ll[[i, i]] = 
         Sqrt[matrin[[i, i]] - 
           Sum[Conjugate[#]*# &@ll[[i, k]], {k, 1, i - 1}]];,
       ll[[i, 
          j]] = (matrin[[i, j]] - 
            Sum[ll[[i, k]]*Conjugate[ll[[j, k]]], {k, 1, j - 1}])/
          Conjugate[ll[[j, j]]];];,
     {i, j, dimens}];,
   {j, dimens}];
  ll
  ]

Moreover, I've seen a "half-functional" implementation, which, however, features a Table function and an outer For loop. So far I have managed to address the need for the Table function by writing:

HalfFunctionalCholesky[matrin_List?PositiveDefiniteMatrixQ] := 
 Module[{dimens,
   uu},
  dimens = Length@matrin;
  uu = ConstantArray[0, {dimens, dimens}];
  Do[uu[[i]] = makerow[matrin, i, uu, dimens], {i, dimens}];
  uu
  ]

makerow[matrin_List, rowindex_Integer, uu_, dimens_Integer] := 
 PadLeft[
  (Join[
      {#},
      offdiagonalelements[matrin, uu, rowindex, dimens, #]]) &@
   diagonalelement[matrin, uu, rowindex],
  dimens]

diagonalelement[matrin_, uu_, rowindex_] := 
 Sqrt[matrin[[rowindex, rowindex]] - Conjugate[#].# &@
   uu[[;; (rowindex- 1), rowindex]]]

offdiagonalelements[matrin_, uu_, rowindex_, dimens_, 
  diag_] := (matrin[[rowindex, (rowindex+ 1) ;; dimens]] -
    Conjugate[
      uu[[;; (rowindex - 1), 
       rowindex]]].uu[[;; (rowindex - 1), (rowindex + 1) ;; 
       dimens]])/diag

I am not satisfied yet. Can I avoid using the outer Do loop by turning it into a Fold or a Nest? The only idea I was able to come up with involved a Nested ReplacePart function, but I think it would be just some pretentious sweeping under the rug. Am I mistaken?

Thank you in advance!

Answer 1

For reference there is a built-in function CholeskyDecomposition.

For improving your existing code Array may be a minor subjective improvement:

HalfFunctionalCholesky2[matrin_List?PositiveDefiniteMatrixQ] := 
 Module[{dimens, uu},
  dimens = Length[matrin]; 
  uu = ConstantArray[0, {dimens, dimens}]; 
  Array[(uu[[#]] = makerow[matrin, #, uu, dimens]) &, dimens];
  uu
 ]

---

After a look through this function, following your instincts I believe Nest may be used:

HalfFunctionalCholesky3[matrin_List?PositiveDefiniteMatrixQ] :=
  Module[{len = Length[matrin], i = 1},
    Nest[Append[#, makerow[matrin, i++, #, len]] &, {}, len]
  ]

---

Upon further examination I believe that the change to Nest allows a further simplification as follows (within Part):

diagonalelement[matrin_, uu_, rowindex_] := 
 Sqrt[ matrin[[rowindex, rowindex]] - Conjugate[#].# & @ uu[[All, rowindex]] ]

offdiagonalelements[matrin_, uu_, rowindex_, len_, diag_] :=
 (matrin[[rowindex, rowindex - len ;;]] - 
    Conjugate[uu[[All, rowindex]]].uu[[All, rowindex - len ;;]]) / diag

Answer 2

This is not a functional implementation (as it stands, it's rather MATLAB-ish), but I'll leave this snippet around and hope somebody could make something purely functional out of this outer product form of Cholesky decomposition:

m = Array[Min, {4, 4}]; (* example matrix *)

Do[
  m[[k, k]] = b = Sqrt[a = m[[k, k]]];
  m[[k, k + 1 ;;]] = (v = m[[k, k + 1 ;;]])/b;
  m[[k + 1 ;;, k + 1 ;;]] -= Outer[Times, v, v]/a;
  , {k, Length[m]}];
UpperTriangularize[m]

{{1, 1, 1, 1}, {0, 1, 1, 1}, {0, 0, 1, 1}, {0, 0, 0, 1}}

---

(added 6/12/2012)

As expected, the method given above can be made purely functional, but to me the functional version looks a lot less elegant:

k = 0; n = Length[m];
MapAt[Sqrt,
 Nest[Function[m, ++k;
   Block[{a = m[[k, k]], v = m[[k, k + 1 ;;]], s},
    s = Sqrt[m[[k, k]]];
    ArrayFlatten[{{ReplacePart[m[[1 ;; k, 1 ;; k]], {k, k} -> s], 
       DiagonalMatrix[Append[ConstantArray[1, k - 1], 1/s]].m[[1 ;; k, k + 1 ;;]]},
                 {0, m[[k + 1 ;;, k + 1 ;;]] - Outer[Times, v, v]/a}}]]],
  m, n - 1],
 {n, n}]