# What's the most “functional” way to do Cholesky decomposition?

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] :=
(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!

-
Welcome to the site! –  belisarius Jun 11 '12 at 14:21
In case you are not aware there is a built-in function CholeskyDecomposition. –  Mr.Wizard Jun 11 '12 at 14:23
Thanks for welcoming me! Anyway, I am aware of the existence of the built-in function and I have used it to test results for both implementations, but for my Computational Physics course I have to implement a Cholesky decomposition by myself. I am free to do it in whatever programming language I want, so I picked Mathematica since it is the one I am most familiar with. I know that the procedural implementation would be more than enough for my course, but I'm trying to learn functional programming as well by myself in the meanwhile. :) –  Andrea Colonna Jun 11 '12 at 14:25
Here's a functional solution: Nest[CholeskyDecomposition, matrix, 1] :P –  The Toad Jun 11 '12 at 14:43
As a tiny note, PositiveDefiniteMatrixQ[] internally computes a Cholesky decomposition to prove the positive-definiteness of a matrix, so in effect, you're doing a Cholesky decomposition twice... :) –  Guess who it is. Jun 12 '12 at 9:20

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

-
It indeed looks better, but wouldn't Append make it actually run slower than the procedural approach? –  Andrea Colonna Jun 11 '12 at 14:59
@Andrea you're not a beginning user, are you? :-) Possibly, but I am working toward elegance, not performance. If your ultimate goal is performance you need to include that in your question. –  Mr.Wizard Jun 11 '12 at 15:06
@Andrea how large are the matrices you will be using? –  Mr.Wizard Jun 11 '12 at 15:08
An alternate to Append is a linked list form. It may speed things up, while retaining a measure of elegance. –  rcollyer Jun 11 '12 at 15:11
@rcollyer I'm not sure that will be true in this case as the intermediate result needs to be used; it cannot be flattened (only) at the end. –  Mr.Wizard Jun 11 '12 at 15:13

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}}


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

-
Sometimes loops are the right choice. Despite everyone's best efforts I did not feel that any answer truly improved on the original. –  Mr.Wizard Jun 12 '12 at 14:06
Well, that question I asked on triangular recursions quite a while back also comes to mind. Still, I feel that outer product Cholesky is much more compact than OP's original route... –  Guess who it is. Jun 12 '12 at 14:09