Essays/Cholesky Decomposition

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

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.

Essays/Cholesky Decomposition

Navigation menu

Search