# 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 [0] 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” [1] and that “a Pascal program is a natural way of expressing the essentially iterative algorithm” [2]. You can judge for yourself whether the algorithm presented here belies these assertions.

## The Algorithm [3]

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 [4], LU decomposition [5], and QR decomposition [6].

## 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 [7] 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. Showing Dyalog 14 features suggest this could be present: p q←(⌈,⌊)n÷2

2. 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.

• 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. 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. 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.