Cholesky decomposition: Difference between revisions

Content added Content deleted

Inline

Latest revision as of 12:57, 17 November 2023

Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:

A=LL^{T}

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

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

{\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 ( $l_{kk}$ ) of $L$ there is a calculation pattern:

l_{11}={\sqrt {a_{11}}}

l_{22}={\sqrt {a_{22}-l_{21}^{2}}}

l_{33}={\sqrt {a_{33}-(l_{31}^{2}+l_{32}^{2})}}

or in general:

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

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

l_{21}={\frac {1}{l_{11}}}a_{21}

l_{31}={\frac {1}{l_{11}}}a_{31}

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

which can also be expressed in a general formula:

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

Task description

The task is to implement a routine which will return a lower Cholesky factor $L$ for every given symmetric, positive definite nxn matrix $A$ . You should then test it on the following two examples and include your output.

Example 1:

25  15  -5                 5   0   0
15  18   0         -->     3   3   0
-5   0  11                -1   1   3

Example 2:

18  22   54   42           4.24264    0.00000    0.00000    0.00000
22  70   86   62   -->     5.18545    6.56591    0.00000    0.00000
54  86  174  134          12.72792    3.04604    1.64974    0.00000
42  62  134  106           9.89949    1.62455    1.84971    1.39262

Note

The Cholesky decomposition of a Pascal upper-triangle matrix is the Identity matrix of the same size.
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]
]

Ada

Works with: Ada 2005

decomposition.ads:

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;

decomposition.adb:

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

Arturo

cholesky: function [m][
    result: array.of: @[size m, size m] 0.0

    loop 0..dec size m\0 'i [
        loop 0..i 'j [
            s: 0.0
            loop 0..j 'k ->
                s: s + result\[i]\[k] * result\[j]\[k]

            result\[i]\[j]: (i = j)? -> sqrt m\[i]\[i] - s
                                     -> (1.0 // result\[j]\[j]) * (m\[i]\[j] - s)
        ]
    ]
    return result
]

printMatrix: function [a]->
    loop a 'b ->
        print to [:string] .format:"8.5f" b

m1: @[
    @[25.0, 15.0, neg 5.0]
    @[15.0, 18.0,  0.0]
    @[neg 5.0,  0.0, 11.0]
]
printMatrix cholesky m1

print ""

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

ATS

%{^
#include <math.h>
#include <float.h>
%}

#include "share/atspre_staload.hats"

macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1

(* The sqrt(3) function made part of the ‘g0float’ typekind series.
   (The ats2-xprelude package will do this for you, but it is easy
   to do if you are not using a lot of math functions. *)
extern fn {tk : tkind} g0float_sqrt : g0float tk -<> g0float tk
overload sqrt with g0float_sqrt
implement g0float_sqrt<fltknd> x = $extfcall (float, "sqrtf", x)
implement g0float_sqrt<dblknd> x = $extfcall (double, "sqrt", x)
implement g0float_sqrt<ldblknd> x = $extfcall (ldouble, "sqrtl", x)

(*------------------------------------------------------------------*)
(* A "very little matrix library"                                   *)

typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (int i1, int j1) -<cloref0>
  [i0, j0 : pos | i0 <= m0; j0 <= n0]
    @(int i0, int j0)

datatype Real_Matrix (tk : tkind,
                      m1 : int, n1 : int,
                      m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
                  int m1, int n1, int m0, int n0,
                  Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
  [m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
  [m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)

extern fn {tk : tkind}
Real_Matrix_make_elt :
  {m0, n0 : pos}
  (int m0, int n0, g0float tk) -< !wrt >
    Real_Matrix (tk, m0, n0, m0, n0)

extern fn {tk : tkind}
Real_Matrix_copy :
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)

extern fn {tk : tkind}
Real_Matrix_copy_to :
  {m1, n1 : pos}
  (Real_Matrix (tk, m1, n1),    (* destination *)
   Real_Matrix (tk, m1, n1)) -< !refwrt >
    void

extern fn {}
Real_Matrix_dimension :
  {tk : tkind}
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)

extern fn {tk : tkind}
Real_Matrix_get_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk

extern fn {tk : tkind}
Real_Matrix_set_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
    void

extern fn {}
Real_Matrix_reflect_lower_triangle :
  (* This operation makes every It is a change in how INDEXING
     works. All the storage is still in the lower triangle. *)
  {tk     : tkind}
  {n1     : pos}
  {m0, n0 : pos}
  Real_Matrix (tk, n1, n1, m0, n0) -<>
    Real_Matrix (tk, n1, n1, m0, n0)  

extern fn {tk : tkind}
Real_Matrix_fprint :
  {m, n : pos}
  (FILEref, Real_Matrix (tk, m, n)) -<1> void

overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload reflect_lower_triangle with
  Real_Matrix_reflect_lower_triangle

(*------------------------------------------------------------------*)
(* Implementation of the "very little matrix library"               *)

implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
  Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
               m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))

implement {}
Real_Matrix_dimension A =
  case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)

implement {tk}
Real_Matrix_get_at (A, i1, j1) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
  end

implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
  end

implement {}
Real_Matrix_reflect_lower_triangle {..} {n1} A =
  let
    typedef t = intBtwe (1, n1)
    val+ Real_Matrix (storage, n1, _, m0, n0, index_map) = A
  in
    Real_Matrix (storage, n1, n1, m0, n0,
                 lam (i, j) =>
                  index_map ((if j <= i then i else j) : t,
                             (if j <= i then j else i) : t))
  end

implement {tk}
Real_Matrix_copy A =
  let
    val @(m1, n1) = dimension A
    val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
    val () = copy_to<tk> (C, A)
  in
    C
  end

implement {tk}
Real_Matrix_copy_to (Dst, Src) =
  let
    val @(m1, n1) = dimension Src
    prval [m1 : int] EQINT () = eqint_make_gint m1
    prval [n1 : int] EQINT () = eqint_make_gint n1

    var i : intGte 1
  in
    for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m1; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n1; j := succ j)
              Dst[i, j] := Src[i, j]
        end
  end

implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
  let
    val @(m, n) = dimension A
    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              let
                typedef FILEstar = $extype"FILE *"
                extern castfn FILEref2star : FILEref -<> FILEstar
                val _ = $extfcall (int, "fprintf", FILEref2star outf,
                                   "%16.6g", A[i, j])
              in
              end;
          fprintln! (outf)
        end
  end

(*------------------------------------------------------------------*)
(* Cholesky-Banachiewicz, in place. See
   https://en.wikipedia.org/w/index.php?title=Cholesky_decomposition&oldid=1149960985#The_Cholesky%E2%80%93Banachiewicz_and_Cholesky%E2%80%93Crout_algorithms

   I would use Cholesky-Crout if my matrices were stored in column
   major order. But it makes little difference. *)

extern fn {tk : tkind}
Real_Matrix_cholesky_decomposition :
  (* Only the lower triangle is considered. *)
  {n : pos}
  Real_Matrix (tk, n, n) -< !refwrt > void

overload cholesky_decomposition with
  Real_Matrix_cholesky_decomposition

implement {tk}
Real_Matrix_cholesky_decomposition {n} A =
  let
    val @(n, _) = dimension A

    (* I arrange the nested loops somewhat differently from how it is
       done in the Wikipedia article's C snippet. *)
    fun
    repeat {i, j : pos | j <= i; i <= n + 1} (* <-- allowed values *)
           .<(n + 1) - i, i - j>. (* <-- proof of termination *)
           (i : int i, j : int j) :<!refwrt> void =
      if i = n + 1 then
        ()                      (* All done. *)
      else
        let
          fun
          _sum {k : pos | k <= j} .<j - k>.
               (x : g0float tk, k : int k) :<!refwrt> g0float tk =
            if k = j then
              x
            else
              _sum (x + (A[i, k] * A[j, k]), succ k)

          val sum = _sum (Zero, 1)
        in
          if j = i then
            begin
              A[i, j] := sqrt (A[i, i] - sum);
              repeat (succ i, 1)
            end
          else
            begin
              A[i, j] := (One / A[j, j]) * (A[i, j] - sum);
              repeat (i, succ j)
            end
        end
  in
    repeat (1, 1)
  end

(*------------------------------------------------------------------*)

fn {tk : tkind}           (* We like Fortran, so COLUMN major here. *)
column_major_list_to_square_matrix
          {n   : pos}
          (n   : int n,
           lst : list (g0float tk, n * n))
    : Real_Matrix (tk, n, n) =
  let
    #define :: list_cons
    prval () = mul_gte_gte_gte {n, n} ()
    val A = Real_Matrix_make_elt (n, n, NAN)
    val lstref : ref (List0 (g0float tk)) = ref lst
    var j : intGte 1
  in
    for* {j : pos | j <= n + 1} .<(n + 1) - j>.
         (j : int j) =>
      (j := 1; j <> succ n; j := succ j)
        let
          var i : intGte 1
        in
          for* {i : pos | i <= n + 1} .<(n + 1) - i>.
               (i : int i) =>
            (i := 1; i <> succ n; i := succ i)
              case- !lstref of
              | hd :: tl =>
                begin
                  A[i, j] := hd;
                  !lstref := tl
                end
        end;
    A
  end

implement
main0 () =
  let
    val _A =
      column_major_list_to_square_matrix
        (3, $list (25.0, 15.0, ~5.0,
                   0.0, 18.0, 0.0,
                   0.0, 0.0, 11.0))
    val A = reflect_lower_triangle _A
    and B = copy _A
    val () =
      begin
        cholesky_decomposition B;
        print! ("\nThe Cholesky decomposition of\n\n");
        Real_Matrix_fprint (stdout_ref, A);
        print! ("is\n");
        Real_Matrix_fprint (stdout_ref, B)
      end

    val _A =
      column_major_list_to_square_matrix
        (4, $list (18.0, 22.0, 54.0, 42.0,
                   0.0, 70.0, 86.0, 62.0,
                   0.0, 0.0, 174.0, 134.0,
                   0.0, 0.0, 0.0, 106.0))
    val A = reflect_lower_triangle _A
    and B = copy _A
    val () =
      begin
        cholesky_decomposition B;
        print! ("\nThe Cholesky decomposition of\n\n");
        Real_Matrix_fprint (stdout_ref, A);
        print! ("is\n");
        Real_Matrix_fprint (stdout_ref, B)
      end
  in
    println! ()
  end

(*------------------------------------------------------------------*)

Output:

$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_GCBDW cholesky_decomposition_task.dats -lgc -lm && ./a.out

The Cholesky decomposition of

              25              15              -5
              15              18               0
              -5               0              11
is
               5               0               0
               3               3               0
              -1               1               3

The Cholesky decomposition of

              18              22              54              42
              22              70              86              62
              54              86             174             134
              42              62             134             106
is
         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

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

Works with: BBC BASIC for Windows

      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 1
CL-USER> (setf a '((25 15 -5) (15 18 0) (-5 0 11)))
((25 15 -5) (15 18 0) (-5 0 11))
CL-USER> (cholesky a)
((5 0 0) (3 3 0) (-1 1 3))
CL-USER> (format t "~{~{~5d~}~%~}" (cholesky a))
    5    0    0
    3    3    0
   -1    1    3
NIL

;; example 2
CL-USER> (setf a '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106)))
((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))
CL-USER> (cholesky a)
((4.2426405 0 0 0) (5.18545 6.565905 0 0) (12.727922 3.0460374 1.6497375 0) (9.899495 1.6245536 1.849715 1.3926151))
CL-USER> (format t "~{~{~10,5f~}~%~}" (cholesky a))
   4.24264   0.00000   0.00000   0.00000
   5.18545   6.56591   0.00000   0.00000
  12.72792   3.04604   1.64974   0.00000
   9.89950   1.62455   1.84971   1.39262
NIL

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 !rows
INTEGER, PARAMETER :: n=3 !cols
COMPLEX, DIMENSION(m,n) :: A 
REAL, DIMENSION(m,n) :: L
REAL :: sum1, sum2
INTEGER i,j,k

! Assign values to the matrix
A(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 values
L(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 do
end do

! write output
do 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())

Print
Cholesky_decomp(m2())
Print_(m2())

' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
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

Frink

Frink's package Matrix.frink contains routines for Cholesky-Crout decomposition of a square Hermitian matrix (which can be real or complex.) This code is adapted from that more powerful class to work on raw 2-dimensional arrays. This also demonstrates Frink's layout routines.

Cholesky[array] :=
{
   n = length[array]
   L = new array[[n,n], 0]

   for j = 0 to n-1
   {
      sum = 0
      for k = 0 to j-1
         sum = sum + (L@j@k)^2

      L@j@j = sqrt[array@j@j - sum]

      for i = j+1 to n-1
      {
         sum = 0
         for k = 0 to j-1
            sum = sum + L@i@k * L@j@k

         L@i@j = (1 / L@j@j * (array@i@j -sum))
      }
   }

   return L
}

A = [[  25, 15, -5],
     [  15, 18,  0],
     [  -5,  0, 11]]

println[formatTable[[[formatMatrix[A], "->", formatMatrix[Cholesky[A]]]]]]

B = [[18,  22,  54,  42],
     [22,  70,  86,  62],
     [54,  86, 174, 134],
     [42,  62, 134, 106]]

println[formatTable[[[formatMatrix[B], "->", formatMatrix[formatFix[Cholesky[B], 1, 5]]]]]]

Output:

┌          ┐    ┌        ┐
│25  15  -5│    │ 5  0  0│
│          │    │        │
│15  18   0│ -> │ 3  3  0│
│          │    │        │
│-5   0  11│    │-1  1  3│
└          ┘    └        ┘
┌                ┐    ┌                                   ┐
│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│
└                ┘    └                                   ┘

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.IArray
import Data.Array.MArray
import Data.Array.Unboxed
import Data.Array.ST

type Idx = (Int,Int)
type Arr = UArray Idx Double

-- Return the (i,j) element of the lower triangular matrix.  (We assume the
-- lower array bound is (0,0).)
get :: Arr -> Arr -> Idx -> Double
get a l (i,j) | i == j = sqrt $ a!(j,j) - dot
              | i  > j = (a!(i,j) - dot) / l!(j,j)
              | otherwise = 0
  where dot = sum [l!(i,k) * l!(j,k) | k <- [0..j-1]]

-- Return the lower triangular matrix of a Cholesky decomposition.  We assume
-- the input is a real, symmetric, positive-definite matrix, with lower array
-- bounds of (0,0).
cholesky :: Arr -> Arr
cholesky a = let n = maxBnd a
             in runSTUArray $ do
               l <- thaw a
               mapM_ (update a l) [(i,j) | i <- [0..n], j <- [0..n]]
               return l
  where maxBnd = fst . snd . bounds
        update a l i = unsafeFreeze l >>= \l' -> writeArray l i (get a l' i)

The main module:

import Data.Array.IArray
import Data.List
import Cholesky

fm _ [] = "" 
fm _ [x] = fst x
fm width ((a,b):xs) = a ++ (take (width - b) $ cycle " ") ++ (fm width xs)

fmt width row (xs,[]) = fm width xs
fmt 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 :: Arr
ex1 = listArray ((0,0),(2,2)) [25, 15, -5, 
                               15, 18,  0, 
                               -5,  0, 11]

ex2 = listArray ((0,0),(3,3)) [18, 22,  54,  42, 
                               22, 70,  86,  62, 
                               54, 86, 174, 134, 
                               42, 62, 134, 106]

main :: IO ()
main = do
  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 R
a = (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 result
end

procedure make_square_array (n)
  result := []
  every (1 to n) do push (result, list(n, 0))
  return result
end

procedure print_array (array)
  every (row := !array) do {
    every writes (!row || " ")
    write ()
  }
end

procedure do_cholesky (array)
  write ("Input:")
  print_array (array)
  result := cholesky (array)
  write ("Result:")
  print_array (result)
end

procedure main ()
  do_cholesky ([[25,15,-5],[15,18,0],[-5,0,11]])
  do_cholesky ([[18,22,54,42],[22,70,86,62],[54,86,174,134],[42,62,134,106]])
end

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 -> Type
Matrix m n t = Vect m (Vect n t)


zeros : (m : Nat) -> (n : Nat) -> Matrix m n Double
zeros m n = replicate m (replicate n 0.0)


indexM : (Fin m, Fin n) -> Matrix m n t -> t
indexM (i, j) a = index j (index i a)


replaceAtM : (Fin m, Fin n) -> t -> Matrix m n t -> Matrix m n t
replaceAtM (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) -> Double
get 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 Double
updateL a l idx = replaceAtM idx (get a l idx) l


cholesky : Matrix m m Double -> Matrix m m Double
cholesky 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 Double
ex1 = cholesky [[25.0, 15.0, -5.0], [15.0, 18.0, 0.0], [-5.0, 0.0, 11.0]] 

ex2 : Matrix 4 4 Double
ex2 = 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 product
h =: +@|:    NB. conjugate transpose

cholesky=: 3 : 0
 n=. #A=. y
 if. 1>:n do.
  assert. (A=|A)>0=A  NB. check for positive definite
  %:A
 else.
  '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 eg2
4.24264       0       0       0
5.18545 6.56591       0       0
12.7279 3.04604 1.64974       0
9.89949 1.62455 1.84971 1.39262

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_ eg2
4.24264       0       0       0
5.18545 6.56591       0       0
12.7279 3.04604 1.64974       0
9.89949 1.62455 1.84971 1.39262

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 truncation
def 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 ;

Task 1:

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

Output:

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

Task 2:

[[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 checked
def 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 lower
def 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

Maxima

/* 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_vec
let test a =
   print_endline "\nin:"; show a;
   print_endline "out:"; show (cholesky a)

let _ =
   test [| [|25.0; 15.0; -5.0|];
           [|15.0; 18.0;  0.0|];
           [|-5.0;  0.0; 11.0|] |];
   test [| [|18.0; 22.0;  54.0;  42.0|];
           [|22.0; 70.0;  86.0;  62.0|];
           [|54.0; 86.0; 174.0; 134.0|];
           [|42.0; 62.0; 134.0; 106.0|] |];

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 niner
hexer =  18  22  54  42,                             /*define a  4x4  matrix. */
         22  70  86  62,
         54  86 174 134,
         42  62 134 106
                           call Cholesky hexer
exit                                   /*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 pprint
from 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)

sub cholesky(@A) {
    my @L = @A »×» 0;
    for ^@A -> \i {
        for 0..i -> \j {
            @L[i;j] = (i == j ?? &sqrt !! 1/@L[j;j] × * )\      # select function
                      (@A[i;j] - [+] (@L[i;*] Z× @L[j;*])[^j])  # provide value
        }
    }
    @L
}

.fmt('%3d').say for cholesky [
    [25],
    [15, 18],
    [-5,  0, 11],
];

say '';

.fmt('%6.3f').say for cholesky [
    [18, 22,  54,  42],       
    [22, 70,  86,  62],
    [54, 86, 174, 134],       
    [42, 62, 134, 106],
];

Output:

  5
  3   3
 -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

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 niner
hexer =  18  22  54  42,                         /*define a  4x4  matrix with elements. */
         22  70  86  62,
         54  86 174 134,
         42  62 134 106
                           call Cholesky hexer
exit                                             /*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 cholesky_factor
    raise ArgumentError, "must provide symmetric matrix" unless symmetric?
    l = Array.new(row_size) {Array.new(row_size, 0)}
    (0 ... row_size).each do |k|
      (0 ... row_size).each do |i|
        if i == k
          sum = (0 .. k-1).inject(0.0) {|sum, j| sum + l[k][j] ** 2}
          val = Math.sqrt(self[k,k] - sum)
          l[k][k] = val
        elsif i > k
          sum = (0 .. k-1).inject(0.0) {|sum, j| sum + l[i][j] * l[k][j]}
          val = (self[k,i] - sum) / l[k][k]
          l[i][k] = val
        end
      end
    end
    Matrix[*l]
  end
end

puts Matrix[[25,15,-5],[15,18,0],[-5,0,11]].cholesky_factor
puts Matrix[[18, 22,  54,  42],
            [22, 70,  86,  62],
            [54, 86, 174, 134],
            [42, 62, 134, 106]].cholesky_factor

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.cholesky
val l2 = a2.cholesky


// Given test results
val 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 Double
Dim 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 = L
End Function

V (Vlang)

Translation of: go

import math

// Symmetric and Lower use a packed representation that stores only
// the Lower triangle.
 
struct Symmetric {
    order int
    ele   []f64
}
 
struct Lower  {
mut:
    order int
    ele   []f64
}
 
// Symmetric.print prints a square matrix from the packed representation,
// printing the upper triange as a transpose of the Lower.
fn (s Symmetric) print() {
    mut row, mut diag := 1, 0
    for i, e in s.ele {
        print("${e:10.5f} ")
        if i == diag {
            for j, col := diag+row, row; col < s.order; j += col {
                print("${s.ele[j]:10.5f} ")
                col++
            }
            println('')
            row++
            diag += row
        }
    }
}
 
// Lower.print prints a square matrix from the packed representation,
// printing the upper triangle as all zeros.
fn (l Lower) print() {
    mut row, mut diag := 1, 0
    for i, e in l.ele {
        print("${e:10.5f} ")
        if i == diag {
            for _ in row..l.order {
                print("${0.0:10.5f} ")
            }
            println('')
            row++
            diag += row
        }
    }
}
 
// cholesky_lower returns the cholesky decomposition of a Symmetric real
// matrix.  The matrix must be positive definite but this is not checked.
fn (a Symmetric) cholesky_lower() Lower {
    mut l := Lower{a.order, []f64{len: a.ele.len}}
    mut row, mut col := 1, 1
    mut dr := 0 // index of diagonal element at end of row
    mut dc := 0 // index of diagonal element at top of column
    for i, e in a.ele {
        if i < dr {
            d := (e - l.ele[i]) / l.ele[dc]
            l.ele[i] = d
            mut ci, mut 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
}
 
fn main() {
    demo(Symmetric{3, [
        f64(25),
        15, 18,
        -5, 0, 11]})
    demo(Symmetric{4, [
        f64(18),
        22, 70,
        54, 86, 174,
        42, 62, 134, 106]})
}
 
fn demo(a Symmetric) {
    println("A:")
    a.print()
    println("L:")
    a.cholesky_lower().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.3926

Wren

Library: Wren-matrix

Library: Wren-fmt

import "./matrix" for Matrix
import "./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|

XPL0

real L(4*4);

func real Cholesky(A, N);
real A;  int N;
real S;
int  I, J, K;
[for I:= 0 to N*N-1 do L(I):= 0.;
for I:= 0 to N-1 do
    for J:= 0 to I do
        [S:= 0.;
        for K:= 0 to J-1 do
            S:= S + L(I*N+K) * L(J*N+K);
        L(I*N+J):= if I = J then sqrt(A(I*N+I) - S)
            else (1.0 / L(J*N+J) * (A(I*N+J) - S));
        ];
return L;
];

proc ShowMatrix(A, N);
real A;  int N;
int  I, J;
[for I:= 0 to N-1 do
    [for J:= 0 to N-1 do
        RlOut(0, A(I*N+J));
    CrLf(0);
    ];
];

int  N;
real M1, C1, M2, C2;
[N:= 3;
M1:=   [25., 15., -5.,
        15., 18.,  0.,
        -5.,  0., 11.];
C1:= Cholesky(M1, N);
ShowMatrix(C1, N);
CrLf(0);

N:= 4;
M2:=   [18., 22.,  54.,  42.,
        22., 70.,  86.,  62.,
        54., 86., 174., 134.,
        42., 62., 134., 106.];
C2:= Cholesky(M2, N);
ShowMatrix(C2, N);
]

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

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 5000
20 LET d=3000: GO SUB 1000: GO SUB 4000: GO SUB 5000
30 STOP 
1000 RESTORE d
1010 READ a,b
1020 DIM m(a,b)
1040 FOR i=1 TO a
1050 FOR j=1 TO b
1060 READ m(i,j)
1070 NEXT j
1080 NEXT i
1090 RETURN 
2000 DATA 3,3,25,15,-5,15,18,0,-5,0,11
3000 DATA 4,4,18,22,54,42,22,70,86,62,54,86,174,134,42,62,134,106
4000 REM Cholesky decomposition
4005 DIM l(a,b)
4010 FOR i=1 TO a
4020 FOR j=1 TO i
4030 LET s=0
4050 FOR k=1 TO j-1
4060 LET s=s+l(i,k)*l(j,k)
4070 NEXT k
4080 IF i=j THEN LET l(i,j)=SQR (m(i,i)-s): GO TO 4100
4090 LET l(i,j)=(m(i,j)-s)/l(j,j)
4100 NEXT j
4110 NEXT i
4120 RETURN 
5000 REM Print
5010 FOR r=1 TO a
5020 FOR c=1 TO b
5030 PRINT l(r,c);" ";
5040 NEXT c
5050 PRINT 
5060 NEXT r
5070 RETURN