14

3

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 `Nest`

ed `ReplacePart`

function, but I think it would be just some pretentious sweeping under the rug. Am I mistaken?

Thank you in advance!

3Welcome to the site! – Dr. belisarius – 2012-06-11T14:21:32.790

1In case you are not aware there is a built-in function

`CholeskyDecomposition`

. – Mr.Wizard – 2012-06-11T14:23:53.0233Thanks 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. :) – Editortoise-Composerpent – 2012-06-11T14:25:33.473

4Here's a functional solution:

`Nest[CholeskyDecomposition, matrix, 1]`

:P – rm -rf – 2012-06-11T14:43:16.220@R.M Superb anti-answer! :) – Dr. belisarius – 2012-06-11T15:03:56.857

@R.M that's just being cheeky. :) – rcollyer – 2012-06-11T15:10:40.140

You might get some useful ideas from the notebook available here

– Daniel Lichtblau – 2012-06-11T15:40:26.640I think one сould also consider LU decomposition as well for completeness. – faleichik – 2012-06-11T20:29:05.343

2As 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... :) – J. M.'s ennui – 2012-06-12T09:20:12.367I tried the original posting which seems to work well with a small set of data, then applied a 100 x 100 covariance matrix which worked on CholeskyDecomposition[covMat]), but applying it to the ProceduralCholesk[covMat] I get "Infinite expression 1/0.x10-4 encountered"? Not sure how to copy and paste my covMat data for you to compare if it matters? – sebastian c. – 2013-01-18T12:50:27.973

I guess it's because my code fails to work whenever a zero has been previously computed on the leading diagonal of the lower triangular matrix. I should have added an exception, but I was probably too lazy to do it at the time. :D I'm just wondering if there are any precision issues involved or if you are actually "supposed to" have a zero in that element, though. – Editortoise-Composerpent – 2013-01-21T18:30:52.177