Cholesky decomposition: Difference between revisions
(add Ruby) |
(→{{header|Go}}: fixed. output was wrong. also added a few comments.) |
||
Line 576: | Line 576: | ||
"math" |
"math" |
||
) |
) |
||
// symmetric and lower use a packed representation that stores only |
|||
// the lower triangle. |
|||
type symmetric struct { |
|||
order int |
|||
ele []float64 |
|||
} |
|||
type lower struct { |
type lower struct { |
||
Line 582: | Line 590: | ||
} |
} |
||
// symmetric.print prints a square matrix from the packed representation, |
|||
// printing the upper triange as a transpose of the lower. |
|||
func (s *symmetric) print() { |
|||
const eleFmt = "%10.5f " |
|||
row, diag := 1, 0 |
|||
for i, e := range s.ele { |
|||
fmt.Printf(eleFmt, e) |
|||
if i == diag { |
|||
for j, col := diag+row, row; col < s.order; j += col { |
|||
fmt.Printf(eleFmt, s.ele[j]) |
|||
col++ |
|||
} |
|||
fmt.Println() |
|||
row++ |
|||
diag += row |
|||
} |
|||
} |
|||
} |
|||
// lower.print prints a square matrix from the packed representation, |
|||
// printing the upper triangle as all zeros. |
|||
func (l *lower) print() { |
func (l *lower) print() { |
||
const eleFmt = "%10.5f " |
const eleFmt = "%10.5f " |
||
Line 598: | Line 627: | ||
} |
} |
||
// choleskyLower returns the cholesky decomposition of a symmetric real |
|||
⚫ | |||
// matrix. The matrix must be positive definite but this is not checked. |
|||
⚫ | |||
l := &lower{a.order, make([]float64, len(a.ele))} |
l := &lower{a.order, make([]float64, len(a.ele))} |
||
row, col := 1, 1 |
row, col := 1, 1 |
||
Line 627: | Line 658: | ||
func main() { |
func main() { |
||
demo(& |
demo(&symmetric{3, []float64{ |
||
25, |
25, |
||
15, 18, |
15, 18, |
||
-5, 0, 11}}) |
-5, 0, 11}}) |
||
demo(& |
demo(&symmetric{4, []float64{ |
||
18, |
18, |
||
22, 70, |
22, 70, |
||
Line 638: | Line 669: | ||
} |
} |
||
func demo(a * |
func demo(a *symmetric) { |
||
fmt.Println("A:") |
fmt.Println("A:") |
||
a.print() |
a.print() |
||
Line 647: | Line 678: | ||
<pre> |
<pre> |
||
A: |
A: |
||
25.00000 |
25.00000 15.00000 -5.00000 |
||
15.00000 18.00000 0.00000 |
15.00000 18.00000 0.00000 |
||
-5.00000 0.00000 11.00000 |
-5.00000 0.00000 11.00000 |
||
Line 655: | Line 686: | ||
-1.00000 1.00000 3.00000 |
-1.00000 1.00000 3.00000 |
||
A: |
A: |
||
18.00000 |
18.00000 22.00000 54.00000 42.00000 |
||
22.00000 70.00000 |
22.00000 70.00000 86.00000 62.00000 |
||
54.00000 86.00000 174.00000 |
54.00000 86.00000 174.00000 134.00000 |
||
42.00000 62.00000 134.00000 106.00000 |
42.00000 62.00000 134.00000 106.00000 |
||
L: |
L: |
Revision as of 04:28, 21 December 2011
You are encouraged to solve this task according to the task description, using any language you may know.
Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:
is called the Cholesky factor of , and can be interpreted as a generalized square root of , as described in Cholesky decomposition.
In a 3x3 example, we have to solve the following system of equations:
We can see that for the diagonal elements () of there is a calculation pattern:
or in general:
For the elements below the diagonal (, where ) there is also a calculation pattern:
which can also be expressed in a general formula:
Task description
The task is to implement a routine which will return a lower Cholesky factor for every given symmetric, positive definite nxn matrix . 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
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
Note: This specimen retains the original C coding style. diff
<lang algol68>#!/usr/local/bin/a68g --script #
MODE FIELD=LONG REAL; PROC (FIELD)FIELD field sqrt = long sqrt; INT field prec = 5; FORMAT field fmt = $g(-(2+1+field prec),field prec)$;
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; FOR j FROM i+1 TO 2 UPB a DO l[i,j]:=0 # Not required if matrix is declared as triangular # OD OD; l
);
PROC print matrix v1 =(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
);
PROC print matrix =(MAT a)VOID:(
FORMAT vector fmt = $"("f(field fmt)n(2 UPB a-2 LWB a)(", " f(field fmt))")"$; FORMAT matrix fmt = $"("f(vector fmt)n( UPB a- LWB a)(","lxf(vector fmt))")"$; printf((matrix fmt, a))
);
main: (
MAT m1 = ((25, 15, -5), (15, 18, 0), (-5, 0, 11)); MAT c1 = cholesky(m1); print 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); print 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>
<lang lisp>;; 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>
<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)
- 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>
<lang lisp>;; case of matrix stored as a list of lists (inner lists are rows of matrix)
- as above, returns the Cholesky decomposition matrix of a square positive-definite, symmetric matrix
(defun cholesky (m)
(let ((l (list (list (sqrt (caar m))))) x (j 0) i) (dolist (cm (cdr m) (mapcar #'(lambda (x) (nconc x (make-list (- (length m) (length x)) :initial-element 0))) l)) (setq x (list (/ (car cm) (caar l))) i 0) (dolist (cl (cdr l)) (setf (cdr (last x)) (list (/ (- (elt cm (incf i)) (*v x cl)) (car (last cl)))))) (setf (cdr (last l)) (list (nconc x (list (sqrt (- (elt cm (incf j)) (*v x x))))))))))
- where *v is the scalar product defined as
(defun *v (v1 v2) (reduce #'+ (mapcar #'* v1 v2)))</lang>
<lang lisp>;; example 1 CL-USER> (setf a '((25 15 -5) (15 18 0) (-5 0 11))) ((25 15 -5) (15 18 0) (-5 0 11)) CL-USER> (cholesky a) ((5 0 0) (3 3 0) (-1 1 3)) CL-USER> (format t "~{~{~5d~}~%~}" (cholesky a))
5 0 0 3 3 0 -1 1 3
NIL</lang>
<lang lisp>;; example 2 CL-USER> (setf a '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))) ((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106)) CL-USER> (cholesky a) ((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)) CL-USER> (format t "~{~{~10,5f~}~%~}" (cholesky a))
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.89950 1.62455 1.84971 1.39262
NIL</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]
DWScript
<lang delphi>function Cholesky(a : array of Float) : array of Float; var
i, j, k, n : Integer; s : Float;
begin
n:=Round(Sqrt(a.Length)); Result:=new Float[n*n]; for i:=0 to n-1 do begin for j:=0 to i do begin s:=0 ; for k:=0 to j-1 do s+=Result[i*n+k] * Result[j*n+k]; if i=j then Result[i*n+j]:=Sqrt(a[i*n+i]-s) else Result[i*n+j]:=1/Result[j*n+j]*(a[i*n+j]-s); end; end;
end;
procedure ShowMatrix(a : array of Float); var
i, j, n : Integer;
begin
n:=Round(Sqrt(a.Length)); for i:=0 to n-1 do begin for j:=0 to n-1 do Print(Format('%2.5f ', [a[i*n+j]])); PrintLn(); end;
end;
var m1 := new Float[9]; m1 := [ 25.0, 15.0, -5.0,
15.0, 18.0, 0.0, -5.0, 0.0, 11.0 ];
var c1 := Cholesky(m1); ShowMatrix(c1);
PrintLn();
var m2 : array of Float := [ 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 ];
var c2 := Cholesky(m2); ShowMatrix(c2);</lang>
Fantom
<lang fantom>
- Cholesky decomposition
class Main {
// create an array of Floats, initialised to 0.0 Float[][] makeArray (Int i, Int j) { Float[][] result := [,] i.times { result.add ([,]) } i.times |Int x| { j.times { result[x].add(0f) } } return result }
// perform the Cholesky decomposition Float[][] cholesky (Float[][] array) { m := array.size Float[][] l := makeArray (m, m) m.times |Int i| { (i+1).times |Int k| { Float sum := (0..<k).toList.reduce (0f) |Float a, Int j -> Float| { a + l[i][j] * l[k][j] } if (i == k) l[i][k] = (array[i][i]-sum).sqrt else l[i][k] = (1.0f / l[k][k]) * (array[i][k] - sum) } } return l }
Void runTest (Float[][] array) { echo (array) echo (cholesky (array)) } Void main () { runTest ([[25f,15f,-5f],[15f,18f,0f],[-5f,0f,11f]]) runTest ([[18f,22f,54f,42f],[22f,70f,86f,62f],[54f,86f,174f,134f],[42f,62f,134f,106f]]) }
} </lang>
Output:
[[25.0, 15.0, -5.0], [15.0, 18.0, 0.0], [-5.0, 0.0, 11.0]] [[5.0, 0.0, 0.0], [3.0, 3.0, 0.0], [-1.0, 1.0, 3.0]] [[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]] [[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]]
Go
<lang go>package main
import (
"fmt" "math"
)
// symmetric and lower use a packed representation that stores only // the lower triangle.
type symmetric struct {
order int ele []float64
}
type lower struct {
order int ele []float64
}
// symmetric.print prints a square matrix from the packed representation, // printing the upper triange as a transpose of the lower. func (s *symmetric) print() {
const eleFmt = "%10.5f " row, diag := 1, 0 for i, e := range s.ele { fmt.Printf(eleFmt, e) if i == diag { for j, col := diag+row, row; col < s.order; j += col { fmt.Printf(eleFmt, s.ele[j]) col++ } fmt.Println() row++ diag += row } }
}
// lower.print prints a square matrix from the packed representation, // printing the upper triangle as all zeros. func (l *lower) print() {
const eleFmt = "%10.5f " row, diag := 1, 0 for i, e := range l.ele { fmt.Printf(eleFmt, e) if i == diag { for j := row; j < l.order; j++ { fmt.Printf(eleFmt, 0.) } fmt.Println() row++ diag += row } }
}
// choleskyLower returns the cholesky decomposition of a symmetric real // matrix. The matrix must be positive definite but this is not checked. func (a *symmetric) choleskyLower() *lower {
l := &lower{a.order, make([]float64, len(a.ele))} row, col := 1, 1 dr := 0 // index of diagonal element at end of row dc := 0 // index of diagonal element at top of column for i, e := range a.ele { if i < dr { d := (e - l.ele[i]) / l.ele[dc] l.ele[i] = d ci, cx := col, dc for j := i + 1; j <= dr; j++ { cx += ci ci++ l.ele[j] += d * l.ele[cx] } col++ dc += col } else { l.ele[i] = math.Sqrt(e - l.ele[i]) row++ dr += row col = 1 dc = 0 } } return l
}
func main() {
demo(&symmetric{3, []float64{ 25, 15, 18, -5, 0, 11}}) demo(&symmetric{4, []float64{ 18, 22, 70, 54, 86, 174, 42, 62, 134, 106}})
}
func demo(a *symmetric) {
fmt.Println("A:") a.print() fmt.Println("L:") a.choleskyLower().print()
}</lang> Output:
A: 25.00000 15.00000 -5.00000 15.00000 18.00000 0.00000 -5.00000 0.00000 11.00000 L: 5.00000 0.00000 0.00000 3.00000 3.00000 0.00000 -1.00000 1.00000 3.00000 A: 18.00000 22.00000 54.00000 42.00000 22.00000 70.00000 86.00000 62.00000 54.00000 86.00000 174.00000 134.00000 42.00000 62.00000 134.00000 106.00000 L: 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
Haskell
We use the Cholesky–Banachiewicz algorithm described in the Wikipedia article.
For more serious numerical analysis there is a Cholesky decomposition function in the hmatrix package.
The Cholesky module: <lang haskell>module Cholesky (Arr, cholesky) where
import Data.Array.IArray import Data.Array.MArray import Data.Array.Unboxed import Data.Array.ST
type Idx = (Int,Int) type Arr = UArray Idx Double
-- Return the (i,j) element of the lower triangular matrix. (We assume the -- lower array bound is (0,0).) get :: Arr -> Arr -> Idx -> Double get a l (i,j) | i == j = sqrt $ a!(j,j) - dot
| i > j = (a!(i,j) - dot) / l!(j,j) | otherwise = 0 where dot = sum [l!(i,k) * l!(j,k) | k <- [0..j-1]]
-- Return the lower triangular matrix of a Cholesky decomposition. We assume -- the input is a real, symmetric, positive-definite matrix, with lower array -- bounds of (0,0). cholesky :: Arr -> Arr cholesky a = let n = maxBnd a
in runSTUArray $ do l <- thaw a mapM_ (update a l) [(i,j) | i <- [0..n], j <- [0..n]] return l where maxBnd = fst . snd . bounds update a l i = unsafeFreeze l >>= \l' -> writeArray l i (get a l' i)</lang>
The main module: <lang haskell>import Data.Array.IArray import Cholesky
ex1, ex2 :: Arr ex1 = listArray ((0,0),(2,2)) [25, 15, -5,
15, 18, 0, -5, 0, 11]
ex2 = listArray ((0,0),(3,3)) [18, 22, 54, 42,
22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106]
main :: IO () main = do
print $ elems $ cholesky ex1 print $ elems $ cholesky ex2</lang>
The resulting matrices are printed as lists, as in the following 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.899494936611665,1.6245538642137891,1.849711005231382,1.3926212476455924]
Icon and Unicon
<lang Icon> procedure cholesky (array)
result := make_square_array (*array) every (i := 1 to *array) do { every (k := 1 to i) do { sum := 0 every (j := 1 to (k-1)) do { sum +:= result[i][j] * result[k][j] } if (i = k) then result[i][k] := sqrt(array[i][i] - sum) else result[i][k] := 1.0 / result[k][k] * (array[i][k] - sum) } } return result
end
procedure make_square_array (n)
result := [] every (1 to n) do push (result, list(n, 0)) return result
end
procedure print_array (array)
every (row := !array) do { every writes (!row || " ") write () }
end
procedure do_cholesky (array)
write ("Input:") print_array (array) result := cholesky (array) write ("Result:") print_array (result)
end
procedure main ()
do_cholesky ([[25,15,-5],[15,18,0],[-5,0,11]]) do_cholesky ([[18,22,54,42],[22,70,86,62],[54,86,174,134],[42,62,134,106]])
end </lang>
Output:
Input: 25 15 -5 15 18 0 -5 0 11 Result: 5.0 0 0 3.0 3.0 0 -1.0 1.0 3.0 Input: 18 22 54 42 22 70 86 62 54 86 174 134 42 62 134 106 Result: 4.242640687 0 0 0 5.185449729 6.565905201 0 0 12.72792206 3.046038495 1.649742248 0 9.899494937 1.624553864 1.849711005 1.392621248
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
<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]]
Mathematica
<lang Mathematica>CholeskyDecomposition[{{25, 15, -5}, {15, 18, 0}, {-5, 0, 11}}] </lang>
OCaml
<lang OCaml>let cholesky inp =
let n = Array.length inp in let res = Array.make_matrix n n 0.0 in let factor i k = let rec sum j = if j = k then 0.0 else res.(i).(j) *. res.(k).(j) +. sum (j+1) in inp.(i).(k) -. sum 0 in for col = 0 to n-1 do res.(col).(col) <- sqrt (factor col col); for row = col+1 to n-1 do res.(row).(col) <- (factor row col) /. res.(col).(col) done done; res
let pr_vec v = Array.iter (Printf.printf " %9.5f") v; print_newline() let show = Array.iter pr_vec let test a =
print_endline "\nin:"; show a; print_endline "out:"; show (cholesky a)
let _ =
test [| [|25.0; 15.0; -5.0|]; [|15.0; 18.0; 0.0|]; [|-5.0; 0.0; 11.0|] |]; test [| [|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:
in: 25.00000 15.00000 -5.00000 15.00000 18.00000 0.00000 -5.00000 0.00000 11.00000 out: 5.00000 0.00000 0.00000 3.00000 3.00000 0.00000 -1.00000 1.00000 3.00000 in: 18.00000 22.00000 54.00000 42.00000 22.00000 70.00000 86.00000 62.00000 54.00000 86.00000 174.00000 134.00000 42.00000 62.00000 134.00000 106.00000 out: 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
Pascal
<lang pascal>Program Cholesky;
type
D2Array = array of array of double;
function cholesky(const A: D2Array): D2Array;
var i, j, k: integer; s: double; begin setlength(cholesky, length(A), length(A)); for i := low(cholesky) to high(cholesky) do for j := 0 to i do begin
s := 0; for k := 0 to j - 1 do s := s + cholesky[i][k] * cholesky[j][k]; if i = j then cholesky[i][j] := sqrt(A[i][i] - s) else
cholesky[i][j] := (A[i][j] - s) / cholesky[j][j]; // save one multiplication compared to the original end; end;
procedure printM(const A: D2Array);
var i, j: integer; begin for i := low(A) to high(A) do begin for j := low(A) to high(A) do write(A[i,j]:8:5); writeln; end; end;
const
m1: array[0..2,0..2] of double = ((25, 15, -5), (15, 18, 0),
(-5, 0, 11));
m2: array[0..3,0..3] of double = ((18, 22, 54, 42), (22, 70, 86, 62),
(54, 86, 174, 134), (42, 62, 134, 106)); var
index: integer; cIn, cOut: D2Array;
begin
setlength(cIn, length(m1), length(m1)); for index := low(m1) to high(m1) do cIn[index] := m1[index]; cOut := cholesky(cIn); printM(cOut);
writeln; setlength(cIn, length(m2), length(m2)); for index := low(m2) to high(m2) do cIn[index] := m2[index]; cOut := cholesky(cIn); printM(cOut);
end.</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
Perl 6
<lang perl6> sub cholesky(@A) {
my @L = @A »*» 0; for ^@A -> $i {
for 0..$i -> $j { @L[$i][$j] = ($i == $j ?? &sqrt !! 1/@L[$j][$j] * * )( @A[$i][$j] - [+] (@L[$i] Z* @L[$j])[0..$j] ); }
} return @L;
} .say for cholesky [
[25], [15, 18], [-5, 0, 11],
];
.say for cholesky [
[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106],
]; </lang>
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]]
Ruby
<lang ruby>require 'matrix'
class Matrix
def symmetric? return false if not square? (0 ... row_size).each do |i| (0 .. i).each do |j| return false if self[i,j] != self[j,i] end end true end
def cholesky_factor raise ArgumentError, "must provide symmetric matrix" unless symmetric? l = Array.new(row_size) {Array.new(row_size, 0)} (0 ... row_size).each do |k| (0 ... row_size).each do |i| if i == k sum = (0 .. k-1).inject(0.0) {|sum, j| sum + l[k][j] ** 2} val = Math.sqrt(self[k,k] - sum) l[k][k] = val elsif i > k sum = (0 .. k-1).inject(0.0) {|sum, j| sum + l[i][j] * l[k][j]} val = (self[k,i] - sum) / l[k][k] l[i][k] = val end end end Matrix[*l] end
end
puts Matrix[[25,15,-5],[15,18,0],[-5,0,11]].cholesky_factor puts Matrix[[18, 22, 54, 42],
[22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]].cholesky_factor</lang>
outputs
Matrix[[5.0, 0, 0], [3.0, 3.0, 0], [-1.0, 1.0, 3.0]] Matrix[[4.242640687119285, 0, 0, 0], [5.185449728701349, 6.565905201197403, 0, 0], [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0], [9.899494936611665, 1.6245538642137891, 1.849711005231382, 1.3926212476455924]]
Tcl
<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}