Cholesky decomposition: Difference between revisions
m
→{{header|Wren}}: Minor tidy
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: disambiguate with '×', clarify function/value setup, add output) |
m (→{{header|Wren}}: Minor tidy) |
||
| (9 intermediate revisions by 8 users not shown) | |||
Line 78:
# The Cholesky decomposition of a [[Pascal matrix generation|Pascal]] upper-triangle matrix is the [[wp:Identity matrix|Identity matrix]] of the same size.
# The Cholesky decomposition of a Pascal symmetric matrix is the Pascal lower-triangle matrix of the same size.
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F cholesky(A)
V l = [[0.0] * A.len] * A.len
L(i) 0 .< A.len
Line 119 ⟶ 118:
]
</pre>
=={{header|Ada}}==
{{works with|Ada 2005}}
decomposition.ads:
<syntaxhighlight lang=
generic
with package Matrix is new Ada.Numerics.Generic_Real_Arrays (<>);
Line 134 ⟶ 132:
decomposition.adb:
<syntaxhighlight lang=
package body Decomposition is
Line 169 ⟶ 167:
Example usage:
<syntaxhighlight lang=
with Ada.Text_IO;
with Decomposition;
Line 236 ⟶ 234:
12.728 3.046 1.650 0.000
9.899 1.625 1.850 1.393</pre>
=={{header|ALGOL 68}}==
{{trans|C}} Note: This specimen retains the original [[#C|C]] coding style. [http://rosettacode.org/mw/index.php?title=Cholesky_decomposition&action=historysubmit&diff=107753&oldid=107752 diff]
Line 242 ⟶ 239:
{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny].}}
{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}}
<syntaxhighlight lang="algol68">#!/usr/local/bin/a68g --script #
MODE FIELD=LONG REAL;
Line 310 ⟶ 307:
(12.72792, 3.04604, 1.64974, 0.00000),
( 9.89949, 1.62455, 1.84971, 1.39262))
</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="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</syntaxhighlight>
{{out}}
<pre> 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</pre>
=={{header|ATS}}==
<syntaxhighlight lang="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
(*------------------------------------------------------------------*)
</syntaxhighlight>
{{out}}
<pre>$ 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
</pre>
=={{header|AutoHotkey}}==
<syntaxhighlight lang=
L := [], n := A.Count()
L[1,1] := Sqrt(A[1,1])
Line 348 ⟶ 757:
return "[" Trim(output, "`n,") "]"
}</syntaxhighlight>
Examples:<syntaxhighlight lang=
, [15, 18, 0]
, [-5, 0 , 11]]
Line 374 ⟶ 783:
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
<syntaxhighlight lang="bbcbasic"> DIM m1(2,2)
m1() = 25, 15, -5, \
\ 15, 18, 0, \
Line 431 ⟶ 840:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|C}}==
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <math.h>
Line 493 ⟶ 901:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|C sharp|C#}}==
<syntaxhighlight lang="csharp">using System;
using System.Collections.Generic;
using System.Linq;
Line 624 ⟶ 1,031:
12.72792, 3.04604, 1.64974, 0.00000,
9.89949, 1.62455, 1.84971, 1.39262,
=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <cassert>
#include <cmath>
#include <iomanip>
Line 750 ⟶ 1,156:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Clojure}}==
{{trans|Python}}
<syntaxhighlight lang="clojure">(defn cholesky
[matrix]
(let [n (count matrix)
Line 765 ⟶ 1,170:
(vec (map vec L))))</syntaxhighlight>
Example:
<syntaxhighlight lang="clojure">(cholesky [[25 15 -5] [15 18 0] [-5 0 11]])
;=> [[ 5.0 0.0 0.0]
; [ 3.0 3.0 0.0]
Line 775 ⟶ 1,180:
; [12.727922061357857 3.0460384954008553 1.6497422479090704 0.0 ]
; [ 9.899494936611667 1.624553864213788 1.8497110052313648 1.3926212476456026]]</syntaxhighlight>
=={{header|Common Lisp}}==
<syntaxhighlight lang="lisp">;; Calculates the Cholesky decomposition matrix L
;; for a positive-definite, symmetric nxn matrix A.
(defun chol (A)
Line 807 ⟶ 1,211:
L))</syntaxhighlight>
<syntaxhighlight lang="lisp">;; Example 1:
(setf A (make-array '(3 3) :initial-contents '((25 15 -5) (15 18 0) (-5 0 11))))
(chol A)
Line 814 ⟶ 1,218:
(-1.0 1.0 3.0))</syntaxhighlight>
<syntaxhighlight lang="lisp">;; Example 2:
(setf B (make-array '(4 4) :initial-contents '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))))
(chol B)
Line 822 ⟶ 1,226:
(9.899495 1.6245536 1.849715 1.3926151))</syntaxhighlight>
<syntaxhighlight lang="lisp">;; case of matrix stored as a list of lists (inner lists are rows of matrix)
;; as above, returns the Cholesky decomposition matrix of a square positive-definite, symmetric matrix
(defun cholesky (m)
Line 834 ⟶ 1,238:
(defun *v (v1 v2) (reduce #'+ (mapcar #'* v1 v2)))</syntaxhighlight>
<syntaxhighlight lang="lisp">;; example 1
CL-USER> (setf a '((25 15 -5) (15 18 0) (-5 0 11)))
((25 15 -5) (15 18 0) (-5 0 11))
Line 845 ⟶ 1,249:
NIL</syntaxhighlight>
<syntaxhighlight lang="lisp">;; example 2
CL-USER> (setf a '((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106)))
((18 22 54 42) (22 70 86 62) (54 86 174 134) (42 62 134 106))
Line 856 ⟶ 1,260:
9.89950 1.62455 1.84971 1.39262
NIL</syntaxhighlight>
=={{header|D}}==
<syntaxhighlight lang="d">import std.stdio, std.math, std.numeric;
T[][] cholesky(T)(in T[][] A) pure nothrow /*@safe*/ {
Line 898 ⟶ 1,301:
=={{header|DWScript}}==
{{Trans|C}}
<syntaxhighlight lang="delphi">function Cholesky(a : array of Float) : array of Float;
var
i, j, k, n : Integer;
Line 944 ⟶ 1,347:
var c2 := Cholesky(m2);
ShowMatrix(c2);</syntaxhighlight>
=={{header|F_Sharp|F#}}==
<syntaxhighlight lang="fsharp">open Microsoft.FSharp.Collections
let cholesky a =
Line 1,004 ⟶ 1,406:
42.00000, 62.00000, 134.00000, 106.00000, 9.89949, 1.62455, 1.84971, 1.39262,
</pre>
=={{header|Fantom}}==
<syntaxhighlight lang="fantom">**
** Cholesky decomposition
**
Line 1,068 ⟶ 1,469:
[[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]]
</pre>
=={{header|Fortran}}==
<syntaxhighlight lang=
! *************************************************!
! LBH @ ULPGC 06/03/2014
Line 1,146 ⟶ 1,546:
-1.0 1.0 3.0
</pre>
=={{header|FreeBASIC}}==
{{trans|BBC BASIC}}
<syntaxhighlight lang="freebasic">' version 18-01-2017
' compile with: fbc -s console
Line 1,224 ⟶ 1,623:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|Frink}}==
Frink's package [https://frinklang.org/fsp/colorize.fsp?f=Matrix.frink 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.
<syntaxhighlight lang="frink">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]]]]]]</syntaxhighlight>
{{out}}
<pre>
┌ ┐ ┌ ┐
│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│
└ ┘ └ ┘
</pre>
=={{header|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.
<syntaxhighlight lang="go">package main
import (
Line 1,356 ⟶ 1,816:
===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.
<syntaxhighlight lang="go">package main
import (
Line 1,465 ⟶ 1,925:
===Library gonum/mat===
<syntaxhighlight lang="go">package main
import (
Line 1,505 ⟶ 1,965:
===Library go.matrix===
<syntaxhighlight lang="go">package main
import (
Line 1,559 ⟶ 2,019:
9.899495, 1.624554, 1.849711, 1.392621}
</pre>
=={{header|Groovy}}==
{{Trans|Java}}
<syntaxhighlight lang="groovy">def decompose = { a ->
assert a.size > 0 && a[0].size == a.size
def m = a.size
Line 1,577 ⟶ 2,036:
}</syntaxhighlight>
Test:
<syntaxhighlight lang="groovy">def test1 = [[25, 15, -5],
[15, 18, 0],
[-5, 0, 11]]
Line 1,599 ⟶ 2,058:
[12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0]
[9.899494936611667, 1.624553864213788, 1.8497110052313648, 1.3926212476456026]</pre>
=={{header|Haskell}}==
We use the [http://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky.E2.80.93Banachiewicz_and_Cholesky.E2.80.93Crout_algorithms Cholesky–Banachiewicz algorithm] described in the Wikipedia article.
Line 1,606 ⟶ 2,064:
The Cholesky module:
<syntaxhighlight lang="haskell">module Cholesky (Arr, cholesky) where
import Data.Array.IArray
Line 1,636 ⟶ 2,094:
update a l i = unsafeFreeze l >>= \l' -> writeArray l i (get a l' i)</syntaxhighlight>
The main module:
<syntaxhighlight lang="haskell">import Data.Array.IArray
import Data.List
import Cholesky
Line 1,678 ⟶ 2,136:
===With Numeric.LinearAlgebra===
<syntaxhighlight lang="haskell">import Numeric.LinearAlgebra
a,b :: Matrix R
Line 1,726 ⟶ 2,184:
, 12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0.0
, 9.899494936611665, 1.6245538642137891, 1.849711005231382, 1.3926212476455904 ]</pre>
=={{header|Icon}} and {{header|Unicon}}==
<syntaxhighlight lang=
result := make_square_array (*array)
every (i := 1 to *array) do {
Line 1,790 ⟶ 2,247:
9.899494937 1.624553864 1.849711005 1.392621248
</pre>
=={{header|Idris}}==
'''works with Idris 0.10'''
'''Solution:'''
<syntaxhighlight lang=
import Data.Vect
Line 1,870 ⟶ 2,326:
[[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]]
</pre>
=={{header|J}}==
'''Solution:'''
<syntaxhighlight lang="j">mp=: +/ . * NB. matrix product
h =: +@|: NB. conjugate transpose
Line 1,890 ⟶ 2,345:
See [[j:Essays/Cholesky Decomposition|Cholesky Decomposition essay]] on the J Wiki.
{{out|Examples}}
<syntaxhighlight lang="j"> eg1=: 25 15 _5 , 15 18 0 ,: _5 0 11
eg2=: 18 22 54 42 , 22 70 86 62 , 54 86 174 134 ,: 42 62 134 106
cholesky eg1
Line 1,902 ⟶ 2,357:
9.89949 1.62455 1.84971 1.39262</syntaxhighlight>
'''Using `math/lapack` addon'''
<syntaxhighlight lang="j"> load 'math/lapack'
load 'math/lapack/potrf'
potrf_jlapack_ eg1
Line 1,913 ⟶ 2,368:
12.7279 3.04604 1.64974 0
9.89949 1.62455 1.84971 1.39262</syntaxhighlight>
=={{header|Java}}==
{{works with|Java|1.5+}}
<syntaxhighlight lang="java5">import java.util.Arrays;
public class Cholesky {
Line 1,950 ⟶ 2,404:
<pre>[[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]]</pre>
=={{header|JavaScript}}==
<syntaxhighlight lang="javascript">
const cholesky = function (array) {
const zeros = [...Array(array.length)].map( _ => Array(array.length).fill(0));
Line 1,979 ⟶ 2,432:
3: (4) [9.899494936611665, 1.6245538642137891, 1.849711005231382, 1.3926212476455924]
</pre>
=={{header|jq}}==
{{Works with|jq|1.4}}
'''Infrastructure''':
<syntaxhighlight lang="jq"># Create an m x n matrix
def matrix(m; n; init):
if m == 0 then []
Line 2,021 ⟶ 2,473:
else [ ., 0, 0, length ] | test
end ;
</syntaxhighlight>'''Cholesky Decomposition''':<syntaxhighlight lang="jq">def cholesky_factor:
if is_symmetric then
length as $length
Line 2,053 ⟶ 2,505:
[42, 62, 134, 106]] | cholesky_factor | neatly(20)
{{Out}}
<syntaxhighlight lang="jq"> 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</syntaxhighlight>
=={{header|Julia}}==
Julia's strong linear algebra support includes Cholesky decomposition.
<syntaxhighlight lang=
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]
Line 2,087 ⟶ 2,538:
9.899494936611667 1.624553864213788 1.8497110052313648 1.3926212476456026]
</pre>
=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.0.6
fun cholesky(a: DoubleArray): DoubleArray {
Line 2,142 ⟶ 2,592:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Lobster}}==
{{trans|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.
<syntaxhighlight lang=
// choleskyLower returns the cholesky decomposition of a symmetric real
Line 2,247 ⟶ 2,696:
9.899494936612 1.624553864214 1.849711005231 1.392621247646
</pre>
=={{header|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.
<syntaxhighlight lang=
[25 15 -5]
[ ]
Line 2,302 ⟶ 2,750:
[ ]
[9.899494934 1.624553864 1.849711006 1.392621248]</syntaxhighlight>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang=
Without the use of built-in functions, making use of memoization:
<syntaxhighlight lang=
Module[{L},
L[k_, k_] := L[k, k] = Sqrt[A[[k, k]] - Sum[L[k, j]^2, {j, 1, k-1}]];
Line 2,312 ⟶ 2,759:
PadRight[Table[L[i, j], {i, Length[A]}, {j, i}]]
]</syntaxhighlight>
=={{header|MATLAB}} / {{header|Octave}}==
The cholesky decomposition chol() is an internal function
<syntaxhighlight lang=
25 15 -5
15 18 0
Line 2,344 ⟶ 2,790:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Maxima}}==
<syntaxhighlight lang="maxima">/* Cholesky decomposition is built-in */
a: hilbert_matrix(4)$
Line 2,358 ⟶ 2,803:
b . transpose(b) - a;
matrix([0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0])</syntaxhighlight>
=={{header|Nim}}==
{{trans|C}}
<syntaxhighlight lang="nim">import math, strutils, strformat
type Matrix[N: static int, T: SomeFloat] = array[N, array[N, T]]
Line 2,403 ⟶ 2,847:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|Objeck}}==
{{trans|C}}
<syntaxhighlight lang="objeck">
class Cholesky {
function : Main(args : String[]) ~ Nil {
Line 2,463 ⟶ 2,906:
9.89949494 1.62455386 1.84971101 1.39262125
</pre>
=={{header|OCaml}}==
<syntaxhighlight lang=
let n = Array.length inp in
let res = Array.make_matrix n n 0.0 in
Line 2,515 ⟶ 2,957:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|ooRexx}}==
{{trans|REXX}}
<syntaxhighlight lang="oorexx">/*REXX program performs the Cholesky decomposition on a square matrix. */
niner = '25 15 -5' , /*define a 3x3 matrix. */
'15 18 0' ,
Line 2,569 ⟶ 3,010:
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.*/</syntaxhighlight>
=={{header|PARI/GP}}==
<syntaxhighlight lang="parigp">cholesky(M) =
{
my (L = matrix(#M,#M));
Line 2,605 ⟶ 3,045:
[9.8995 1.6246 1.8497 1.3926]
</pre>
=={{header|Pascal}}==
<syntaxhighlight lang="pascal">program CholeskyApp;
type
Line 2,686 ⟶ 3,125:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Perl}}==
<syntaxhighlight lang="perl">sub cholesky {
my $matrix = shift;
my $chol = [ map { [(0) x @$matrix ] } @$matrix ];
Line 2,727 ⟶ 3,165:
9.8995 1.6246 1.8497 1.3926
</pre>
=={{header|Phix}}==
{{trans|Sidef}}
<!--<syntaxhighlight lang=
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cholesky</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">matrix</span><span style="color: #0000FF;">)</span>
Line 2,766 ⟶ 3,203:
{9.899494937,1.624553864,1.849711005,1.392621248}}
</pre>
=={{header|PicoLisp}}==
<syntaxhighlight lang=
(load "@lib/math.l")
Line 2,786 ⟶ 3,222:
(prinl) ) ) )</syntaxhighlight>
Test:
<syntaxhighlight lang=
'((25.0 15.0 -5.0) (15.0 18.0 0) (-5.0 0 11.0)) )
Line 2,806 ⟶ 3,242:
12.72792 3.04604 1.64974 0.00000
9.89949 1.62455 1.84971 1.39262</pre>
=={{header|PL/I}}==
<syntaxhighlight lang=
decompose: procedure options (main); /* 31 October 2013 */
declare a(*,*) float controlled;
Line 2,875 ⟶ 3,310:
9.89950 1.62455 1.84971 1.39262
</pre>
=={{header|PowerShell}}==
<syntaxhighlight lang=
function cholesky ($a) {
$l = @()
Line 2,954 ⟶ 3,388:
9.89949493661167 1.62455386421379 1.84971100523138 1.39262124764559
</pre>
=={{header|Python}}==
===Python2.X version===
<syntaxhighlight lang="python">from __future__ import print_function
from pprint import pprint
Line 2,997 ⟶ 3,430:
Factors out accesses to <code>A[i], L[i], and L[j]</code> by creating <code>Ai, Li and Lj</code> respectively as well as using <code>enumerate</code> instead of <code>range(len(some_array))</code>.
<syntaxhighlight lang="python">def cholesky(A):
L = [[0.0] * len(A) for _ in range(len(A))]
for i, (Ai, Li) in enumerate(zip(A, L)):
Line 3,008 ⟶ 3,441:
{{out}}
(As above)
=={{header|q}}==
<syntaxhighlight lang="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]}
Line 3,038 ⟶ 3,470:
9.899495 1.624554 1.849711 1.392621
</pre>
=={{header|R}}==
<syntaxhighlight lang="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
Line 3,052 ⟶ 3,483:
# [3,] 12.727922 3.046038 1.649742 0.000000
# [4,] 9.899495 1.624554 1.849711 1.392621</syntaxhighlight>
=={{header|Racket}}==
<syntaxhighlight lang="racket">
#lang racket
(require math)
Line 3,085 ⟶ 3,515:
</syntaxhighlight>
Output:
<syntaxhighlight lang="racket">
'#(#(5 0 0)
#(3 3 0)
Line 3,095 ⟶ 3,525:
</syntaxhighlight>
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>sub cholesky(@A) {
my @L = @A »×» 0;
for ^@A -> \i {
Line 3,117 ⟶ 3,546:
say '';
.fmt('%6.3f').say for cholesky [
[18, 22, 54, 42],
[22, 70, 86, 62],
Line 3,132 ⟶ 3,561:
12.728 3.046 1.650 0.000
9.899 1.625 1.850 1.393</pre>
=={{header|REXX}}==
If trailing zeros are wanted after the decimal point for values of zero (0), the <big><big>'''/ 1'''</big></big> (a division by unity performs
<br>REXX number normalization) can be removed from the line (number 40) which contains the source statement:
::::: <b> z=z right( format(@.row.col, , dPlaces) / 1, width) </b>
<syntaxhighlight lang="rexx">/*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' ,
Line 3,216 ⟶ 3,644:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Ring}}==
<syntaxhighlight lang="ring">
# Project : Cholesky decomposition
Line 3,273 ⟶ 3,700:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|Ruby}}==
<syntaxhighlight lang="ruby">require 'matrix'
class Matrix
def cholesky_factor
raise ArgumentError, "must provide symmetric matrix" unless symmetric?
Line 3,325 ⟶ 3,741:
{{trans|C}}
<syntaxhighlight lang="rust">fn cholesky(mat: Vec<f64>, n: usize) -> Vec<f64> {
let mut res = vec![0.0; mat.len()];
for i in 0..n {
Line 3,378 ⟶ 3,794:
9.8995 1.6246 1.8497 1.3926
</pre>
=={{header|Scala}}==
<syntaxhighlight lang="scala">case class Matrix( val matrix:Array[Array[Double]] ) {
// Assuming matrix is positive-definite, symmetric and not empty...
Line 3,441 ⟶ 3,856:
assert(math.round( result * 100000 ) * 0.00001 == math.round( test * 100000 ) * 0.00001)
}</syntaxhighlight>
=={{header|Scilab}}==
The Cholesky decomposition is builtin, and an upper triangular matrix is returned, such that $A=L^TL$.
<syntaxhighlight lang="scilab">a = [25 15 -5; 15 18 0; -5 0 11];
chol(a)
ans =
Line 3,465 ⟶ 3,879:
0. 0. 1.6497422 1.849711
0. 0. 0. 1.3926212</syntaxhighlight>
=={{header|Seed7}}==
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i";
include "float.s7i";
include "math.s7i";
Line 3,539 ⟶ 3,952:
9.89950 1.62455 1.84971 1.39262
</pre>
=={{header|Sidef}}==
{{trans|Perl}}
<syntaxhighlight lang="ruby">func cholesky(matrix) {
var chol = matrix.len.of { matrix.len.of(0) }
for row in ^matrix {
Line 3,557 ⟶ 3,969:
Examples:
<syntaxhighlight lang="ruby">var example1 = [ [ 25, 15, -5 ],
[ 15, 18, 0 ],
[ -5, 0, 11 ] ];
Line 3,588 ⟶ 4,000:
9.8995 1.6246 1.8497 1.3926
</pre>
=={{header|Smalltalk}}==
<syntaxhighlight lang=
FloatMatrix>>#cholesky
| l |
Line 3,614 ⟶ 4,025:
^l
</syntaxhighlight>
=={{header|Stata}}==
See [http://www.stata.com/help.cgi?mf_cholesky Cholesky square-root decomposition] in Stata help.
<syntaxhighlight lang="stata">mata
: a=25,15,-5\15,18,0\-5,0,11
Line 3,657 ⟶ 4,067:
4 | 9.899494937 1.624553864 1.849711005 1.392621248 |
+---------------------------------------------------------+</syntaxhighlight>
=={{header|Swift}}==
{{trans|Rust}}
<syntaxhighlight lang="swift">func cholesky(matrix: [Double], n: Int) -> [Double] {
var res = [Double](repeating: 0, count: matrix.count)
Line 3,721 ⟶ 4,130:
12.727922061357857 3.0460384954008553 1.6497422479090704 0.0
9.899494936611667 1.624553864213788 1.8497110052313648 1.3926212476456026</pre>
=={{header|Tcl}}==
{{trans|Java}}
<syntaxhighlight lang="tcl">proc cholesky a {
set m [llength $a]
set n [llength [lindex $a 0]]
Line 3,744 ⟶ 4,152:
}</syntaxhighlight>
Demonstration code:
<syntaxhighlight lang="tcl">set test1 {
{25 15 -5}
{15 18 0}
Line 3,762 ⟶ 4,170:
{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}
</pre>
=={{header|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 <code>Option Base 0</code> 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 <code>Byte</code>-type variables to <code>Long</code>. 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 "<code>=Cholesky(A1:E5)</code>" and htting <code>Ctrl-Shift-Enter</code> will populate the target cells with the lower Cholesky decomposition.
<syntaxhighlight lang="vb">Function Cholesky(Mat As Range) As Variant
Dim A() As Double, L() As Double, sum As Double, sum2 As Double
Line 3,819 ⟶ 4,226:
End Function
</syntaxhighlight>
=={{header|V (Vlang)}}==
{{trans|go}}
<syntaxhighlight lang="v (vlang)">import math
// Symmetric and Lower use a packed representation that stores only
Line 3,947 ⟶ 4,353:
{{libheader|Wren-matrix}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang=
import "./fmt" for Fmt
var arrays = [
Line 3,992 ⟶ 4,398:
|12.72792 3.04604 1.64974 0.00000|
| 9.89949 1.62455 1.84971 1.39262|
</pre>
=={{header|XPL0}}==
<syntaxhighlight lang "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);
]</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|zkl}}==
Using the GNU Scientific Library:
<syntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (GNU Scientific Library)
fcn lowerCholesky(m){ // trans: C
rows:=m.rows;
Line 4,032 ⟶ 4,497:
Or, using lists:
{{trans|C}}
<syntaxhighlight lang="zkl">fcn cholesky(mat){
rows:=mat.len();
r:=(0).pump(rows,List().write, (0).pump(rows,List,0.0).copy); // matrix of zeros
Line 4,042 ⟶ 4,507:
r
}</syntaxhighlight>
<syntaxhighlight lang="zkl">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("-----------------");
Line 4,050 ⟶ 4,515:
L(42.0, 62.0, 134.0, 106.0,) );
printM(cholesky(ex2));</syntaxhighlight>
<syntaxhighlight lang="zkl">fcn printM(m){ m.pump(Console.println,rowFmt) }
fcn rowFmt(row){ ("%9.5f "*row.len()).fmt(row.xplode()) }</syntaxhighlight>
{{out}}
Line 4,063 ⟶ 4,528:
9.89949 1.62455 1.84971 1.39262
</pre>
=={{header|ZX Spectrum Basic}}==
{{trans|BBC_BASIC}}
<syntaxhighlight lang="zxbasic">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
| |||