(*:Version: Mathematica 2.0 *) (*:Name: LinearAlgebra`Cholesky` *) (*:Title: Symmetric Matrix Routines *) (*:Author: Jerry B. Keiper (Wolfram Research Inc.) Revised by Hon Wah Tam (Wolfram Research Inc.) *) (*:Keywords: symmetric matrices, Cholesky decomposition *) (*:Requirements: The input matrix is assumed to be symmetric positive definite. Cholesky.m does not test for symmetry. *) (*:Warnings: Does not check whether A is symmetric positive definite, or even complex. *) (*:Sources: *) (*:Summary: For a given symmetric positive definite matrix A, CholeskyDecomposition[A] returns an upper-triangular matrix U such that A = Transpose[U] . U. This function uses the row-oriented Cholesky decomposition method to compute U row by row. That A is symmetric positive definite stipulates that A is real. The answers are totally garbage if A is complex. If A is not positive definite, U will contain Indeterminate entries. *) BeginPackage["LinearAlgebra`Cholesky`"] CholeskyDecomposition::usage = "For a symmetric positive definite matrix A, CholeskyDecomposition[A] returns an upper-triangular matrix U such that A = Transpose[U] . U." Begin["LinearAlgebra`Cholesky`private`"] column[matrix_List, index_Integer]:= Map[(#[[index]])&, matrix] diagonalprod[m_List, u_List, index_]:= Module[{partial = Take[column[u, index], index-1]}, Sqrt[m[[index, index]]-partial . partial]] offdiagonalprod[m_List, u_List, partial1_List, index1_, index2_, diag_]:= Module[{ partial2 = Take[column[u, index2], index1-1]}, (m[[index1, index2]]-partial1 . partial2)/diag] makerow[m_List, u_List, rowindex_]:= Module[{part1 = Table[0, {rowindex-1}], partial1 = Take[column[u, rowindex], rowindex-1 ], diag = diagonalprod[m, u, rowindex], part2, iter}, part2 = Table[offdiagonalprod[m, u, partial1, rowindex, iter, diag], {iter, rowindex+1, Length[m]}]; Join[part1, {diag}, part2]] CholeskyDecomposition[m_List] := Module[{len = Length[m], u, iter1}, u = DiagonalMatrix[Table[0, {len}]]; For[iter1=1, iter1<= len, ++iter1, u[[iter1]] = makerow[m, u, iter1]]; u] End[] (* LinearAlgebra`Cholesky`Private` *) Protect[CholeskyDecomposition] EndPackage[] (* LinearAlgebra`Cholesky` *) (*:Limitations: Does not decompose complex matrices. *) (*:Examples: CholeskyDecomposition[ Table[ 1/(i+j),{i,1,5},{j,1,5}] ] *)