Essays/Cholesky Decomposition
< Essays
Jump to navigation
Jump to search
Given a positive definite matrix A , the Cholesky decomposition computes a lower triangular matrix L such that
mp=: +/ . * NB. matrix product h =: +@|: NB. conjugate transpose A -: L mp h L
L is a sort of "square root" of the matrix A .
A recursive solution for the Cholesky decomposition obtains by considering A as a 2-by-2 matrix of matrices and substituting appropriate matrix operations for scalar ones in the Cholesky method for a 2-by-2 matrix of scalars.
Cholesky=: 3 : 0 n=.#A=.y if. 1>:n do. assert. (A=|A)>0=A NB. check for positive definite %:A else. 'X Y t Z'=. , (;~n$(>.-:n){.1) <;.1 A L0=.Cholesky X L1=.Cholesky Z-(T=.(h Y) mp %.X) mp Y L0,(T mp L0),.L1 end. )
The verb is now applied to a few examples:
A =: _4]\ 33 7j_8 7j_10 3j_4, 7j8 28 2j4 _10j_11, 7j10 2j_4 22 3j3, 3j4 _10j11 3j_3 16 A 33 7j_8 7j_10 3j_4 7j8 28 2j4 _10j_11 7j10 2j_4 22 3j3 3j4 _10j11 3j_3 16 L=: Cholesky A L 5.74456 0 0 0 1.21854j1.39262 4.95739 0 0 1.21854j1.74078 _0.3851j_0.892453 4.06695j_2.72987e_17 0 0.522233j0.696311 _2.34116j2.19446 0.543008j_0.00121275 2.15659j1.02056e_16 0~:L 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1 *./, (>:/~i.#L) >: 0~:L NB. L is lower triangular 1 A -: L mp h L 1 ] B=: 3 3$ 1 4 6 4 0 3 6 3 2 1 4 6 4 0 3 6 3 2 Cholesky B NB. B is not positive definite |assertion failure: Cholesky | (A=|A)>0=A A1=: (+/ .* |:) _10+10 20 ?.@$ 20 $A1 10 10 L1=: Cholesky A1 A1 -: L1 mp h L1 1 *./, (>:/~i.#L1) >: 0~:L1 1
See also
- Cholesky Decomposition
- LU Decomposition
- QR Decomposition
- Matrix Inverse
- Triangular Matrix Inverse
- Determinant
- Minors
- Hilbert Matrix
- Block Matrix Inverse
- Kronecker Product
Contributed by Roger Hui. Substantially the same text previously appeared as part of A Note on Programming Style in Vector, Volume 12, Number 3, January 1996.