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_{21}+l_{32}l_{22}\\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

C

<lang c>#include <stdio.h>

include <stdlib.h>
include <math.h>

double *cholesky(double *A, int n) {

   double *L = (double*)calloc(n * n, sizeof(double));
   if (L == NULL)
       exit(EXIT_FAILURE);

   for (int i = 0; i < n; i++)
       for (int j = 0; j < (i+1); j++) {
           double s = 0;
           for (int k = 0; k < j; k++)
               s += L[i * n + k] * L[j * n + k];
           L[i * n + j] = (i == j) ?
                          sqrt(A[i * n + i] - s) :
                          (1.0 / L[j * n + j] * (A[i * n + j] - s));
       }

   return L;

}

void show_matrix(double *A, int n) {

   for (int i = 0; i < n; i++) {
       for (int j = 0; j < n; j++)
           printf("%2.5f ", A[i * n + j]);
       printf("\n");
   }

}

int main() {

   int n = 3;
   double m1[] = {25, 15, -5,
                  15, 18,  0,
                  -5,  0, 11};
   double *c1 = cholesky(m1, n);
   show_matrix(c1, n);
   printf("\n");
   free(c1);

   n = 4;
   double m2[] = {18, 22,  54,  42,
                  22, 70,  86,  62,
                  54, 86, 174, 134,
                  42, 62, 134, 106};
   double *c2 = cholesky(m2, n);
   show_matrix(c2, n);
   free(c2);

   return 0;

}</lang> Output:

5.00000 0.00000 0.00000
3.00000 3.00000 0.00000
-1.00000 1.00000 3.00000

4.24264 0.00000 0.00000 0.00000
5.18545 6.56591 0.00000 0.00000
12.72792 3.04604 1.64974 0.00000
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>

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]

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>

Java

Works with: Java version 1.5+

<lang java5>import java.util.Arrays;

public class Cholesky { public static double[][] chol(double[][] a){ int m = a.length; int n = a[0].length; double[][] l = new double[m][n]; //automatically initialzed to 0's for(int i = 0; i< m;i++){ for(int k = 0; k < (i+1); k++){ double sum = 0; for(int j = 0; j < k; j++){ sum += l[i][j] * l[k][j]; } l[i][k] = (i == k) ? Math.sqrt(a[i][i] - sum) : (1.0 / l[k][k] * (a[i][k] - sum)); } } return l; }

public static void main(String[] args){ double[][] test1 = {{25, 15, -5}, {15, 18, 0}, {-5, 0, 11}}; System.out.println(Arrays.deepToString(chol(test1))); double[][] test2 = {{18, 22, 54, 42}, {22, 70, 86, 62}, {54, 86, 174, 134}, {42, 62, 134, 106}}; System.out.println(Arrays.deepToString(chol(test2))); } }</lang> Output:

[[5.0, 0.0, 0.0], [3.0, 3.0, 0.0], [-1.0, 1.0, 3.0]]
[[4.242640687119285, 0.0, 0.0, 0.0], [5.185449728701349, 6.565905201197403, 0.0, 0.0], [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0.0], [9.899494936611667, 1.624553864213788, 1.8497110052313648, 1.3926212476456026]]

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

Tcl

Translation of: Java

<lang tcl>proc cholesky a {

   set m [llength $a]
   set n [llength [lindex $a 0]]
   set l [lrepeat $m [lrepeat $n 0.0]]
   for {set i 0} {$i < $m} {incr i} {

for {set k 0} {$k < $i+1} {incr k} { set sum 0.0 for {set j 0} {$j < $k} {incr j} { set sum [expr {$sum + [lindex $l $i $j] * [lindex $l $k $j]}] } lset l $i $k [expr { $i == $k ? sqrt([lindex $a $i $i] - $sum) : (1.0 / [lindex $l $k $k] * ([lindex $a $i $k] - $sum)) }] }

   }
   return $l

}</lang> Demonstration code: <lang tcl>set test1 {

   {25 15 -5}
   {15 18  0}
   {-5  0 11}

} puts [cholesky $test1] set test2 {

   {18 22  54  42}
   {22 70  86  62}
   {54 86 174 134}
   {42 62 134 106}

} puts [cholesky $test2]</lang> Output:

{5.0 0.0 0.0} {3.0 3.0 0.0} {-1.0 1.0 3.0}
{4.242640687119285 0.0 0.0 0.0} {5.185449728701349 6.565905201197403 0.0 0.0} {12.727922061357857 3.0460384954008553 1.6497422479090704 0.0} {9.899494936611667 1.624553864213788 1.8497110052313648 1.3926212476456026}