Cholesky Decomposition

Morten was visiting Dyalog clients and forwarded a request: Can we have the 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:

⎕io←0 ⋄ ⎕rl←7*5 ⋄ ⎕pp←6

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.

Some writers asserted that an algorithm for the Cholesky decomposition “cannot be expressed without a loop”  and that “a Pascal program is a natural way of expressing the essentially iterative algorithm” . You can judge for yourself whether the algorithm presented here belies these assertions.

The Algorithm 

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←{
⍝ Cholesky decomposition of a Hermitian positive-definite matrix
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 , LU decomposition , and QR decomposition .

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

⎕io←0 ⋄ ⎕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  and was the connection (through Professor Shlomo Moran) by which I acquired an Erdős number of 2.

References

1. Wikipedia, Cholesky decomposition, 2014-11-25.
2. Thomson, Norman, J-ottings 7, The Education Vector, Volume 12, Number 2, 1995, pp. 21-25.
3. Muller, Antje, Tineke van Woudenberg, and Alister Young, Two Numerical Algorithms in J, The Education Vector, Volume 12, Number 2, 1995, pp. 26-30.
4. Hui, Roger, Cholesky Decomposition, J Wiki Essay, 2005-10-14.
5. Hui, Roger, Triangular Matrix Inverse, J Wiki Essay, 2005-10-27.
6. Hui, Roger, LU Decomposition, J Wiki Essay, 2005-10-31.
7. Hui, Roger, QR Decomposition, J Wiki Essay, 2005-10-30.
8. Ibarra, Oscar, Shlomo Moran, and Roger Hui, A Generalization of the Fast LUP Matrix Decomposition Algorithm and Applications, Journal of Algorithms 3, 1982, pp. 45-56.

5 thoughts on “Cholesky Decomposition”

1. Randy A MacDonald on said:

Showing Dyalog 14 features suggest this could be present: p q←(⌈,⌊)n÷2

2. Roger Hui on said:

You are correct. In this case I prefer the

p←⌈n÷2
q←⌊n÷2

formulation as I wanted to have the symmetry when the expressions are vertically aligned.

• Veli-Matti on said:

Hmm.. That v14/atop was the first idea that came to my mind, too, but could you as well do simple
(p q)←⌈1 ¯1×0.5×n
and then just
X←p p↑⍵⊣Y←p q↑⍵⊣Z←q q↑⍵
..and later
(p n↑L0)⍪(T+.×L0),L1

Yours -wm

3. Roger Hui on said:

0. As I said, I do want the expressions to line up vertically. My eyes can more readily see the symmetry in

p←⌈n÷2
q←⌊n÷2

and hence the relationship between p and q than in

(p q)←⌈1 ¯1×0.5×n

1. I prefer to let p and q have the same sign.

2. I have a personal antipathy against strand notation in most situations. I prefer (p,p)↑⍵ over p p↑⍵. I know some people would prefer the latter over the former.

4. Vered Madar on said:

Actually Roger’s method can be further developed.
There is a new paper I was submitting to Stat and Prob letters in which I uncover every entry of Cholesky decomposition.
Actually there are two forms one that uses semi-partial correlations and a second form that uses successive ratios of differences between sub-determinants.
See http://arxiv.org/abs/1412.1181v2 for axiv version.