Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
3  
Welcome to the site! –  belisarius Jun 11 '12 at 14:21
1  
In case you are not aware there is a built-in function CholeskyDecomposition. –  Mr.Wizard Jun 11 '12 at 14:23
3  
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
4  
Here's a functional solution: Nest[CholeskyDecomposition, matrix, 1] :P –  rm -rf Jun 11 '12 at 14:43
2  
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... :) –  J. M. Jun 12 '12 at 9:20
show 6 more comments

2 Answers

up vote 7 down vote accepted

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
share|improve this answer
    
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
1  
@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
show 2 more comments

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}]
share|improve this answer
    
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... –  J. M. Jun 12 '12 at 14:09
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.