<<   >>

41. Cholesky Decomposition

If A is a Hermitian, positive-definite matrix, its Cholesky decomposition is a lower-triangular matrix L such that A ≡ L +.× +⍉L . The matrix L is a sort of “square root” of the matrix A .

For example:

   ⎕rl←7*5

   A←t+.×⍉t←¯10+?5 5⍴20
   A
231   42  ¯63  16  26
 42  199 ¯127 ¯68  53
¯63 ¯127  245  66 ¯59
 16  ¯68   66 112 ¯75
 26   53  ¯59 ¯75  75

   L←Cholesky A
   L
15.1987   0        0        0       0
 2.7634  13.8334   0        0       0
¯4.1451  ¯8.35263 12.5719   0       0
 1.05272 ¯5.12592  2.1913   8.93392 0
 1.71067  3.48957 ¯1.81055 ¯6.15028 4.33502

   A ≡ L +.× +⍉L
1

For real matrices, “Hermitian” reduces to symmetric and the conjugate transpose +⍉ to transpose ⍉ . The symmetry arises in solving least-squares problems.

The Algorithm [111]

A recursive solution for the Cholesky decomposition obtains by considering A as a 2-by-2 matrix of matrices. It is algorithmically interesting but not necessarily the best with respect to numerical stability.

Cholesky←{
  1≥n←≢⍵:⍵*0.5
  p←⌈n÷2
  q←⌊n÷2
  X←(p,p)↑⍵ ⊣ Y←(p,-q)↑⍵ ⊣ Z←(-q,q)↑⍵
  L0←∇ X
  L1←∇ Z-T+.×Y ⊣ T←(+⍉Y)+.×⌹X
  ((p,n)↑L0)⍪(T+.×L0),L1
}

The recursive block matrix technique can be used for triangular matrix inversion [112], LU decomposition [113], and QR decomposition [114].

Proof of Correctness

The algorithm can be stated as a block matrix equation:

  ┌───┬───┐          ┌──────────────┬──────────────┐
  │ X │ Y │          │   L0 ← ∇ X   │       0      │
∇ ├───┼───┤  ←→  L ← ├──────────────┼──────────────┤
  │+⍉Y│ Z │          │    T+.×L0    │L1 ← ∇ Z-T+.×Y│
  └───┴───┘          └──────────────┴──────────────┘

where T←(+⍉Y)+.×⌹X . To verify that the result is correct, we need to show that A≡L+.×+⍉L and that L is lower triangular. For the first, we need to show:

┌───┬───┐     ┌──────┬───────┐     ┌────────┬────────┐
│ X │ Y │     │  L0  │   0   │     │  +⍉L0  │+⍉T+.×L0│
├───┼───┤  ≡  ├──────┼───────┤ +.× ├────────┼────────┤
│+⍉Y│ Z │     │T+.×L0│   L1  │     │    0   │  +⍉L1  │
└───┴───┘     └──────┴───────┘     └────────┴────────┘

that is:

(a)  X     ≡ L0 +.× +⍉L0
(b)  Y     ≡ L0 +.× +⍉ T+.×L0
(c)  (+⍉Y) ≡ (T+.×L0) +.× +⍉L0
(d)  Z     ≡ ((T+.×L0) +.× (+⍉T+.×L0)) + (L1+.×+⍉L1)

(a) holds because L0 is the Cholesky decomposition of X .

(b) is seen to be true as follows:

L0 +.× +⍉ T+.×L0
L0 +.× +⍉ ((+⍉Y)+.×⌹X)+.×L0  definition of T
L0 +.× (+⍉L0)+.×(+⍉⌹X)+.×Y +⍉A+.×B ←→ (+⍉B)+.×+⍉A and
+⍉+⍉Y ←→ Y
(L0+.×+⍉L0)+.×(+⍉⌹X)+.×Y +.× is associative
X+.×(+⍉⌹X)+.×Y (a)
X+.×(⌹X)+.×Y X and hence ⌹X are Hermitian
I+.×Y associativity; matrix inverse
Y identity matrix

(c) follows from (b) by application of +⍉ to both sides of the equation.

(d) turns on that L1 is the Cholesky decomposition of Z-T+.×Y :

((T+.×L0)+.×(+⍉T+.×L0)) + (L1+.×+⍉L1)
((T+.×L0)+.×(+⍉T+.×L0)) + Z-T+.×Y
((T+.×L0)+.×(+⍉L0)+.×+⍉T) + Z-T+.×Y
(T+.×X+.×+⍉T) + Z-T+.×Y
(T+.×X+.×+⍉(+⍉Y)+.×⌹X) + Z-T+.×Y
(T+.×X+.×(+⍉⌹X)+.×Y) + Z-T+.×Y
(T+.×X+.×(⌹X)+.×Y) + Z-T+.×Y
(T+.×I+.×Y) + Z-T+.×Y
(T+.×Y) + Z-T+.×Y
Z

Finally, L is lower triangular if L0 and L1 are lower triangular, and they are by induction.

A Complex Example

   ⎕rl←7*5

   A←t+.×+⍉t←(¯10+?5 5⍴20)+0j1ׯ10+?5 5⍴20
   A
382        17J131  ¯91J¯124 ¯43J0107  20J0035
 17J¯131  314     ¯107J0005 ¯60J¯154  26J¯137
¯91J0124 ¯107J¯05  379       49J0034  20J0137
¯43J¯107  ¯60J154   49J¯034 272       35J0103
 20J¯035   26J137   20J¯137  35J¯103 324

   L←Cholesky A

   A ≡ L +.× +⍉L
1
   0≠L
1 0 0 0 0
1 1 0 0 0
1 1 1 0 0
1 1 1 1 0
1 1 1 1 1

A Personal Note

This way of computing the Cholesky decomposition was one of the topics of [115] and was the connection (through Professor Shlomo Moran) by which I acquired an Erdős number of 2.


Appeared in [116, 117].