Cholesky decomposition

Cholesky decomposition
You are encouraged to solve this task according to the task description, using any language you may know.

Method description

Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:

$A=LL^{T}$

$L$ is called the Cholesky factor of $A$, and can be interpreted as a generalized square root of $A$, as described in Cholesky decomposition.

In a 3x3 example, we have to solve the following system of equations:

$A={\begin{pmatrix}a_{11}&a_{21}&a_{31}\\a_{21}&a_{22}&a_{32}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\\\end{pmatrix}}{\begin{pmatrix}l_{11}&l_{21}&l_{31}\\0&l_{22}&l_{32}\\0&0&l_{33}\end{pmatrix}}={\begin{pmatrix}l_{11}^{2}&l_{21}l_{11}&l_{31}l_{11}\\l_{21}l_{11}&l_{21}^{2}+l_{22}^{2}&\\l_{31}l_{11}&l_{31}l_{21}+l_{32}l_{22}&l_{31}^{2}+l_{32}^{2}+l_{33}^{2}\end{pmatrix}}=LL^{T}$

We can see that for the diagonal elements ($l_{kk}$) of $L$ there is a calculation pattern:

$l_{11}={\sqrt {a_{11}}}$

$l_{22}={\sqrt {a_{22}-l_{21}^{2}}}$

$l_{33}={\sqrt {a_{33}-(l_{31}^{2}+l_{32}^{2})}}$

or in general:

$l_{kk}={\sqrt {a_{kk}-\sum _{j=1}^{k-1}l_{kj}^{2}}}$

For the elements below the diagonal ($l_{ik}$, where $i>k$) there is also a calculation pattern:

$l_{21}={\frac {1}{l_{11}}}a_{21}$

$l_{31}={\frac {1}{l_{11}}}a_{31}$

$l_{32}={\frac {1}{l_{22}}}(a_{32}-l_{31}l_{21})$

which can also be expressed in a general formula:

$l_{ik}={\frac {1}{l_{kk}}}\left(a_{ik}-\sum _{j=1}^{k-1}l_{ij}l_{kj}\right)$

Task description

The task is to implement a routine which will return a lower Cholesky factor $L$ for every given symmetric, positive definite nxn matrix $A$. You should then test it on the following two examples and include your output.

Example 1:

25  15  -5                 5   0   0
15  18   0         -->     3   3   0
-5   0  11                -1   1   3

Example 2:

18  22   54   42           4.24264    0.00000    0.00000    0.00000
22  70   86   62   -->     5.18545    6.56591    0.00000    0.00000
54  86  174  134          12.72792    3.04604    1.64974    0.00000
42  62  134  106           9.89949    1.62455    1.84971    1.39262

Common Lisp

<lang lisp>;; Calculates the Cholesky decomposition matrix L

for a positive-definite, symmetric nxn matrix A.

(defun chol (A)

 (let* ((n (car (array-dimensions A)))
        (L (make-array `(,n ,n) :initial-element 0)))

   (do ((k 0 (incf k))) ((> k (- n 1)) nil)
       ;; First, calculate diagonal elements L_kk.
       (setf (aref L k k)
             (sqrt (- (aref A k k)
                      (do* ((j 0 (incf j))
                            (sum (expt (aref L k j) 2) 
                                 (incf sum (expt (aref L k j) 2))))
                           ((> j (- k 1)) sum)))))

       ;; Then, all elements below a diagonal element, L_ik, i=k+1..n.
       (do ((i (+ k 1) (incf i)))
           ((> i (- n 1)) nil)

           (setf (aref L i k)
                 (/ (- (aref A i k)
                       (do* ((j 0 (incf j))
                             (sum (* (aref L i j) (aref L k j))
                                  (incf sum (* (aref L i j) (aref L k j)))))
                            ((> j (- k 1)) sum)))
                    (aref L k k)))))

   ;; Return the calculated matrix L.
   L))</lang>

<lang lisp>

Example 1

(setf A (make-array '(3 3) :initial-contents '((25 15 -5) (15 18 0) (-5 0 11)))) (chol A)

1. 2A((5.0 0 0)

   (3.0 3.0 0)
   (-1.0 1.0 3.0))

</lang>

<lang lisp>

Example 2

(setf B (make-array '(4 4) :initial-contents '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106)))) (chol B)

1. 2A((4.2426405 0 0 0)

   (5.18545 6.565905 0 0)
   (12.727922 3.0460374 1.6497375 0)
   (9.899495 1.6245536 1.849715 1.3926151))

</lang>

J

Solution: <lang j>mp=: +/ . * NB. matrix product h =: +@|: NB. conjugate transpose

cholesky=: 3 : 0

n=. #A=. y
if. 1>:n do.
 assert. (A=|A)>0=A  NB. check for positive definite
 %:A
else.
 p=. >.n%2 [ q=. <.n%2
 X=. (p,p) {.A [ Y=. (p,-q){.A [ Z=. (-q,q){.A
 L0=. cholesky X
 L1=. cholesky Z-(T=.(h Y) mp %.X) mp Y
 L0,(T mp L0),.L1
end.

)</lang> See Choleski Decomposition essay on the J Wiki.

Examples: <lang j> eg1=: 25 15 _5 , 15 18 0 ,: _5 0 11

  eg2=: 18 22 54 42 , 22 70 86 62 , 54 86 174 134 ,: 42 62 134 106
  Choleski eg1
5 0 0
3 3 0

_1 1 3

  Choleski eg2

4.24264 0 0 0 5.18545 6.56591 0 0 12.7279 3.04604 1.64974 0 9.89949 1.62455 1.84971 1.39262</lang>

D

<lang d>import std.stdio, std.math;

/// Calculates the Cholesky decomposition matrix L /// for a positive-definite, symmetric nxn matrix A. pure T[][] cholesky(T)(const T[][] A) {

   int n = A.length;
   auto L = new T[][](n, n);
   foreach (r, row; L)
       row[r+1 .. $] = 0;

   foreach (i; 0 .. n)
       foreach (j; 0 .. i+1) {
           T sum = 0;
           foreach (k; 0 .. j)
               sum += L[i][k] * L[j][k];
           L[i][j] = (i == j) ? sqrt(A[i][i] - sum) :
                                (1.0 / L[j][j] * (A[i][j] - sum));
       }

   return L;

}

void main() {

   double[][] m1 = [[25, 15, -5],
                    [15, 18,  0],
                    [-5,  0, 11]];
   foreach (row; cholesky(m1))
       writeln(row);
   writeln();

   double[][] m2 = [[18, 22,  54,  42],
                    [22, 70,  86,  62],
                    [54, 86, 174, 134],
                    [42, 62, 134, 106]];
   foreach (row; cholesky(m2))
       writeln(row);

}</lang> Output:

[5, 0, 0]
[3, 3, 0]
[-1, 1, 3]

[4.24264, 0, 0, 0]
[5.18545, 6.56591, 0, 0]
[12.7279, 3.04604, 1.64974, 0]
[9.89949, 1.62455, 1.84971, 1.39262]

Python

<lang python>import math, pprint

def cholesky(A):

   L = [[0.0] * len(A) for _ in xrange(len(A))]
   for i in xrange(len(A)):
       for j in xrange(i+1):
           s = sum(L[i][k] * L[j][k] for k in xrange(j))
           L[i][j] = math.sqrt(A[i][i] - s) if (i == j) else \
                     (1.0 / L[j][j] * (A[i][j] - s))
   return L

m1 = [[25, 15, -5],

     [15, 18,  0],
     [-5,  0, 11]]

pprint.pprint(cholesky(m1)) print

m2 = [[18, 22, 54, 42],

     [22, 70,  86,  62],
     [54, 86, 174, 134],
     [42, 62, 134, 106]]

pprint.pprint(cholesky(m2))</lang> Output:

[[5.0, 0.0, 0.0], [3.0, 3.0, 0.0], [-1.0, 1.0, 3.0]]

[[4.2426406871192848, 0.0, 0.0, 0.0],
 [5.1854497287013492, 6.5659052011974026, 0.0, 0.0],
 [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0.0],
 [9.8994949366116671,
  1.624553864213788,
  1.8497110052313648,
  1.3926212476456026]]