Cholesky decomposition

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>

Example 1

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

2A((5.0 0 0)

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

</lang>

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)

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>

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]]