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.

Follow

Quicksort is a classic sorting algorithm invented by C.A.R. Hoare in 1961 [0, 1]. It has been known for some time that quicksort has a terse rendition in APL [2]. To get right to it, here is the code:

Q←{1≥≢⍵:⍵ ⋄ S←{⍺⌿⍨⍺ ⍺⍺ ⍵} ⋄ ⍵((∇<S)⍪=S⍪(∇>S))⍵⌷⍨?≢⍵}

The “pivot” ⍵⌷⍨?≢⍵ is randomly chosen. ((∇<S)⍪=S⍪(∇>S)) is a fork, selecting the parts of ⍵ which are less than the pivot, equal to the pivot, and greater than the pivot. The function is recursively applied to the first and the last of these three parts.

      ⎕io←0 ⋄ ⎕rl←7*5

      ⎕←x←?13⍴20
3 2 19 16 11 4 18 17 9 17 7 3 1
      Q x
1 2 3 3 4 7 9 11 16 17 17 18 19

The variant Q1 obtains by enclosing each of the three parts. Its result exhibits an interesting structure. The middle item of each triplet is the value of the pivot at each recursion. Since the pivot is randomly chosen, the result of Q1 can be different on the same argument, as illustrated below:

      Q1←{1≥≢⍵:⍵ ⋄ S←{⍺⌿⍨⍺ ⍺⍺ ⍵} ⋄ ⍵((⊂∘∇<S)⍪(⊂=S)⍪(⊂∘∇>S))⍵⌷⍨?≢⍵}

      Q1 x
┌──────┬───┬────────────────────────────────────┐
│┌─┬─┬┐│3 3│┌─┬─┬──────────────────────────────┐│
││1│2│││   ││4│7│┌┬─┬─────────────────────────┐││
│└─┴─┴┘│   ││ │ │││9│┌──┬──┬─────────────────┐│││
│      │   ││ │ │││ ││11│16│┌─────────┬──┬──┐││││
│      │   ││ │ │││ ││  │  ││┌┬─────┬┐│18│19│││││
│      │   ││ │ │││ ││  │  ││││17 17│││  │  │││││
│      │   ││ │ │││ ││  │  ││└┴─────┴┘│  │  │││││
│      │   ││ │ │││ ││  │  │└─────────┴──┴──┘││││
│      │   ││ │ │││ │└──┴──┴─────────────────┘│││
│      │   ││ │ │└┴─┴─────────────────────────┘││
│      │   │└─┴─┴──────────────────────────────┘│
└──────┴───┴────────────────────────────────────┘

      Q1 x
┌───────────────────────┬─┬─────────────────────────┐
│┌──────────────────┬─┬┐│9│┌──┬──┬─────────────────┐│
││┌┬─┬─────────────┐│7│││ ││11│16│┌───────────┬──┬┐││
││││1│┌┬─┬────────┐││ │││ ││  │  ││┌┬─────┬──┐│19││││
││││ │││2│┌┬───┬─┐│││ │││ ││  │  ││││17 17│18││  ││││
││││ │││ │││3 3│4││││ │││ ││  │  ││└┴─────┴──┘│  ││││
││││ │││ │└┴───┴─┘│││ │││ ││  │  │└───────────┴──┴┘││
││││ │└┴─┴────────┘││ │││ │└──┴──┴─────────────────┘│
││└┴─┴─────────────┘│ │││ │                         │
│└──────────────────┴─┴┘│ │                         │
└───────────────────────┴─┴─────────────────────────┘

The enlist of the result of Q1 x is the same as Q x, the sort of x:

      ∊ Q1 x
1 2 3 3 4 7 9 11 16 17 17 18 19
      Q x
1 2 3 3 4 7 9 11 16 17 17 18 19

This note is meant to explore the workings of a classical algorithm. To actually sort data in Dyalog, it is more convenient and more efficient to use {⍵[⍋⍵]}. Earlier versions of this text appeared in [3, 4].

References

1. Hoare, C.A.R, Algorithm 63: Partition, Communications of the ACM, Volume 4, Number 7, 1961-07.
2. Hoare, C.A.R, Algorithm 64: Quicksort, Communications of the ACM, Volume 4, Number 7, 1961-07.
3. Hui, Roger K.W., and Kenneth E. Iverson, J Introduction and Dictionary, 1991-2014; if. entry.
4. Hui, Roger K.W., Quicksort, J Wiki Essay, 2005-09-28.
5. Hui, Roger K.W. Sixteen APL Amuse-Bouches, 2014-11-02.

Follow