Given a positive definite matrix A , the Choleski decomposition computes a lower triangular matrix L such that
mp=: +/ . * NB. matrix product h =: +@|: NB. conjugate transpose A -: L mp h L
A recursive solution reveals itself by considering A as a 2-by-2 matrix of matrices and substituting appropriate matrix operations for scalar ones in the Choleski method for a 2-by-2 matrix of scalars.
Choleski=: 3 : 0
n=.#A=.y
if. 1>:n do.
assert. (A=|A)>0=A NB. check for positive definite
%:A
else.
p=.>.n%2 [ q=.<.n%2
X=.(p,p){.A [ Y=.(p,-q){.A [ Z=.(-q,q){.A
L0=.Choleski X
L1=.Choleski 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=: Choleski 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
Choleski B NB. B is not positive definite
|assertion failure: Choleski
| (A=|A)>0=A
A1=: (+/ .* |:) _10+10 20 ?.@$ 20
$A1
10 10
L1=: Choleski A1
A1 -: L1 mp h L1
1
*./, (>:/~i.#L1) >: 0~:L1
1
See also
Contributed by RogerHui. Substantially the same text previously appeared as part of A Note on Programming Style in Vector, Volume 12, Number 3, January 1996.