What's the most "functional" way to do Cholesky decomposition?
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
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}]