Cholesky decomposition: Difference between revisions

Content added Content deleted

Inline

Revision as of 11:18, 19 May 2011

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

Ada

Works with: Ada 2005

decomposition.ads: <lang Ada>with Ada.Numerics.Generic_Real_Arrays; generic

  with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);

package Decomposition is

  -- decompose a square matrix A by A = L * Transpose (L)
  procedure Decompose (A : Matrix.Real_Matrix; L : out Matrix.Real_Matrix);

end Decomposition;</lang>

decomposition.adb: <lang Ada>with Ada.Numerics.Generic_Elementary_Functions;

package body Decomposition is

  package Math is new Ada.Numerics.Generic_Elementary_Functions
    (Matrix.Real);

  procedure Decompose (A : Matrix.Real_Matrix; L : out Matrix.Real_Matrix) is
     use type Matrix.Real_Matrix, Matrix.Real;
     Order : constant Positive := A'Length (1);
     S     : Matrix.Real;
  begin
     L := (others => (others => 0.0));
     for I in 0 .. Order - 1 loop
        for K in 0 .. I loop
           S := 0.0;
           for J in 0 .. K - 1 loop
              S := S +
                L (L'First (1) + I, L'First (2) + J) *
                L (L'First (1) + K, L'First (2) + J);
           end loop;
           -- diagonals
           if K = I then
              L (L'First (1) + K, L'First (2) + K) :=
                Math.Sqrt (A (A'First (1) + K, A'First (2) + K) - S);
           else
              L (L'First (1) + I, L'First (2) + K) :=
                1.0 / L (L'First (1) + K, L'First (2) + K) *
                (A (A'First (1) + I, A'First (2) + K) - S);
           end if;
        end loop;
     end loop;
  end Decompose;

end Decomposition;</lang>

Example usage: <lang Ada>with Ada.Numerics.Real_Arrays; with Ada.Text_IO; with Decomposition; procedure Decompose_Example is

  package Real_Decomposition is new Decomposition
    (Matrix => Ada.Numerics.Real_Arrays);

  package Real_IO is new Ada.Text_IO.Float_IO (Float);

  procedure Print (M : Ada.Numerics.Real_Arrays.Real_Matrix) is
  begin
     for Row in M'Range (1) loop
        for Col in M'Range (2) loop
           Real_IO.Put (M (Row, Col), 4, 3, 0);
        end loop;
        Ada.Text_IO.New_Line;
     end loop;
  end Print;

  Example_1 : constant Ada.Numerics.Real_Arrays.Real_Matrix :=
    ((25.0, 15.0, -5.0),
     (15.0, 18.0, 0.0),
     (-5.0, 0.0, 11.0));
  L_1 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_1'Range (1),
                                              Example_1'Range (2));
  Example_2 : constant Ada.Numerics.Real_Arrays.Real_Matrix :=
    ((18.0, 22.0, 54.0, 42.0),
     (22.0, 70.0, 86.0, 62.0),
     (54.0, 86.0, 174.0, 134.0),
     (42.0, 62.0, 134.0, 106.0));
  L_2 : Ada.Numerics.Real_Arrays.Real_Matrix (Example_2'Range (1),
                                              Example_2'Range (2));

begin

  Real_Decomposition.Decompose (A => Example_1,
                                L => L_1);
  Real_Decomposition.Decompose (A => Example_2,
                                L => L_2);
  Ada.Text_IO.Put_Line ("Example 1:");
  Ada.Text_IO.Put_Line ("A:"); Print (Example_1);
  Ada.Text_IO.Put_Line ("L:"); Print (L_1);
  Ada.Text_IO.New_Line;
  Ada.Text_IO.Put_Line ("Example 2:");
  Ada.Text_IO.Put_Line ("A:"); Print (Example_2);
  Ada.Text_IO.Put_Line ("L:"); Print (L_2);

end Decompose_Example;</lang>

Output:

Example 1:
A:
  25.000  15.000  -5.000
  15.000  18.000   0.000
  -5.000   0.000  11.000
L:
   5.000   0.000   0.000
   3.000   3.000   0.000
  -1.000   1.000   3.000

Example 2:
A:
  18.000  22.000  54.000  42.000
  22.000  70.000  86.000  62.000
  54.000  86.000 174.000 134.000
  42.000  62.000 134.000 106.000
L:
   4.243   0.000   0.000   0.000
   5.185   6.566   0.000   0.000
  12.728   3.046   1.650   0.000
   9.899   1.625   1.850   1.393

ALGOL 68

Translation of: C

Note: This specimen retains the original C coding style. diff

Works with: ALGOL 68 version Revision 1 - no extensions to language used.

Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny.

<lang algol68>MODE FIELD=LONG REAL; PROC (FIELD)FIELD field sqrt = long sqrt; INT field prec = 5;

MODE MAT = [0,0]FIELD;

PROC cholesky = (MAT a) MAT:(

   [UPB a, 2 UPB a]FIELD l;

   FOR i FROM LWB a TO UPB a DO
       FOR j FROM 2 LWB a TO i DO
           FIELD s := 0;
           FOR k FROM 2 LWB a TO j-1 DO
               s +:= l[i,k] * l[j,k]
           OD;
           l[i,j] := IF i = j 
                     THEN field sqrt(a[i,i] - s) 
                     ELSE 1.0 / l[j,j] * (a[i,j] - s) FI
       OD
   OD;
   l

);

PROC show matrix=(MAT a)VOID:(

   FOR i FROM LWB a TO UPB a DO
       FOR j FROM 2 LWB a TO 2 UPB a DO
           printf(($g(-(2+1+field prec),field prec)$, a[i,j]))
       OD;
       printf($l$)
   OD

);

main: (

   MAT m1 = ((25, 15, -5),
             (15, 18,  0),
             (-5,  0, 11));
   MAT c1 = cholesky(m1);
   show matrix(c1);
   printf($l$);

   MAT m2 = ((18, 22,  54,  42),
             (22, 70,  86,  62),
             (54, 86, 174, 134),
             (42, 62, 134, 106));
   MAT c2 = cholesky(m2);
   show matrix(c2)

)</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

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, std.numeric;

T[][] cholesky(T)(in T[][] A) {

   auto L = new T[][](A.length, A.length);
   foreach (r, row; L)
       row[r+1 .. $] = 0;
   foreach (i; 0 .. A.length)
       foreach (j; 0 .. i+1) {
           T t = dotProduct(L[i][0..j], L[j][0..j]);
           L[i][j] = (i == j) ? (A[i][i] - t) ^^ 0.5 :
                                (1.0 / L[j][j] * (A[i][j] - t));
       }
   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
  cholesky eg1
5 0 0
3 3 0

_1 1 3

  cholesky 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; double[][] l = new double[m][m]; //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]]

PicoLisp

<lang PicoLisp>(scl 9) (load "@lib/math.l")

(de cholesky (A)

  (let L (mapcar '(() (need (length A) 0)) A)
     (for (I . R) A
        (for J I
           (let S (get R J)
              (for K (inc J)
                 (dec 'S (*/ (get L I K) (get L J K) 1.0)) )
              (set (nth L I J)
                 (if (= I J)
                    (sqrt (* 1.0 S))
                    (*/ S 1.0 (get L J J)) ) ) ) ) )
     (for R L
        (for N R (prin (align 9 (round N 5))))
        (prinl) ) ) )</lang>

Test: <lang PicoLisp>(cholesky

  '((25.0 15.0 -5.0) (15.0 18.0 0) (-5.0 0 11.0)) )

(prinl)

(cholesky

  (quote
     (18.0  22.0   54.0   42.0)
     (22.0  70.0   86.0   62.0)
     (54.0  86.0  174.0  134.0)
     (42.0  62.0  134.0  106.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

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}