Cholesky decomposition
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]]
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]]
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 a = Array.iter pr_vec a 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
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
<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}