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

- Wikipedia,
*Cholesky decomposition*, 2014-11-25. - Thomson, Norman,
*J-ottings 7*, The Education Vector, Volume 12, Number 2, 1995, pp. 21-25. - Muller, Antje, Tineke van Woudenberg, and Alister Young,
*Two Numerical Algorithms in J*, The Education Vector, Volume 12, Number 2, 1995, pp. 26-30. - Hui, Roger,
*Cholesky Decomposition*, J Wiki Essay, 2005-10-14. - Hui, Roger,
*Triangular Matrix Inverse*, J Wiki Essay, 2005-10-27. - Hui, Roger,
*LU Decomposition*, J Wiki Essay, 2005-10-31. - Hui, Roger,
*QR Decomposition*, J Wiki Essay, 2005-10-30. - 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.

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

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.

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.