I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Cholesky decomposition

(Redirected from Cholesky Decomposition)
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:

${\displaystyle A=LL^{T}}$

${\displaystyle L}$ is called the Cholesky factor of ${\displaystyle A}$, and can be interpreted as a generalized square root of ${\displaystyle A}$, as described in Cholesky decomposition.

In a 3x3 example, we have to solve the following system of equations:

{\displaystyle {\begin{aligned}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}}\equiv LL^{T}\\&={\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}}\end{aligned}}}

We can see that for the diagonal elements (${\displaystyle l_{kk}}$) of ${\displaystyle L}$ there is a calculation pattern:

${\displaystyle l_{11}={\sqrt {a_{11}}}}$
${\displaystyle l_{22}={\sqrt {a_{22}-l_{21}^{2}}}}$
${\displaystyle l_{33}={\sqrt {a_{33}-(l_{31}^{2}+l_{32}^{2})}}}$

or in general:

${\displaystyle l_{kk}={\sqrt {a_{kk}-\sum _{j=1}^{k-1}l_{kj}^{2}}}}$

For the elements below the diagonal (${\displaystyle l_{ik}}$, where ${\displaystyle i>k}$) there is also a calculation pattern:

${\displaystyle l_{21}={\frac {1}{l_{11}}}a_{21}}$
${\displaystyle l_{31}={\frac {1}{l_{11}}}a_{31}}$
${\displaystyle l_{32}={\frac {1}{l_{22}}}(a_{32}-l_{31}l_{21})}$

which can also be expressed in a general formula:

${\displaystyle l_{ik}={\frac {1}{l_{kk}}}\left(a_{ik}-\sum _{j=1}^{k-1}l_{ij}l_{kj}\right)}$

The task is to implement a routine which will return a lower Cholesky factor ${\displaystyle L}$ for every given symmetric, positive definite nxn matrix ${\displaystyle 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


Note
1. The Cholesky decomposition of a Pascal upper-triangle matrix is the Identity matrix of the same size.
2. The Cholesky decomposition of a Pascal symmetric matrix is the Pascal lower-triangle matrix of the same size.

## 11l

Translation of: Python
F cholesky(A)   V l = [[0.0] * A.len] * A.len   L(i) 0 .< A.len      L(j) 0 .. i         V s = sum((0 .< j).map(k -> @l[@i][k] * @l[@j][k]))         l[i][j] = I (i == j) {sqrt(A[i][i] - s)} E (1.0 / l[j][j] * (A[i][j] - s))   R l F pprint(m)   print(‘[’)   L(row) m      print(row)   print(‘]’) V m1 = [[25, 15, -5],        [15, 18,  0],        [-5,  0, 11]]print(cholesky(m1))print() V m2 = [[18, 22,  54,  42],        [22, 70,  86,  62],        [54, 86, 174, 134],        [42, 62, 134, 106]]pprint(cholesky(m2))
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]
]


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;

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;

Example usage:

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;
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.
#!/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))
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))


## AutoHotkey

Cholesky_Decomposition(A){    L := [], n := A.Count()    L[1,1] := Sqrt(A[1,1])    loop % n {        k := A_Index        loop % n-1 {            i := A_Index+1             Sigma := 0, j := 0            while (++j <= k-1)                Sigma += L[i, j] * L[k, j]            L[i, k] := (A[i, k] - Sigma) / L[k, k]             Sigma := 0, j := 0            while (++j <= k-1)                Sigma += (L[k, j])**2                L[k, k] := Sqrt(A[k, k] - Sigma)        }    }    loop % n{        k := A_Index        loop % n            L[k, A_Index] := L[k, A_Index] ? L[k, A_Index] : 0    }    return L}ShowMatrix(L){    for r, obj in L{        row := ""        for c, v in obj            row .= Format("{:.3f}", v) ", "        output .= "[" trim(row, ", ") "]n,"    }    return "[" Trim(output, "n,") "]"}
Examples:
A := [[25, 15, -5]    , [15, 18,  0]    , [-5, 0 , 11]]L1 := Cholesky_Decomposition(A) A := [[18, 22,  54,  42]    , [22, 70,  86,  62]    , [54, 86, 174, 134]    , [42, 62, 134, 106]]L2 := Cholesky_Decomposition(A) MsgBox % Result := ShowMatrix(L1) "n----n" ShowMatrix(L2) "n----"return
Output:
[[5.000, 0.000, 0.000]
,[3.000, 3.000, 0.000]
,[-1.000, 1.000, 3.000]]
----
[[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]]
----

## BBC BASIC

      DIM m1(2,2)      m1() = 25, 15, -5, \      \      15, 18,  0, \      \      -5,  0, 11      PROCcholesky(m1())      PROCprint(m1())      PRINT       @% = &2050A      DIM m2(3,3)      m2() = 18, 22,  54,  42, \      \      22, 70,  86,  62, \      \      54, 86, 174, 134, \      \      42, 62, 134, 106      PROCcholesky(m2())      PROCprint(m2())      END       DEF PROCcholesky(a())      LOCAL i%, j%, k%, l(), s      DIM l(DIM(a(),1),DIM(a(),2))      FOR i% = 0 TO DIM(a(),1)        FOR j% = 0 TO i%          s = 0          FOR k% = 0 TO j%-1            s += l(i%,k%) * l(j%,k%)          NEXT          IF i% = j% THEN            l(i%,j%) = SQR(a(i%,i%) - s)          ELSE            l(i%,j%) = (a(i%,j%) - s) / l(j%,j%)          ENDIF        NEXT j%      NEXT i%      a() = l()      ENDPROC       DEF PROCprint(a())      LOCAL row%, col%      FOR row% = 0 TO DIM(a(),1)        FOR col% = 0 TO DIM(a(),2)          PRINT a(row%,col%);        NEXT        PRINT      NEXT row%      ENDPROC

Output:

         5         0         0
3         3         0
-1         1         3

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

#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;}
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#

using System;using System.Collections.Generic;using System.Linq;using System.Text; namespace Cholesky{    class Program    {        /// <summary>        /// This is example is written in C#, and compiles with .NET Framework 4.0        /// </summary>        /// <param name="args"></param>        static void Main(string[] args)        {            double[,] test1 = new double[,]            {                {25, 15, -5},                {15, 18, 0},                {-5, 0, 11},            };             double[,] test2 = new double[,]            {                {18, 22, 54, 42},                {22, 70, 86, 62},                {54, 86, 174, 134},                {42, 62, 134, 106},            };             double[,] chol1 = Cholesky(test1);            double[,] chol2 = Cholesky(test2);             Console.WriteLine("Test 1: ");            Print(test1);            Console.WriteLine("");            Console.WriteLine("Lower Cholesky 1: ");            Print(chol1);            Console.WriteLine("");            Console.WriteLine("Test 2: ");            Print(test2);            Console.WriteLine("");            Console.WriteLine("Lower Cholesky 2: ");            Print(chol2);         }         public static void Print(double[,] a)        {            int n = (int)Math.Sqrt(a.Length);             StringBuilder sb = new StringBuilder();            for (int r = 0; r < n; r++)            {                string s = "";                for (int c = 0; c < n; c++)                {                    s += a[r, c].ToString("f5").PadLeft(9) + ",";                }                sb.AppendLine(s);            }             Console.WriteLine(sb.ToString());        }         /// <summary>        /// Returns the lower Cholesky Factor, L, of input matrix A.         /// Satisfies the equation: L*L^T = A.        /// </summary>        /// <param name="a">Input matrix must be square, symmetric,         /// and positive definite. This method does not check for these properties,        /// and may produce unexpected results of those properties are not met.</param>        /// <returns></returns>        public static double[,] Cholesky(double[,] a)        {            int n = (int)Math.Sqrt(a.Length);             double[,] ret = new double[n, n];            for (int r = 0; r < n; r++)                for (int c = 0; c <= r; c++)                {                    if (c == r)                    {                        double sum = 0;                        for (int j = 0; j < c; j++)                        {                            sum += ret[c, j] * ret[c, j];                        }                        ret[c, c] = Math.Sqrt(a[c, c] - sum);                    }                    else                    {                        double sum = 0;                        for (int j = 0; j < c; j++)                            sum += ret[r, j] * ret[c, j];                        ret[r, c] = 1.0 / ret[c, c] * (a[r, c] - sum);                    }                }             return ret;        }    }}
Output:

Test 1:

25.00000, 15.00000, -5.00000,
15.00000, 18.00000,  0.00000,
-5.00000,  0.00000, 11.00000,


Lower Cholesky 1:

 5.00000,  0.00000,  0.00000,
3.00000,  3.00000,  0.00000,
-1.00000,  1.00000,  3.00000,


Test 2:

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,


Lower Cholesky 2:

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

#include <cassert>#include <cmath>#include <iomanip>#include <iostream>#include <vector> template <typename scalar_type> class matrix {public:    matrix(size_t rows, size_t columns)        : rows_(rows), columns_(columns), elements_(rows * columns) {}     matrix(size_t rows, size_t columns, scalar_type value)        : rows_(rows), columns_(columns), elements_(rows * columns, value) {}     matrix(size_t rows, size_t columns,        const std::initializer_list<std::initializer_list<scalar_type>>& values)        : rows_(rows), columns_(columns), elements_(rows * columns) {        assert(values.size() <= rows_);        size_t i = 0;        for (const auto& row : values) {            assert(row.size() <= columns_);            std::copy(begin(row), end(row), &elements_[i]);            i += columns_;        }    }     size_t rows() const { return rows_; }    size_t columns() const { return columns_; }     const scalar_type& operator()(size_t row, size_t column) const {        assert(row < rows_);        assert(column < columns_);        return elements_[row * columns_ + column];    }    scalar_type& operator()(size_t row, size_t column) {        assert(row < rows_);        assert(column < columns_);        return elements_[row * columns_ + column];    }private:    size_t rows_;    size_t columns_;    std::vector<scalar_type> elements_;}; template <typename scalar_type>void print(std::ostream& out, const matrix<scalar_type>& a) {    size_t rows = a.rows(), columns = a.columns();    out << std::fixed << std::setprecision(5);    for (size_t row = 0; row < rows; ++row) {        for (size_t column = 0; column < columns; ++column) {            if (column > 0)                out << ' ';            out << std::setw(9) << a(row, column);        }        out << '\n';    }} template <typename scalar_type>matrix<scalar_type> cholesky_factor(const matrix<scalar_type>& input) {    assert(input.rows() == input.columns());    size_t n = input.rows();    matrix<scalar_type> result(n, n);    for (size_t i = 0; i < n; ++i) {        for (size_t k = 0; k < i; ++k) {            scalar_type value = input(i, k);            for (size_t j = 0; j < k; ++j)                value -= result(i, j) * result(k, j);            result(i, k) = value/result(k, k);        }        scalar_type value = input(i, i);        for (size_t j = 0; j < i; ++j)            value -= result(i, j) * result(i, j);        result(i, i) = std::sqrt(value);    }    return result;} void print_cholesky_factor(const matrix<double>& matrix) {    std::cout << "Matrix:\n";    print(std::cout, matrix);    std::cout << "Cholesky factor:\n";    print(std::cout, cholesky_factor(matrix));} int main() {    matrix<double> matrix1(3, 3,       {{25, 15, -5},        {15, 18, 0},        {-5, 0, 11}});    print_cholesky_factor(matrix1);     matrix<double> matrix2(4, 4,       {{18, 22, 54, 42},        {22, 70, 86, 62},        {54, 86, 174, 134},        {42, 62, 134, 106}});    print_cholesky_factor(matrix2);     return 0;}
Output:
Matrix:
25.00000  15.00000  -5.00000
15.00000  18.00000   0.00000
-5.00000   0.00000  11.00000
Cholesky factor:
5.00000   0.00000   0.00000
3.00000   3.00000   0.00000
-1.00000   1.00000   3.00000
Matrix:
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
Cholesky factor:
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


## Clojure

Translation of: Python
(defn cholesky  [matrix]  (let [n (count matrix)        A (to-array-2d matrix)        L (make-array Double/TYPE n n)]    (doseq [i (range n) j (range (inc i))]      (let [s (reduce + (for [k (range j)] (* (aget L i k) (aget L j k))))]        (aset L i j (if (= i j)                      (Math/sqrt (- (aget A i i) s))                      (* (/ 1.0 (aget L j j)) (- (aget A i j) s))))))    (vec (map vec L))))

Example:

(cholesky [[25 15 -5] [15 18 0] [-5 0 11]]);=> [[ 5.0 0.0 0.0];    [ 3.0 3.0 0.0];    [-1.0 1.0 3.0]] (cholesky [[18 22 54 42] [22 70 86 62] [54 86 174 134] [42 62 134 106]]);=> [[ 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]]

## Common 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))
;; 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))
;; 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))
;; 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)))
;; example 1CL-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    3NIL
;; example 2CL-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.39262NIL

## D

import std.stdio, std.math, std.numeric; T[][] cholesky(T)(in T[][] A) pure nothrow /*@safe*/ {    auto L = new T[][](A.length, A.length);    foreach (immutable r, row; L)        row[r + 1 .. $] = 0; foreach (immutable i; 0 .. A.length) foreach (immutable j; 0 .. i + 1) { auto 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() { immutable double[][] m1 = [[25, 15, -5], [15, 18, 0], [-5, 0, 11]]; writefln("%(%(%2.0f %)\n%)\n", m1.cholesky); immutable double[][] m2 = [[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]]; writefln("%(%(%2.3f %)\n%)", m2.cholesky);} Output:  5 0 0 3 3 0 -1 1 3 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 ## Delphi See Pascal. ## DWScript Translation of: C 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); ## F# open Microsoft.FSharp.Collections let cholesky a = let calc (a: float[,]) (l: float[,]) i j = let c1 j = let sum = List.sumBy (fun k -> l.[j, k] ** 2.0) [0..j - 1] sqrt (a.[j, j] - sum) let c2 i j = let sum = List.sumBy (fun k -> l.[i, k] * l.[j, k]) [0..j - 1] (1.0 / l.[j, j]) * (a.[i, j] - sum) if j > i then 0.0 else if i = j then c1 j else c2 i j let l = Array2D.zeroCreate (Array2D.length1 a) (Array2D.length2 a) Array2D.iteri (fun i j _ -> l.[i, j] <- calc a l i j) l l let printMat a = let arrow = (Array2D.length2 a |> float) / 2.0 |> int let c = cholesky a for row in 0..(Array2D.length1 a) - 1 do for col in 0..(Array2D.length2 a) - 1 do printf "%.5f,\t" a.[row, col] printf (if arrow = row then "--> \t" else "\t\t") for col in 0..(Array2D.length2 c) - 1 do printf "%.5f,\t" c.[row, col] printfn "" let ex1 = array2D [ [25.0; 15.0; -5.0]; [15.0; 18.0; 0.0]; [-5.0; 0.0; 11.0]] let ex2 = array2D [ [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]] printfn "ex1:"printMat ex1 printfn "ex2:"printMat ex2  Output: ex1: 25.00000, 15.00000, -5.00000, 5.00000, 0.00000, 0.00000, 15.00000, 18.00000, 0.00000, --> 3.00000, 3.00000, 0.00000, -5.00000, 0.00000, 11.00000, -1.00000, 1.00000, 3.00000, ex2: 18.00000, 22.00000, 54.00000, 42.00000, 4.24264, 0.00000, 0.00000, 0.00000, 22.00000, 70.00000, 86.00000, 62.00000, 5.18545, 6.56591, 0.00000, 0.00000, 54.00000, 86.00000, 174.00000, 134.00000, --> 12.72792, 3.04604, 1.64974, 0.00000, 42.00000, 62.00000, 134.00000, 106.00000, 9.89949, 1.62455, 1.84971, 1.39262,  ## 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]]) }} 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]]  ## Fortran Program Cholesky_decomp! *************************************************!! LBH @ ULPGC 06/03/2014! Compute the Cholesky decomposition for a matrix A! after the attached ! http://rosettacode.org/wiki/Cholesky_decomposition! note that the matrix A is complex since there might! be values, where the sqrt has complex solutions.! Here, only the real values are taken into account!*************************************************! implicit none INTEGER, PARAMETER :: m=3 !rowsINTEGER, PARAMETER :: n=3 !colsCOMPLEX, DIMENSION(m,n) :: A REAL, DIMENSION(m,n) :: LREAL :: sum1, sum2INTEGER i,j,k ! Assign values to the matrixA(1,:)=(/ 25, 15, -5 /) A(2,:)=(/ 15, 18, 0 /) A(3,:)=(/ -5, 0, 11 /)! !!!!!!!!!!!another example!!!!!!!! A(1,:) = (/ 18, 22, 54, 42 /) ! A(2,:) = (/ 22, 70, 86, 62 /) ! A(3,:) = (/ 54, 86, 174, 134 /) ! A(4,:) = (/ 42, 62, 134, 106 /) ! Initialize valuesL(1,1)=real(sqrt(A(1,1)))L(2,1)=A(2,1)/L(1,1)L(2,2)=real(sqrt(A(2,2)-L(2,1)*L(2,1)))L(3,1)=A(3,1)/L(1,1)! for greater order than m,n=3 add initial row value! for instance if m,n=4 then add the following line! L(4,1)=A(4,1)/L(1,1) do i=1,n do k=1,i sum1=0 sum2=0 do j=1,k-1 if (i==k) then sum1=sum1+(L(k,j)*L(k,j)) L(k,k)=real(sqrt(A(k,k)-sum1)) elseif (i > k) then sum2=sum2+(L(i,j)*L(k,j)) L(i,k)=(1/L(k,k))*(A(i,k)-sum2) else L(i,k)=0 end if end do end doend do ! write outputdo i=1,m print "(3(1X,F6.1))",L(i,:)end do End program Cholesky_decomp Output:  5.0 0.0 0.0 3.0 3.0 0.0 -1.0 1.0 3.0  ## FreeBASIC Translation of: BBC BASIC ' version 18-01-2017' compile with: fbc -s console Sub Cholesky_decomp(array() As Double) Dim As Integer i, j, k Dim As Double s, l(UBound(array), UBound(array, 2)) For i = 0 To UBound(array) For j = 0 To i s = 0 For k = 0 To j -1 s += l(i, k) * l(j, k) Next If i = j Then l(i, j) = Sqr(array(i, i) - s) Else l(i, j) = (array(i, j) - s) / l(j, j) End If Next Next For i = 0 To UBound(array) For j = 0 To UBound(array, 2) Swap array(i, j), l(i, j) Next Next End Sub Sub Print_(array() As Double) Dim As Integer i, j For i = 0 To UBound(array) For j = 0 To UBound(array, 2) Print Using "###.#####";array(i,j); Next Print Next End Sub ' ------=< MAIN >=------ Dim m1(2,2) As Double => {{25, 15, -5}, _ {15, 18, 0}, _ {-5, 0, 11}} Dim m2(3, 3) As Double => {{18, 22, 54, 42}, _ {22, 70, 86, 62}, _ {54, 86, 174, 134}, _ {42, 62, 134, 106}} Cholesky_decomp(m1())Print_(m1()) PrintCholesky_decomp(m2())Print_(m2()) ' empty keyboard bufferWhile Inkey <> "" : WendPrint : Print "hit any key to end program"SleepEnd 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 ## Go ### Real This version works with real matrices, like most other solutions on the page. The representation is packed, however, storing only the lower triange of the input symetric matrix and the output lower matrix. The decomposition algorithm computes rows in order from top to bottom but is a little different thatn Cholesky–Banachiewicz. 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()} 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  ### Hermitian This version handles complex Hermitian matricies as described on the WP page. The matrix representation is flat, and storage is allocated for all elements, not just the lower triangles. The decomposition algorithm is Cholesky–Banachiewicz. package main import ( "fmt" "math/cmplx") type matrix struct { stride int ele []complex128} func like(a *matrix) *matrix { return &matrix{a.stride, make([]complex128, len(a.ele))}} func (m *matrix) print(heading string) { if heading > "" { fmt.Print("\n", heading, "\n") } for e := 0; e < len(m.ele); e += m.stride { fmt.Printf("%7.2f ", m.ele[e:e+m.stride]) fmt.Println() }} func (a *matrix) choleskyDecomp() *matrix { l := like(a) // Cholesky-Banachiewicz algorithm for r, rxc0 := 0, 0; r < a.stride; r++ { // calculate elements along row, up to diagonal x := rxc0 for c, cxc0 := 0, 0; c < r; c++ { sum := a.ele[x] for k := 0; k < c; k++ { sum -= l.ele[rxc0+k] * cmplx.Conj(l.ele[cxc0+k]) } l.ele[x] = sum / l.ele[cxc0+c] x++ cxc0 += a.stride } // calcualate diagonal element sum := a.ele[x] for k := 0; k < r; k++ { sum -= l.ele[rxc0+k] * cmplx.Conj(l.ele[rxc0+k]) } l.ele[x] = cmplx.Sqrt(sum) rxc0 += a.stride } return l} func main() { demo("A:", &matrix{3, []complex128{ 25, 15, -5, 15, 18, 0, -5, 0, 11, }}) demo("A:", &matrix{4, []complex128{ 18, 22, 54, 42, 22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106, }}) // one more example, from the Numpy manual, with a non-real demo("A:", &matrix{2, []complex128{ 1, -2i, 2i, 5, }})} func demo(heading string, a *matrix) { a.print(heading) a.choleskyDecomp().print("Cholesky factor L:")} Output: A: [( 25.00 +0.00i) ( 15.00 +0.00i) ( -5.00 +0.00i)] [( 15.00 +0.00i) ( 18.00 +0.00i) ( 0.00 +0.00i)] [( -5.00 +0.00i) ( 0.00 +0.00i) ( 11.00 +0.00i)] Cholesky factor L: [( 5.00 +0.00i) ( 0.00 +0.00i) ( 0.00 +0.00i)] [( 3.00 +0.00i) ( 3.00 +0.00i) ( 0.00 +0.00i)] [( -1.00 +0.00i) ( 1.00 +0.00i) ( 3.00 +0.00i)] A: [( 18.00 +0.00i) ( 22.00 +0.00i) ( 54.00 +0.00i) ( 42.00 +0.00i)] [( 22.00 +0.00i) ( 70.00 +0.00i) ( 86.00 +0.00i) ( 62.00 +0.00i)] [( 54.00 +0.00i) ( 86.00 +0.00i) ( 174.00 +0.00i) ( 134.00 +0.00i)] [( 42.00 +0.00i) ( 62.00 +0.00i) ( 134.00 +0.00i) ( 106.00 +0.00i)] Cholesky factor L: [( 4.24 +0.00i) ( 0.00 +0.00i) ( 0.00 +0.00i) ( 0.00 +0.00i)] [( 5.19 +0.00i) ( 6.57 +0.00i) ( 0.00 +0.00i) ( 0.00 +0.00i)] [( 12.73 +0.00i) ( 3.05 +0.00i) ( 1.65 +0.00i) ( 0.00 +0.00i)] [( 9.90 +0.00i) ( 1.62 +0.00i) ( 1.85 +0.00i) ( 1.39 +0.00i)] A: [( 1.00 +0.00i) ( 0.00 -2.00i)] [( 0.00 +2.00i) ( 5.00 +0.00i)] Cholesky factor L: [( 1.00 +0.00i) ( 0.00 +0.00i)] [( 0.00 +2.00i) ( 1.00 +0.00i)]  ### Library gonum/mat package main import ( "fmt" "gonum.org/v1/gonum/mat") func cholesky(order int, elements []float64) fmt.Formatter { var c mat.Cholesky c.Factorize(mat.NewSymDense(order, elements)) return mat.Formatted(c.LTo(nil))} func main() { fmt.Println(cholesky(3, []float64{ 25, 15, -5, 15, 18, 0, -5, 0, 11, })) fmt.Printf("\n%.5f\n", cholesky(4, []float64{ 18, 22, 54, 42, 22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106, }))} Output: ⎡ 5 0 0⎤ ⎢ 3 3 0⎥ ⎣-1 1 3⎦ ⎡ 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⎦  ### Library go.matrix package main import ( "fmt" mat "github.com/skelterjohn/go.matrix") func main() { demo(mat.MakeDenseMatrix([]float64{ 25, 15, -5, 15, 18, 0, -5, 0, 11, }, 3, 3)) demo(mat.MakeDenseMatrix([]float64{ 18, 22, 54, 42, 22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106, }, 4, 4))} func demo(m *mat.DenseMatrix) { fmt.Println("A:") fmt.Println(m) l, err := m.Cholesky() if err != nil { fmt.Println(err) return } fmt.Println("L:") fmt.Println(l)} Output: A: {25, 15, -5, 15, 18, 0, -5, 0, 11} L: { 5, 0, 0, 3, 3, 0, -1, 1, 3} A: { 18, 22, 54, 42, 22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106} L: { 4.242641, 0, 0, 0, 5.18545, 6.565905, 0, 0, 12.727922, 3.046038, 1.649742, 0, 9.899495, 1.624554, 1.849711, 1.392621}  ## Groovy Translation of: Java def decompose = { a -> assert a.size > 0 && a[0].size == a.size def m = a.size def l = [].withEagerDefault { [].withEagerDefault { 0 } } (0..<m).each { i -> (0..i).each { k -> Number s = (0..<k).sum { j -> l[i][j] * l[k][j] } ?: 0 l[i][k] = (i == k) ? Math.sqrt(a[i][i] - s) : (1.0 / l[k][k] * (a[i][k] - s)) } } l} Test: def test1 = [[25, 15, -5], [15, 18, 0], [-5, 0, 11]] def test2 = [[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]]; [test1,test2]. each { test -> println() decompose(test).each { println it[0..<(test.size)] }} Output: [5.0, 0, 0] [3.0, 3.0, 0] [-1.0, 1.0, 3.0] [4.242640687119285, 0, 0, 0] [5.185449728701349, 6.565905201197403, 0, 0] [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 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: module Cholesky (Arr, cholesky) where import Data.Array.IArrayimport Data.Array.MArrayimport Data.Array.Unboxedimport 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 -> Doubleget 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 -> Arrcholesky 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) The main module: import Data.Array.IArrayimport Data.Listimport Cholesky fm _ [] = "" fm _ [x] = fst xfm width ((a,b):xs) = a ++ (take (width - b)$ cycle " ") ++ (fm width xs) fmt width row (xs,[]) = fm width xsfmt width row (xs,ys) = fm width xs  ++ "\n" ++ fmt width row (splitAt row ys) showMatrice row xs   = ys where  vs = map (\s -> let sh = show s in (sh,length sh)) xs  width = (maximum $snd$ unzip vs) + 1  ys = fmt width row (splitAt row vs) ex1, ex2 :: Arrex1 = 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  putStrLn $showMatrice 3$ elems $cholesky ex1 putStrLn$ showMatrice 4 $elems$ cholesky ex2

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


### With Numeric.LinearAlgebra

import Numeric.LinearAlgebra a,b :: Matrix Ra = (3><3)   [25, 15, -5  ,15, 18, 0  ,-5,  0, 11] b = (4><4)  [ 18, 22, 54, 42  , 22, 70, 86, 62  , 54, 86,174,134  , 42, 62,134,106] main = do  let sa = sym a      sb = sym b  print sa  print $chol sa print sb print$ chol sb  print $tr$ chol sb
Output:
Herm (3><3)
[ 25.0, 15.0, -5.0
, 15.0, 18.0,  0.0
, -5.0,  0.0, 11.0 ]
(3><3)
[ 5.0, 3.0, -1.0
, 0.0, 3.0,  1.0
, 0.0, 0.0,  3.0 ]
Herm (4><4)
[ 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><4)
[ 4.242640687119285, 5.185449728701349, 12.727922061357857,  9.899494936611665
,               0.0, 6.565905201197403, 3.0460384954008553, 1.6245538642137891
,               0.0,               0.0, 1.6497422479090704,  1.849711005231382
,               0.0,               0.0,                0.0, 1.3926212476455904 ]
(4><4)
[  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.3926212476455904 ]

## Icon and Unicon

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 resultend procedure make_square_array (n)  result := []  every (1 to n) do push (result, list(n, 0))  return resultend 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
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


## Idris

works with Idris 0.10

Solution:

module Main import Data.Vect Matrix : Nat -> Nat -> Type -> TypeMatrix m n t = Vect m (Vect n t)  zeros : (m : Nat) -> (n : Nat) -> Matrix m n Doublezeros m n = replicate m (replicate n 0.0)  indexM : (Fin m, Fin n) -> Matrix m n t -> tindexM (i, j) a = index j (index i a)  replaceAtM : (Fin m, Fin n) -> t -> Matrix m n t -> Matrix m n treplaceAtM (i, j) e a = replaceAt i (replaceAt j e (index i a)) a  get : Matrix m m Double -> Matrix m m Double -> (Fin m, Fin m) -> Doubleget a l (i, j) {m} =  if i == j then sqrt $indexM (j, j) a - dot else if i > j then (indexM (i, j) a - dot) / indexM (j, j) l else 0.0 where -- Obtain indicies 0 to j -1 ks : List (Fin m) ks = case (findIndices (\_ => True) a) of [] => [] (x::xs) => init (x::xs) dot : Double dot = sum [(indexM (i, k) l) * (indexM (j, k) l) | k <- ks] updateL : Matrix m m Double -> Matrix m m Double -> (Fin m, Fin m) -> Matrix m m DoubleupdateL a l idx = replaceAtM idx (get a l idx) l cholesky : Matrix m m Double -> Matrix m m Doublecholesky a {m} = foldl (\l',i => foldl (\l'',j => updateL a l'' (i, j)) l' (js i)) l is where l = zeros m m is : List (Fin m) is = findIndices (\_ => True) a js : Fin m -> List (Fin m) js n = filter (<= n) is ex1 : Matrix 3 3 Doubleex1 = cholesky [[25.0, 15.0, -5.0], [15.0, 18.0, 0.0], [-5.0, 0.0, 11.0]] ex2 : Matrix 4 4 Doubleex2 = cholesky [[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]] main : IO () main = do print ex1 putStrLn "\n" print ex2 putStrLn "\n"  Output: [[5, 0, 0], [3, 3, 0], [-1, 1, 3]] [[4.242640687119285, 0, 0, 0], [5.185449728701349, 6.565905201197403, 0, 0], [12.72792206135786, 3.046038495400855, 1.64974224790907, 0], [9.899494936611665, 1.624553864213789, 1.849711005231382, 1.392621247645587]]  ## J Solution: mp=: +/ . * NB. matrix producth =: [email protected]|: 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. 'X Y t Z'=. , (;~n$(>.-:n){.1) <;.1 A  L0=. cholesky X  L1=. cholesky Z-(T=.(h Y) mp %.X) mp Y  L0,(T mp L0),.L1 end.)

See Cholesky Decomposition essay on the J Wiki.

Examples:
   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 eg24.24264       0       0       05.18545 6.56591       0       012.7279 3.04604 1.64974       09.89949 1.62455 1.84971 1.39262

Using math/lapack addon

   load 'math/lapack'   load 'math/lapack/potrf'   potrf_jlapack_ eg1 5 0 0 3 3 0_1 1 3   potrf_jlapack_ eg24.24264       0       0       05.18545 6.56591       0       012.7279 3.04604 1.64974       09.89949 1.62455 1.84971 1.39262

## Java

Works with: Java version 1.5+
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)));	}}
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]]

## JavaScript

 const cholesky = function (array) {	const zeros = [...Array(array.length)].map( _ => Array(array.length).fill(0));	const L = zeros.map((row, r, xL) => row.map((v, c) => {		const sum = row.reduce((s, _, i) => i < c ? s + xL[r][i] * xL[c][i] : s, 0);		return xL[r][c] = c < r + 1 ? r === c ? Math.sqrt(array[r][r] - sum) : (array[r][c] - sum) / xL[c][c] : v;	}));	return L;} let arr3 = [[25, 15, -5], [15, 18, 0], [-5, 0, 11]];console.log(cholesky(arr3));let arr4 = [[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]];console.log(cholesky(arr4));
Output:
0: (3) [5, 0, 0]
1: (3) [3, 3, 0]
2: (3) [-1, 1, 3]

0: (4) [4.242640687119285, 0, 0, 0]
1: (4) [5.185449728701349, 6.565905201197403, 0, 0]
2: (4) [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0]
3: (4) [9.899494936611665, 1.6245538642137891, 1.849711005231382, 1.3926212476455924]


## jq

Works with: jq version 1.4

Infrastructure:

# Create an m x n matrix def matrix(m; n; init):   if m == 0 then []   elif m == 1 then [range(0; n)] | map(init)   elif m > 0 then     matrix(1; n; init) as $row | [range(0; m)] | map($row )   else error("matrix\(m);_;_) invalid")   end ; # Print a matrix neatly, each cell ideally occupying n spaces,# but without truncationdef neatly(n):  def right: tostring | ( " " * (n-length) + .);  . as $in | length as$length  | reduce range (0; $length) as$i      (""; . + reduce range(0; $length) as$j      (""; "\(.) \($in[$i][$j] | right )" ) + "\n" ) ; def is_square: type == "array" and (map(type == "array") | all) and length == 0 or ( (.[0]|length) as$l | map(length == $l) | all) ; # This implementation of is_symmetric/0 uses a helper function that circumvents# limitations of jq 1.4:def is_symmetric: # [matrix, i,j, len] def test: if .[1] > .[3] then true elif .[1] == .[2] then [ .[0], .[1] + 1, 0, .[3]] | test elif .[0][.[1]][.[2]] == .[0][.[2]][.[1]] then [ .[0], .[1], .[2]+1, .[3]] | test else false end; if is_square|not then false else [ ., 0, 0, length ] | test end ;  Cholesky Decomposition: def cholesky_factor: if is_symmetric then length as$length    | . as $self | reduce range(0;$length) as $k ( matrix(length; length; 0); # the matrix that will hold the answer reduce range(0;$length) as $i (.; if$i == $k then (. as$lower                     | reduce range(0; $k) as$j                         (0; . + ($lower[$k][$j] | .*.) )) as$sum                 | .[$k][$k] = (($self[$k][$k] -$sum) | sqrt)             elif $i >$k               then (. as $lower | reduce range(0;$k) as $j (0; . +$lower[$i][$j] * $lower[$k][$j])) as$sum                 | .[$i][$k] = (($self[$k][$i] -$sum) / .[$k][$k] )             else .             end ))  else error( "cholesky_factor: matrix is not symmetric" )  end ;

[[25,15,-5],[15,18,0],[-5,0,11]] | cholesky_factor

Output:
[[5,0,0],[3,3,0],[-1,1,3]]


[[18, 22,  54,  42],
[22, 70,  86,  62],
[54, 86, 174, 134],
[42, 62, 134, 106]] | cholesky_factor | neatly(20)

Output:
    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

## Julia

Julia's strong linear algebra support includes Cholesky decomposition.

 a = [25 15 5; 15 18 0; -5 0 11]b = [18 22 54 22; 22 70 86 62; 54 86 174 134; 42 62 134 106] println(a, "\n => \n", chol(a, :L))println(b, "\n => \n", chol(b, :L))
Output:
[25 15 5
15 18 0
-5 0 11]
=>
[5.0 0.0 0.0
3.0 3.0 0.0
-1.0 1.0 3.0]
[18 22 54 22
22 70 86 62
54 86 174 134
42 62 134 106]
=>
[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]


## Kotlin

Translation of: C
// version 1.0.6 fun cholesky(a: DoubleArray): DoubleArray {    val n = Math.sqrt(a.size.toDouble()).toInt()    val l = DoubleArray(a.size)     var s: Double    for (i in 0 until n)         for (j in 0 .. i) {            s = 0.0            for (k in 0 until j) s += l[i * n + k] * l[j * n + k]            l[i * n + j] = when {                (i == j) -> Math.sqrt(a[i * n + i] - s)                else     -> 1.0 / l[j * n + j] * (a[i * n + j] - s)            }        }    return l} fun showMatrix(a: DoubleArray) {    val n = Math.sqrt(a.size.toDouble()).toInt()    for (i in 0 until n) {        for (j in 0 until n) print("%8.5f ".format(a[i * n + j]))        println()    }}  fun main(args: Array<String>) {    val m1 = doubleArrayOf(25.0, 15.0, -5.0,                           15.0, 18.0,  0.0,                           -5.0,  0.0, 11.0)    val c1 = cholesky(m1)    showMatrix(c1)    println()    val m2 = doubleArrayOf(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)    val c2 = cholesky(m2)    showMatrix(c2) }
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


## Lobster

Translation of: Go

Translated from the Go Real version: This version works with real matrices, like most other solutions on the page. The representation is packed, however, storing only the lower triange of the input symetric matrix and the output lower matrix. The decomposition algorithm computes rows in order from top to bottom but is a little different than Cholesky–Banachiewicz.

import std // choleskyLower returns the cholesky decomposition of a symmetric real// matrix.  The matrix must be positive definite but this is not checkeddef choleskyLower(order, a) -> [float]:    let l = map(a.length): 0.0    var row, col = 1, 1    var dr = 0 // index of diagonal element at end of row    var dc = 0 // index of diagonal element at top of column    for(a) e, i:        if i < dr:            let d = (e - l[i]) / l[dc]            l[i] = d            var ci, cx = col, dc            var j = i + 1            while j <= dr:                cx += ci                ci += 1                l[j] += d * l[cx]                j += 1            col += 1            dc += col        else:            l[i] = sqrt(e - l[i])            row += 1            dr += row            col = 1            dc = 0    return l // symmetric.print prints a square matrix from the packed representation,// printing the upper triange as a transpose of the lowerdef print_symmetric(order, s):    //const eleFmt = "%10.5f "    var str = ""    var row, diag = 1, 0    for(s) e, i:        str += e + " " // format?        if i == diag:            var j, col = diag+row, row            while col < order:                str += s[j] + " " // format?                col++                j += col            print(str); str = ""            row += 1            diag += row // lower.print prints a square matrix from the packed representation,// printing the upper triangle as all zeros.def print_lower(order, l):    //const eleFmt = "%10.5f "    var str = ""    var row, diag = 1, 0    for(l) e, i:        str += e + " " // format?        if i == diag:            var j = row            while j < order:                str += 0.0 + " " // format?                j += 1            print(str); str = ""            row += 1            diag += row def demo(order, a):    print("A:")    print_symmetric(order, a)    print("L:")    print_lower(order, choleskyLower(order, a)) demo(3, [25.0,         15.0, 18.0,         -5.0, 0.0, 11.0]) demo(4, [18.0,         22.0, 70.0,         54.0, 86.0, 174.0,         42.0, 62.0, 134.0, 106.0])
Output:
A:
25.0 15.0 -5.0
15.0 18.0 0.0
-5.0 0.0 11.0
L:
5.0 0.0 0.0
3.0 3.0 0.0
-1.0 1.0 3.0
A:
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:
4.242640687119 0.0 0.0 0.0
5.185449728701 6.565905201197 0.0 0.0
12.72792206135 3.046038495401 1.649742247909 0.0
9.899494936612 1.624553864214 1.849711005231 1.392621247646


## Maple

The Cholesky decomposition is obtained by passing the method = Cholesky' option to the LUDecomposition procedure in the LinearAlgebra pacakge. This is illustrated below for the two requested examples. The first is computed exactly; the second is also, but the subsequent application of evalf' to the result produces a matrix with floating point entries which can be compared with the expected output in the problem statement.

> A := << 25, 15, -5; 15, 18, 0; -5, 0, 11 >>;                              [25    15    -5]                              [              ]                         A := [15    18     0]                              [              ]                              [-5     0    11] > B := << 18, 22, 54, 42; 22, 70, 86, 62; 54, 86, 174, 134; 42, 62, 134, 106>>;                          [18    22     54     42]                          [                      ]                          [22    70     86     62]                     B := [                      ]                          [54    86    174    134]                          [                      ]                          [42    62    134    106] > use LinearAlgebra in>       LUDecomposition( A, method = Cholesky );>       LUDecomposition( B, method = Cholesky );>       evalf( % );> end use;                             [ 5    0    0]                             [            ]                             [ 3    3    0]                             [            ]                             [-1    1    3]              [   1/2                                      ]             [3 2           0            0            0   ]             [                                            ]             [    1/2        1/2                          ]             [11 2       2 97                             ]             [-------    -------         0            0   ]             [   3          3                             ]             [                                            ]             [                1/2          1/2            ]             [   1/2     30 97       2 6402               ]             [9 2        --------    ---------        0   ]             [              97          97                ]             [                                            ]             [                1/2           1/2        1/2]             [   1/2     16 97       74 6402       8 33   ]             [7 2        --------    ----------    -------]             [              97          3201         33   ]        [4.242640686        0.             0.             0.     ]       [                                                        ]       [5.185449728    6.565905202        0.             0.     ]       [                                                        ]       [12.72792206    3.046038495    1.649742248        0.     ]       [                                                        ]       [9.899494934    1.624553864    1.849711006    1.392621248]

## Mathematica / Wolfram Language

CholeskyDecomposition[{{25, 15, -5}, {15, 18, 0}, {-5, 0, 11}}]

Without the use of built-in functions, making use of memoization:

chol[A_] := Module[{L},  L[k_, k_] := L[k, k] = Sqrt[A[[k, k]] - Sum[L[k, j]^2, {j, 1, k-1}]];  L[i_, k_] := L[i, k] = L[k, k]^-1 (A[[i, k]] - Sum[L[i, j] L[k, j], {j, 1, k-1}]);  PadRight[Table[L[i, j], {i, Length[A]}, {j, i}]] ]

## MATLAB / Octave

The cholesky decomposition chol() is an internal function

  A = [  25  15  -5  15  18   0  -5   0  11 ];   B  = [   18  22   54   42  22  70   86   62  54  86  174  134  42  62  134  106   ];   [L] = chol(A,'lower')    [L] = chol(B,'lower')
Output:
   >   [L] = chol(A,'lower')
L =

5   0   0
3   3   0
-1   1   3

>   [L] = chol(B,'lower')
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


/* Cholesky decomposition is built-in */ a: hilbert_matrix(4)$b: cholesky(a);/* matrix([1, 0, 0, 0 ], [1/2, 1/(2*sqrt(3)), 0, 0 ], [1/3, 1/(2*sqrt(3)), 1/(6*sqrt(5)), 0 ], [1/4, 3^(3/2)/20, 1/(4*sqrt(5)), 1/(20*sqrt(7))]) */ b . transpose(b) - a;matrix([0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]) ## Nim Translation of: C import math, strutils, strformat type Matrix[N: static int, T: SomeFloat] = array[N, array[N, T]] proc cholesky[Matrix](a: Matrix): Matrix = for i in 0 ..< a[0].len: for j in 0 .. i: var s = 0.0 for k in 0 ..< j: s += result[i][k] * result[j][k] result[i][j] = if i == j: sqrt(a[i][i]-s) else: 1.0 / result[j][j] * (a[i][j] - s) proc $(a: Matrix): string =  result = ""  for b in a:    var line = ""    for c in b:      line.addSep(" ", 0)      line.add fmt"{c:8.5f}"    result.add line & '\n' let m1 = [[25.0, 15.0, -5.0],          [15.0, 18.0,  0.0],          [-5.0,  0.0, 11.0]]echo cholesky(m1) let m2 = [[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]]echo cholesky(m2)
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

## Objeck

Translation of: C
 class Cholesky {  function : Main(args : String[]) ~ Nil {    n := 3;    m1 := [25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0];    c1 := Cholesky(m1, n);    ShowMatrix(c1, n);     IO.Console->PrintLine();     n := 4;    m2 := [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];    c2 := Cholesky(m2, n);    ShowMatrix(c2, n);  }   function : ShowMatrix(A : Float[], n : Int) ~ Nil {    for (i := 0; i < n; i+=1;) {      for (j := 0; j < n; j+=1;) {        IO.Console->Print(A[i * n + j])->Print('\t');      };      IO.Console->PrintLine();    };  }   function : Cholesky(A : Float[], n : Int) ~ Float[] {    L := Float->New[n * n];     for (i := 0; i < n; i+=1;) {      for (j := 0; j < (i+1); j+=1;) {        s := 0.0;        for (k := 0; k < j; k+=1;) {          s += L[i * n + k] * L[j * n + k];        };        L[i * n + j] := (i = j) ?          (A[i * n + i] - s)->SquareRoot() :          (1.0 / L[j * n + j] * (A[i * n + j] - s));      };    };     return L;  }}
5       0       0
3       3       0
-1      1       3

4.24264069      0               0               0
5.18544973      6.5659052       0               0
12.7279221      3.0460385       1.64974225      0
9.89949494      1.62455386      1.84971101      1.39262125


## 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_veclet 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|] |];
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

## ooRexx

Translation of: REXX
/*REXX program performs the  Cholesky  decomposition  on a square matrix.     */niner =  '25  15  -5' ,                              /*define a  3x3  matrix. */         '15  18   0' ,         '-5   0  11'                           call Cholesky ninerhexer =  18  22  54  42,                             /*define a  4x4  matrix. */         22  70  86  62,         54  86 174 134,         42  62 134 106                           call Cholesky hexerexit                                   /*stick a fork in it,  we're all done. *//*----------------------------------------------------------------------------*/Cholesky: procedure;  parse arg mat;   say;   say;   call tell 'input matrix',mat             do    r=1  for ord                do c=1  for r; d=0;  do i=1  for c-1; d=d+!.r.i*!.c.i; end /*i*/                if r=c  then !.r.r=sqrt(!.r.r-d)                        else !.r.c=1/!.c.c*(a.r.c-d)                end   /*c*/             end      /*r*/          call tell  'Cholesky factor',,!.,'-'          return/*----------------------------------------------------------------------------*/err:   say;  say;  say '***error***!';    say;  say arg(1);  say;  say;  exit 13/*----------------------------------------------------------------------------*/tell:  parse arg hdr,x,y,sep;   n=0;             if sep==''  then sep='-'       dPlaces= 5                    /*n decimal places past the decimal point*/       width  =10                    /*width of field used to display elements*/       if y==''  then !.=0                 else do row=1  for ord; do col=1  for ord; x=x !.row.col; end; end       w=words(x)           do ord=1  until ord**2>=w; end  /*a fast way to find matrix's order*/       say       if ord**2\==w  then call err  "matrix elements don't form a square matrix."       say center(hdr, ((width+1)*w)%ord, sep)       say               do   row=1  for ord;       z=''                 do col=1  for ord;       n=n+1                                    a.row.col=word(x,n)                 if col<=row  then  !.row.col=a.row.col                 z=z  right( format(a.row.col,, dPlaces) / 1,   width)                 end   /*col*/               say z               end        /*row*/       return/*----------------------------------------------------------------------------*/sqrt:  procedure; parse arg x;   if x=0  then return 0;  d=digits();  i=''; m.=9       numeric digits 9; numeric form; h=d+6; if x<0  then  do; x=-x; i='i'; end       parse value format(x,2,1,,0) 'E0'  with  g 'E' _ .;       g=g*.5'e'_%2          do j=0  while h>9;      m.j=h;              h=h%2+1;        end  /*j*/          do k=j+5  to 0  by -1;  numeric digits m.k; g=(g+x/g)*.5;   end  /*k*/       numeric digits d;     return (g/1)i            /*make complex if X < 0.*/

## PARI/GP

cholesky(M) ={  my (L = matrix(#M,#M));   for (i = 1, #M,    for (j = 1, i,      s = sum (k = 1, j-1, L[i,k] * L[j,k]);      L[i,j] = if (i == j, sqrt(M[i,i] - s), (M[i,j] - s) / L[j,j])    )  );  L}

Output: (set displayed digits with: \p 5)

gp > cholesky([25,15,-5;15,18,0;-5,0,11])

[ 5.0000      0      0]

[ 3.0000 3.0000      0]

[-1.0000 1.0000 3.0000]

gp > cholesky([18,22,54,42;22,70,86,62;54,86,174,134;42,62,134,106])

[4.2426      0      0      0]

[5.1854 6.5659      0      0]

[12.728 3.0460 1.6497      0]

[9.8995 1.6246 1.8497 1.3926]


## Pascal

program CholeskyApp; type  D2Array = array of array of double; function cholesky(const A: D2Array): D2Array;var  i, j, k: integer;  s: double;begin  setlength(Result, length(A), length(A));  for i := low(Result) to high(Result) do    for j := 0 to i do    begin      s := 0;      for k := 0 to j - 1 do        s := s + Result[i][k] * Result[j][k];      if i = j then        Result[i][j] := sqrt(A[i][i] - s)      else        Result[i][j] := (A[i][j] - s) / Result[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, i: integer;  cIn, cOut: D2Array; begin  setlength(cIn, length(m1), length(m1));  for index := low(m1) to high(m1) do  begin    SetLength(cIn[index], length(m1[index]));    for i := 0 to High(m1[Index]) do      cIn[index][i] := m1[index][i];  end;  cOut := cholesky(cIn);  printM(cOut);   writeln;   setlength(cIn, length(m2), length(m2));  for index := low(m2) to high(m2) do  begin    SetLength(cIn[index], length(m2[Index]));    for i := 0 to High(m2[Index]) do      cIn[index][i] := m2[index][i];  end;  cOut := cholesky(cIn);  printM(cOut);end.
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

sub cholesky {   my $matrix = shift; my$chol = [ map { [(0) x @$matrix ] } @$matrix ];   for my $row (0..@$matrix-1) {     for my $col (0..$row) {       my $x = $$matrix[row][col]; x -=$$chol[$row][$_]*$$chol[col][_] for 0..col;$$chol[$row][$col] =$row == $col ? sqrt$x : $x/$$chol[$col][$col]; } } return$chol; }  my $example1 = [ [ 25, 15, -5 ], [ 15, 18, 0 ], [ -5, 0, 11 ] ]; print "Example 1:\n"; print +(map { sprintf "%7.4f\t",$_ } @$_), "\n" for @{ cholesky$example1 };  my $example2 = [ [ 18, 22, 54, 42], [ 22, 70, 86, 62], [ 54, 86, 174, 134], [ 42, 62, 134, 106] ]; print "\nExample 2:\n"; print +(map { sprintf "%7.4f\t",$_ } @$_), "\n" for @{ cholesky$example2 };
Output:
Example 1:
5.0000	 0.0000	 0.0000
3.0000	 3.0000	 0.0000
-1.0000	 1.0000	 3.0000

Example 2:
4.2426	 0.0000	 0.0000	 0.0000
5.1854	 6.5659	 0.0000	 0.0000
12.7279	 3.0460	 1.6497	 0.0000
9.8995	 1.6246	 1.8497	 1.3926


## Phix

Translation of: Sidef
with javascript_semantics
function cholesky(sequence matrix)
integer l = length(matrix)
sequence chol = repeat(repeat(0,l),l)
for row=1 to l do
for col=1 to row do
atom x = matrix[row][col]
for i=1 to col do
x -= chol[row][i] * chol[col][i]
end for
chol[row][col] = iff(row == col ? sqrt(x) : x/chol[col][col])
end for
end for
return chol
end function

ppOpt({pp_Nest,1})
pp(cholesky({{ 25, 15, -5 },
{ 15, 18,  0 },
{ -5,  0, 11 }}))
pp(cholesky({{ 18, 22,  54,  42},
{ 22, 70,  86,  62},
{ 54, 86, 174, 134},
{ 42, 62, 134, 106}}))

Output:
{{5,0,0},
{3,3,0},
{-1,1,3}}
{{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}}


## 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 S 1.0)                     (*/ S 1.0 (get L J J)) ) ) ) ) )      (for R L         (for N R (prin (align 9 (round N 5))))         (prinl) ) ) )

Test:

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

## PL/I

(subscriptrange):decompose: procedure options (main);   /* 31 October 2013 */   declare a(*,*) float controlled;    allocate a(3,3) initial (25, 15, -5,                            15, 18,  0,                            -5,  0, 11);    put skip list ('Original matrix:');    put edit (a) (skip, 3 f(4));     call cholesky(a);    put skip list ('Decomposed matrix');    put edit (a) (skip, 3 f(4));    free a;    allocate a(4,4) initial (18, 22,  54,  42,                             22, 70,  86,  62,                             54, 86, 174, 134,                             42, 62, 134, 106);    put skip list ('Original matrix:');    put edit (a) (skip, (hbound(a,1)) f(12) );    call cholesky(a);    put skip list ('Decomposed matrix');    put edit (a) (skip, (hbound(a,1)) f(12,5) ); cholesky: procedure(a);   declare a(*,*) float;   declare L(hbound(a,1), hbound(a,2)) float;   declare s float;   declare (i, j, k) fixed binary;    L = 0;   do i = lbound(a,1) to hbound(a,1);      do j = lbound(a,2) to i;         s = 0;         do k = lbound(a,2) to j-1;            s = s + L(i,k) * L(j,k);         end;         if i = j then            L(i,j) = sqrt(a(i,i) - s);         else            L(i,j) = (a(i,j) - s) / L(j,j);      end;   end;   a = L;end cholesky; end decompose;

ACTUAL RESULTS:-

Original matrix:
25  15  -5
15  18   0
-5   0  11
Decomposed matrix
5   0   0
3   3   0
-1   1   3
Original matrix:
18          22          54          42
22          70          86          62
54          86         174         134
42          62         134         106
Decomposed matrix
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


## PowerShell

 function cholesky ($a) {$l = @()    if ($a) {$n = $a.count$end = $n - 1$l = 1..$n | foreach {$row = @(0) * $n; ,$row}        foreach ($k in 0..$end) {            $m =$k - 1            $sum = 0 if(0 -lt$k) {                foreach ($j in 0..$m) {$sum +=$l[$k][$j]*$l[$k][$j]} }$l[$k][$k] = [Math]::Sqrt($a[$k][$k] -$sum)            if ($k -lt$end) {                foreach ($i in ($k+1)..$end) {$sum = 0                    if (0 -lt $k) { foreach ($j in 0..$m) {$sum += $l[$i][$j]*$l[$k][$j]}                    }                    $l[$i][$k] = ($a[$i][$k] - $sum)/$l[$k][$k]                }            }        }    }    $l} function show($a) {$a | foreach {"$_"}} $a1 = @(@(25, 15, -5),@(15, 18, 0),@(-5, 0, 11))"a1 ="show$a1"""l1 ="show (cholesky $a1)""$a2 = @(@(18, 22, 54, 42),@(22, 70, 86, 62),@(54, 86, 174, 134),@(42, 62, 134, 106))"a2 ="show $a2"""l2 ="show (cholesky$a2)

Output:

a1 =
25 15 -5
15 18 0
-5 0 11

l1 =
5 0 0
3 3 0
-1 1 3

a2 =
18 22 54 42
22 70 86 62
54 86 174 134
42 62 134 106

l2 =
4.24264068711928 0 0 0
5.18544972870135 6.5659052011974 0 0
12.7279220613579 3.04603849540086 1.64974224790907 0
9.89949493661167 1.62455386421379 1.84971100523138 1.39262124764559


## Python

### Python2.X version

from __future__ import print_function from pprint import pprintfrom math import sqrt  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] = sqrt(A[i][i] - s) if (i == j) else \                      (1.0 / L[j][j] * (A[i][j] - s))    return L if __name__ == "__main__":    m1 = [[25, 15, -5],          [15, 18,  0],          [-5,  0, 11]]    pprint(cholesky(m1))    print()     m2 = [[18, 22,  54,  42],          [22, 70,  86,  62],          [54, 86, 174, 134],          [42, 62, 134, 106]]    pprint(cholesky(m2), width=120)
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]]


### Python3.X version using extra Python idioms

Factors out accesses to A[i], L[i], and L[j] by creating Ai, Li and Lj respectively as well as using enumerate instead of range(len(some_array)).

def cholesky(A):    L = [[0.0] * len(A) for _ in range(len(A))]    for i, (Ai, Li) in enumerate(zip(A, L)):        for j, Lj in enumerate(L[:i+1]):            s = sum(Li[k] * Lj[k] for k in range(j))            Li[j] = sqrt(Ai[i] - s) if (i == j) else \                      (1.0 / Lj[j] * (Ai[j] - s))    return L
Output:

(As above)

## q

solve:{[A;B] $[0h>type A;B%A;inv[A] mmu B]}ak:{[m;k] (),/:m[;k]til k:k-1}akk:{[m;k] m[k;k:k-1]}transpose:{$[0h=type x;flip x;enlist each x]}mult:{[A;B]$[0h=type A;A mmu B;A*B]} cholesky:{[A] {[A;L;n] l_k:solve[L;ak[A;n]]; l_kk:first over sqrt[akk[A;n] - mult[transpose l_k;l_k]]; ({$[0h<type x;enlist x;x]}L,'0f),enlist raze transpose[l_k],l_kk		}[A]/[sqrt A[0;0];1_1+til count first A]	} show cholesky (25 15 -5f;15 18 0f;-5 0 11f)-1"";show cholesky (18 22 54 42f;22 70 86 62f;54 86 174 134f;42 62 134 106f)
Output:
5  0 0
3  3 0
-1 1 3

4.242641 0        0        0
5.18545  6.565905 0        0
12.72792 3.046038 1.649742 0
9.899495 1.624554 1.849711 1.392621


## R

t(chol(matrix(c(25, 15, -5, 15, 18, 0, -5, 0, 11), nrow=3, ncol=3)))#      [,1] [,2] [,3]# [1,]    5    0    0# [2,]    3    3    0# [3,]   -1    1    3 t(chol(matrix(c(18, 22, 54, 42, 22, 70, 86, 62, 54, 86, 174, 134, 42, 62, 134, 106), nrow=4, ncol=4)))#           [,1]     [,2]     [,3]     [,4]# [1,]  4.242641 0.000000 0.000000 0.000000# [2,]  5.185450 6.565905 0.000000 0.000000# [3,] 12.727922 3.046038 1.649742 0.000000# [4,]  9.899495 1.624554 1.849711 1.392621

## Racket

 #lang racket(require math) (define (cholesky A)  (define mref matrix-ref)  (define n (matrix-num-rows A))  (define L (for/vector ([_ n]) (for/vector ([_ n]) 0)))  (define (set L i j x) (vector-set! (vector-ref L i) j x))  (define (ref L i j) (vector-ref (vector-ref L i) j))  (for* ([i n] [k n])    (set L i k         (cond            [(= i k)             (sqrt (- (mref A i i) (for/sum ([j k]) (sqr (ref L k j)))))]           [(> i k)             (/ (- (mref A i k) (for/sum ([j k]) (* (ref L i j) (ref L k j))))               (ref L k k))]           [else 0])))  L) (cholesky (matrix [[25 15 -5]                   [15 18  0]                   [-5  0 11]])) (cholesky (matrix [[18 22  54 42]                   [22 70  86 62]                   [54 86 174 134]                   [42 62 134 106]]))

Output:

 '#(#(5 0 0)    #(3 3 0)   #(-1 1 3))'#(#(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))

## Raku

(formerly Perl 6)

Works with: Rakudo version 2015.12
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;*])[^$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],];

## REXX

If trailing zeros are wanted after the decimal point for values of zero (0),   the     / 1     (a division by unity performs
REXX number normalization)   can be removed from the line   (number 40)   which contains the source statement:

z=z   right( format(@.row.col, ,   dPlaces) / 1,     width)
/*REXX program performs the Cholesky decomposition on a square matrix & displays results*/niner =  '25  15  -5' ,                          /*define a  3x3  matrix with elements. */         '15  18   0' ,         '-5   0  11'                           call Cholesky ninerhexer =  18  22  54  42,                         /*define a  4x4  matrix with elements. */         22  70  86  62,         54  86 174 134,         42  62 134 106                           call Cholesky hexerexit                                             /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/Cholesky: procedure;  parse arg mat;   say;   say;   call tell 'input array',mat              do    r=1  for ord                 do c=1  for r; $=0; do i=1 for c-1;$= $+ !.r.i * !.c.i; end /*i*/ if r=c then !.r.r= sqrt(!.r.r -$)                          else !.r.c= 1 / !.c.c * (@.r.c - $) end /*c*/ end /*r*/ call tell 'Cholesky factor',,!.,'─' return/*──────────────────────────────────────────────────────────────────────────────────────*/err: say; say; say '***error***!'; say; say arg(1); say; say; exit 13/*──────────────────────────────────────────────────────────────────────────────────────*/tell: parse arg hdr,x,y,sep; #=0; if sep=='' then sep= '═' dPlaces= 5 /*# dec. places past the decimal point.*/ width =10 /*field width used to display elements.*/ if y=='' then !.=0 else do row=1 for ord; do col=1 for ord; x=x !.row.col; end; end w=words(x) do ord=1 until ord**2>=w; end /*a fast way to find the matrix's order*/ say if ord**2\==w then call err "matrix elements don't form a square matrix." say center(hdr, ((width + 1) * w) % ord, sep) say do row=1 for ord; z= do col=1 for ord; #= # + 1 @.row.col= word(x, #) if col<=row then !.row.col= @.row.col z=z right( format(@.row.col, , dPlaces) / 1, width) end /*col*/ /* ↑↑↑ */ say z /* └┴┴──◄──normalization for zero*/ end /*row*/ return/*──────────────────────────────────────────────────────────────────────────────────────*/sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6 numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_ %2 do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/; return g/1 output: ═══════════input matrix══════════ 25 15 -5 15 18 0 -5 0 11 ─────────Cholesky factor───────── 5 0 0 3 3 0 -1 1 3 ════════════════input matrix════════════════ 18 22 54 42 22 70 86 62 54 86 174 134 42 62 134 106 ──────────────Cholesky factor─────────────── 4.24264 0 0 0 5.18545 6.56591 0 0 12.72792 3.04604 1.64974 0 9.89949 1.62455 1.84971 1.39262  ## Ring  # Project : Cholesky decomposition load "stdlib.ring"decimals(5)m1 = [[25, 15, -5], [15, 18, 0], [-5, 0, 11]]cholesky(m1)printarray(m1)see nl m2 = [[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]]cholesky(m2)printarray(m2) func cholesky(a)l = newlist(len(a), len(a))for i = 1 to len(a) for j = 1 to i s = 0 for k = 1 to j s = s + l[i][k] * l[j][k] next if i = j l[i][j] = sqrt(a[i][i] - s) else l[i][j] = (a[i][j] - s) / l[j][j] ok next next a = l func printarray(a) for row = 1 to len(a) for col = 1 to len(a) see "" + a[row][col] + " " next see nl next  Output: 5 0 0 3 3 0 -1 1 3 4.24264 0 0 0 5.18545 6.56591 0 0 12.72792 3.04604 1.64974 0 9.89949 1.62455 1.84971 1.39262  ## 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] endend puts Matrix[[25,15,-5],[15,18,0],[-5,0,11]].cholesky_factorputs Matrix[[18, 22, 54, 42], [22, 70, 86, 62], [54, 86, 174, 134], [42, 62, 134, 106]].cholesky_factor Output: 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]]  ## Rust Translation of: C fn cholesky(mat: Vec<f64>, n: usize) -> Vec<f64> { let mut res = vec![0.0; mat.len()]; for i in 0..n { for j in 0..(i+1){ let mut s = 0.0; for k in 0..j { s += res[i * n + k] * res[j * n + k]; } res[i * n + j] = if i == j { (mat[i * n + i] - s).sqrt() } else { (1.0 / res[j * n + j] * (mat[i * n + j] - s)) }; } } res} fn show_matrix(matrix: Vec<f64>, n: usize){ for i in 0..n { for j in 0..n { print!("{:.4}\t", matrix[i * n + j]); } println!(""); } println!("");} fn main(){ let dimension = 3 as usize; let m1 = vec![25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0]; let res1 = cholesky(m1, dimension); show_matrix(res1, dimension); let dimension = 4 as usize; let m2 = vec![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]; let res2 = cholesky(m2, dimension); show_matrix(res2, dimension);}  Output: 5.0000 0.0000 0.0000 3.0000 3.0000 0.0000 -1.0000 1.0000 3.0000 4.2426 0.0000 0.0000 0.0000 5.1854 6.5659 0.0000 0.0000 12.7279 3.0460 1.6497 0.0000 9.8995 1.6246 1.8497 1.3926  ## Scala case class Matrix( val matrix:Array[Array[Double]] ) { // Assuming matrix is positive-definite, symmetric and not empty... val rows,cols = matrix.size def getOption( r:Int, c:Int ) : Option[Double] = Pair(r,c) match { case (r,c) if r < rows && c < rows => Some(matrix(r)(c)) case _ => None } def isLowerTriangle( r:Int, c:Int ) : Boolean = { c <= r } def isDiagonal( r:Int, c:Int ) : Boolean = { r == c} override def toString = matrix.map(_.mkString(", ")).mkString("\n") /** * Perform Cholesky Decomposition of this matrix */ lazy val cholesky : Matrix = { val l = Array.ofDim[Double](rows*cols) for( i <- (0 until rows); j <- (0 until cols) ) yield { val s = (for( k <- (0 until j) ) yield { l(i*rows+k) * l(j*rows+k) }).sum l(i*rows+j) = (i,j) match { case (r,c) if isDiagonal(r,c) => scala.math.sqrt(matrix(i)(i) - s) case (r,c) if isLowerTriangle(r,c) => (1.0 / l(j*rows+j) * (matrix(i)(j) - s)) case _ => 0 } } val m = Array.ofDim[Double](rows,cols) for( i <- (0 until rows); j <- (0 until cols) ) m(i)(j) = l(i*rows+j) Matrix(m) }} // A little test...val a1 = Matrix(Array[Array[Double]](Array(25,15,-5),Array(15,18,0),Array(-5,0,11)))val a2 = Matrix(Array[Array[Double]](Array(18,22,54,42), Array(22,70,86,62), Array(54,86,174,134), Array(42,62,134,106))) val l1 = a1.choleskyval l2 = a2.cholesky // Given test resultsval r1 = Array[Double](5,0,0,3,3,0,-1,1,3)val r2 = Array[Double](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) // Verify assertions (l1.matrix.flatten.zip(r1)).foreach{ case (result,test) => assert(math.round( result * 100000 ) * 0.00001 == math.round( test * 100000 ) * 0.00001) } (l2.matrix.flatten.zip(r2)).foreach{ case (result,test) => assert(math.round( result * 100000 ) * 0.00001 == math.round( test * 100000 ) * 0.00001) } ## Scilab The Cholesky decomposition is builtin, and an upper triangular matrix is returned, such that$A=L^TL$. a = [25 15 -5; 15 18 0; -5 0 11];chol(a) ans = 5. 3. -1. 0. 3. 1. 0. 0. 3. a = [18 22 54 42; 22 70 86 62; 54 86 174 134; 42 62 134 106]; chol(a) ans = 4.2426407 5.1854497 12.727922 9.8994949 0. 6.5659052 3.0460385 1.6245539 0. 0. 1.6497422 1.849711 0. 0. 0. 1.3926212 ## Seed7 $ include "seed7_05.s7i";  include "float.s7i";  include "math.s7i"; const type: matrix is array array float; const func matrix: cholesky (in matrix: a) is func  result    var matrix: cholesky is 0 times 0 times 0.0;  local    var integer: i is 0;    var integer: j is 0;    var integer: k is 0;    var float: sum is 0.0;  begin    cholesky := length(a) times length(a) times 0.0;    for key i range cholesky do      for j range 1 to i do	sum := 0.0;	for k range 1 to j do	  sum +:= cholesky[i][k] * cholesky[j][k];        end for;	if i = j then	  cholesky[i][i] := sqrt(a[i][i] - sum)	else          cholesky[i][j] := (a[i][j] - sum) / cholesky[j][j];        end if;      end for;    end for;  end func;   const proc: writeMat (in matrix: a) is func  local    var integer: i is 0;    var float: num is 0.0;  begin    for key i range a do      for num range a[i] do        write(num digits 5 lpad 8);      end for;      writeln;    end for;  end func; const matrix: m1 is [] (    [] (25.0, 15.0, -5.0),    [] (15.0, 18.0,  0.0),    [] (-5.0,  0.0, 11.0));const matrix: m2 is [] (    [] (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)); const proc: main is func  begin    writeMat(cholesky(m1));    writeln;    writeMat(cholesky(m2));  end func;

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.89950 1.62455 1.84971 1.39262


## Sidef

Translation of: Perl
func cholesky(matrix) {    var chol = matrix.len.of { matrix.len.of(0) }    for row in ^matrix {        for col in (0..row) {            var x = matrix[row][col]            for i in (0..col) {                x -= (chol[row][i] * chol[col][i])            }            chol[row][col] = (row == col ? x.sqrt : x/chol[col][col])        }    }    return chol}

Examples:

var example1 = [ [ 25, 15, -5 ],                 [ 15, 18,  0 ],                 [ -5,  0, 11 ] ]; say "Example 1:";cholesky(example1).each { |row|    say row.map {'%7.4f' % _}.join(' ');} var example2 = [ [ 18, 22,  54,  42],                 [ 22, 70,  86,  62],                 [ 54, 86, 174, 134],                 [ 42, 62, 134, 106] ]; say "\nExample 2:";cholesky(example2).each { |row|    say row.map {'%7.4f' % _}.join(' ');}
Output:
Example 1:
5.0000  0.0000  0.0000
3.0000  3.0000  0.0000
-1.0000  1.0000  3.0000

Example 2:
4.2426  0.0000  0.0000  0.0000
5.1854  6.5659  0.0000  0.0000
12.7279  3.0460  1.6497  0.0000
9.8995  1.6246  1.8497  1.3926


## Smalltalk

 FloatMatrix>>#cholesky	| l |	l := FloatMatrix zero: numRows.	1 to: numRows do: [:i | 		1 to: i do: [:k | | rowSum lkk factor aki partialSum |			i = k				ifTrue: [					rowSum := (1 to: k - 1) sum: [:j | | lkj |						lkj := l at: j @ k.						lkj squared].					lkk := (self at: k @ k) - rowSum.					lkk := lkk sqrt.					l at: k @ k put: lkk]				ifFalse: [					factor := l at: k @ k.					aki := self at: k @ i.					partialSum := (1 to: k - 1) sum: [:j | | ljk lji |						lji := l at: j @ i.						ljk := l at: j @ k.						lji * ljk].					l at: k @ i put: aki - partialSum * factor reciprocal]]].	^l

## Stata

See Cholesky square-root decomposition in Stata help.

mata: a=25,15,-5\15,18,0\-5,0,11 : a[symmetric]        1    2    3    +----------------+  1 |  25            |  2 |  15   18       |  3 |  -5    0   11  |    +----------------+ : cholesky(a)        1    2    3    +----------------+  1 |   5    0    0  |  2 |   3    3    0  |  3 |  -1    1    3  |    +----------------+ : a=18,22,54,42\22,70,86,62\54,86,174,134\42,62,134,106 : a[symmetric]         1     2     3     4    +-------------------------+  1 |   18                    |  2 |   22    70              |  3 |   54    86   174        |  4 |   42    62   134   106  |    +-------------------------+ : cholesky(a)                 1             2             3             4    +---------------------------------------------------------+  1 |  4.242640687             0             0             0  |  2 |  5.185449729   6.565905201             0             0  |  3 |  12.72792206   3.046038495   1.649742248             0  |  4 |  9.899494937   1.624553864   1.849711005   1.392621248  |    +---------------------------------------------------------+

## Swift

Translation of: Rust
func cholesky(matrix: [Double], n: Int) -> [Double] {  var res = [Double](repeating: 0, count: matrix.count)   for i in 0..<n {    for j in 0..<i+1 {      var s = 0.0       for k in 0..<j {        s += res[i * n + k] * res[j * n + k]      }       if i == j {        res[i * n + j] = (matrix[i * n + i] - s).squareRoot()      } else {        res[i * n + j] = (1.0 / res[j * n + j] * (matrix[i * n + j] - s))      }    }  }   return res} func printMatrix(_ matrix: [Double], n: Int) {  for i in 0..<n {    for j in 0..<n {      print(matrix[i * n + j], terminator: " ")    }     print()  }} let res1 = cholesky(  matrix: [25.0, 15.0, -5.0,           15.0, 18.0,  0.0,           -5.0,  0.0, 11.0],  n: 3) let res2 = cholesky(  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],  n: 4) printMatrix(res1, n: 3)print()printMatrix(res2, n: 4)
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

## Tcl

Translation of: Java
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} Demonstration code: 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]
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}


## VBA

This function returns the lower Cholesky decomposition of a square matrix fed to it. It does not check for positive semi-definiteness, although it does check for squareness. It assumes that Option Base 0 is set, and thus the matrix entry indices need to be adjusted if Base is set to 1. It also assumes a matrix of size less than 256x256. To handle larger matrices, change all Byte-type variables to Long. It takes the square matrix range as an input, and can be implemented as an array function on the same sized square range of cells as output. For example, if the matrix is in cells A1:E5, highlighting cells A10:E14, typing "=Cholesky(A1:E5)" and htting Ctrl-Shift-Enter will populate the target cells with the lower Cholesky decomposition.

Function Cholesky(Mat As Range) As Variant Dim A() As Double, L() As Double, sum As Double, sum2 As DoubleDim m As Byte, i As Byte, j As Byte, k As Byte 'Ensure matrix is square    If Mat.Rows.Count <> Mat.Columns.Count Then        MsgBox ("Correlation matrix is not square")        Exit Function    End If     m = Mat.Rows.Count 'Initialize and populate matrix A of values and matrix L which will be the lower Cholesky    ReDim A(0 To m - 1, 0 To m - 1)    ReDim L(0 To m - 1, 0 To m - 1)    For i = 0 To m - 1        For j = 0 To m - 1            A(i, j) = Mat(i + 1, j + 1).Value2            L(i, j) = 0        Next j    Next i 'Handle the simple cases explicitly to save time    Select Case m        Case Is = 1            L(0, 0) = Sqr(A(0, 0))         Case Is = 2            L(0, 0) = Sqr(A(0, 0))            L(1, 0) = A(1, 0) / L(0, 0)            L(1, 1) = Sqr(A(1, 1) - L(1, 0) * L(1, 0))         Case Else            L(0, 0) = Sqr(A(0, 0))            L(1, 0) = A(1, 0) / L(0, 0)            L(1, 1) = Sqr(A(1, 1) - L(1, 0) * L(1, 0))            For i = 2 To m - 1                sum2 = 0                For k = 0 To i - 1                    sum = 0                    For j = 0 To k                        sum = sum + L(i, j) * L(k, j)                    Next j                    L(i, k) = (A(i, k) - sum) / L(k, k)                    sum2 = sum2 + L(i, k) * L(i, k)                Next k                L(i, i) = Sqr(A(i, i) - sum2)            Next i    End Select    Cholesky = LEnd Function

## Wren

Library: Wren-matrix
Library: Wren-fmt
import "/matrix" for Matriximport "/fmt" for Fmt var arrays = [    [ [25, 15, -5],      [15, 18,  0],      [-5,  0, 11] ],     [ [18, 22,  54,  42],      [22, 70,  86,  62],      [54, 86, 174, 134],      [42, 62, 134, 106] ]] for (array in arrays) {    System.print("Original:")    Fmt.mprint(array, 3, 0)    System.print("\nLower Cholesky factor:")    Fmt.mprint(Matrix.new(array).cholesky(), 8, 5)    System.print()}
Output:
Original:
| 25  15  -5|
| 15  18   0|
| -5   0  11|

Lower Cholesky factor:
| 5.00000  0.00000  0.00000|
| 3.00000  3.00000  0.00000|
|-1.00000  1.00000  3.00000|

Original:
| 18  22  54  42|
| 22  70  86  62|
| 54  86 174 134|
| 42  62 134 106|

Lower Cholesky factor:
| 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|


## zkl

Using the GNU Scientific Library:

var [const] GSL=Import("zklGSL");	// libGSL (GNU Scientific Library)fcn lowerCholesky(m){  // trans: C   rows:=m.rows;   lcm:=GSL.Matrix(rows,rows);	// zero filled   foreach i,j in (rows,i+1){       s:=(0).reduce(j,'wrap(s,k){ s + lcm[i,k]*lcm[j,k] },0.0);      lcm[i,j]=( if(i==j)(m[i,i] - s).sqrt()	         else     1.0/lcm[j,j]*(m[i,j] - s) );   }   lcm}
Output:
lowerCholesky(GSL.Matrix(3,3).set(25, 15, -5, 	// example 1
15, 18,  0,
-5,  0, 11))
.format(6).println();
5.00,  0.00,  0.00
3.00,  3.00,  0.00
-1.00,  1.00,  3.00

Output:
lowerCholesky(GSL.Matrix(4,4).set(	// example 2
18, 22,  54,  42,
22, 70,  86,  62,
54, 86, 174, 134,
42, 62, 134, 106) )
.format(8,4).println();
4.2426,  0.0000,  0.0000,  0.0000
5.1854,  6.5659,  0.0000,  0.0000
12.7279,  3.0460,  1.6497,  0.0000
9.8995,  1.6246,  1.8497,  1.3926


Or, using lists:

Translation of: C
fcn cholesky(mat){   rows:=mat.len();   r:=(0).pump(rows,List().write, (0).pump(rows,List,0.0).copy); // matrix of zeros   foreach i,j in (rows,i+1){       s:=(0).reduce(j,'wrap(s,k){ s + r[i][k]*r[j][k] },0.0);      r[i][j]=( if(i==j)(mat[i][i] - s).sqrt()	        else    1.0/r[j][j]*(mat[i][j] - s) );   }   r}
ex1:=L( L(25.0,15.0,-5.0), L(15.0,18.0,0.0), L(-5.0,0.0,11.0) );printM(cholesky(ex1));println("-----------------");ex2:=L( L(18.0, 22.0,  54.0,  42.0,),         L(22.0, 70.0,  86.0,  62.0,), 	L(54.0, 86.0, 174.0, 134.0,), 	L(42.0, 62.0, 134.0, 106.0,) );printM(cholesky(ex2));
fcn printM(m){ m.pump(Console.println,rowFmt) }fcn rowFmt(row){ ("%9.5f "*row.len()).fmt(row.xplode()) }
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


## ZX Spectrum Basic

Translation of: BBC_BASIC
10 LET d=2000: GO SUB 1000: GO SUB 4000: GO SUB 500020 LET d=3000: GO SUB 1000: GO SUB 4000: GO SUB 500030 STOP 1000 RESTORE d1010 READ a,b1020 DIM m(a,b)1040 FOR i=1 TO a1050 FOR j=1 TO b1060 READ m(i,j)1070 NEXT j1080 NEXT i1090 RETURN 2000 DATA 3,3,25,15,-5,15,18,0,-5,0,113000 DATA 4,4,18,22,54,42,22,70,86,62,54,86,174,134,42,62,134,1064000 REM Cholesky decomposition4005 DIM l(a,b)4010 FOR i=1 TO a4020 FOR j=1 TO i4030 LET s=04050 FOR k=1 TO j-14060 LET s=s+l(i,k)*l(j,k)4070 NEXT k4080 IF i=j THEN LET l(i,j)=SQR (m(i,i)-s): GO TO 41004090 LET l(i,j)=(m(i,j)-s)/l(j,j)4100 NEXT j4110 NEXT i4120 RETURN 5000 REM Print5010 FOR r=1 TO a5020 FOR c=1 TO b5030 PRINT l(r,c);" ";5040 NEXT c5050 PRINT 5060 NEXT r5070 RETURN