# Cholesky simple example

**URL:** <https://discourse.julialang.org/t/cholesky-simple-example/69480>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [October 9, 2021, 6:19pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480 "2021-10-09T18:19:48Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)\
**Post date:** [October 9, 2021, 6:19pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/1 "2021-10-09T18:19:48Z")

</div>

Hello I try to implement from basis (I need it for GPU kernel (just small 3 by 3 matrix one of multiple things that a kernel do - no point in using cublas …) so -can not using Linear algebra package for a moment , also for some unclear reason the julia linear algebra implementation is using inversion which purpose is hidden from me- and from my perspective much complicating thing)  
Hence I try to adapt python code from  
[Cholesky Decomposition : Matrix Decomposition - GeeksforGeeks](https://www.geeksforgeeks.org/cholesky-decomposition-matrix-decomposition/)

and I know that it is very simple task, yet for some reason I am stuck for hours and I can not make it work

python working code

```julia
# Python3 program to decompose

# a matrix using Cholesky

# Decomposition

import math

MAX = 100;

 

def Cholesky_Decomposition(matrix, n):

 

    lower = [[0 for x in range(n + 1)]

                for y in range(n + 1)];

 

    # Decomposing a matrix

    # into Lower Triangular

    for i in range(n):

        for j in range(i + 1):

            sum1 = 0;

 

            # summation for diagonals

            if (j == i):

                for k in range(j):

                    sum1 += pow(lower[j][k], 2);

                lower[j][j] = int(math.sqrt(matrix[j][j] - sum1));

            else:

                 

                # Evaluating L(i, j)

                # using L(j, j)

                for k in range(j):

                    sum1 += (lower[i][k] *lower[j][k]);

                if(lower[j][j] > 0):

                    lower[i][j] = int((matrix[i][j] - sum1) /

                                               lower[j][j]);

 

    # Displaying Lower Triangular

    # and its Transpose

    print("Lower Triangular\t\tTranspose");

    for i in range(n):

         

        # Lower Triangular

        for j in range(n):

            print(lower[i][j], end = "\t");

        print("", end = "\t");

         

        # Transpose of

        # Lower Triangular

        for j in range(n):

            print(lower[j][i], end = "\t");

        print("");

 

# Driver Code

n = 3;

matrix = [[4, 12, -16],

          [12, 37, -43],

          [-16, -43, 98]];

Cholesky_Decomposition(matrix, n);

 

# This code is contributed by mits

```

my Julia implementation that do not give correct results

```julia
using LinearAlgebra

A = [6 15 55;15 55 225;55 225 979]
L = zeros(3,3)

L = zeros(3,3)

lower =L
matrix = A
n=3

for i in 0:(n-1)

    for j in 0:i

        sum1 = 0;

        # summation for diagonals

        if (j == i)

            for k in 0:j-1

            sum1 += lower[j+1,k+1]*lower[j+1,k+1]

            lower[j+1,j+1] = sqrt(matrix[j+1,j+1] - sum1)

            end#for

        else

             

            # Evaluating L(i, j)

            # using L(j, j)

            for k in 0:j-1

                sum1 += (lower[i+1,k+1] *lower[j+1,k+1]);

                if(lower[j+1,j+1] > 0)

                    lower[i+1,j+1] = (matrix[i+1,j+1] - sum1) /lower[j+1,j+1]

                end#if

            end#for

        end#if else

    end#for

end#for
#this gives diffrent results
L
cholesky([6 15 55;15 55 225;55 225 979])

```

Where is my mistake?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 9, 2021, 8:08pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/2 "2021-10-09T20:08:43Z")

</div>

Can you provide a runnable example of the complete code? And best yet with what you expect to be the correct solution?

Take a look at: [Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [October 9, 2021, 8:13pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/3 "2021-10-09T20:13:58Z")

</div>

> [@Jakub\_Mitura](#):
>
> ```julia
> matrix = [[4, 12, -16],
> [12, 37, -43],
> [-16, -43, 98]];
> 
> ```

Just in case, this is not a matrix in Julia but a 3-element Vector{Vector{Int64}}.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [October 9, 2021, 8:40pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/4 "2021-10-09T20:40:21Z")

</div>

_ **NB:** there seems to be a very old entry on this topic in [Rosetta code](https://rosettacode.org/wiki/Cholesky_decomposition#Julia), that a competent Julia user might want to update._

---

<div class="post-metadata">

**Author:** ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)\
**Post date:** [October 10, 2021, 4:15am UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/5 "2021-10-10T04:15:14Z")

</div>

You are right! I updated the code

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.julialang.org/u/jroon)\
**Post date:** [October 10, 2021, 7:08am UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/6 "2021-10-10T07:08:40Z")

</div>

Hi. I found two small errors.

The first is indexing k, it should be ` for k in 0:j` that is without the -1.  
The second is the location of the `lower[j,j] =` statement. It should not be inside the `for k ` loop. Remember Python code is indentation sensitive (yuck!) - so in these Python lines:

```julia
for k in range(j):
    sum1 += pow(lower[j][k], 2);
lower[j][j] = int(math.sqrt(matrix[j][j] - sum1));

```

→ the `lower[j][j]` is **after** the for loop not within it.

Here is your code modified to fix above and removing the Python-esque indexing:

```julia
for i in 1:n
    for j in 1:i
        sum1 = 0;
        # summation for diagonals
        if (j == i)
            for k in 1:j
                sum1 += (lower[j,k])^2
            end#for
            lower[j,j] = sqrt(matrix[j,j] - sum1)
        else
            # Evaluating L(i, j) using L(j, j)
            for k in 1:j
                sum1 += (lower[i,k] * lower[j,k]);
            end#for
            if(lower[j,j] > 0)
                lower[i,j] = (matrix[i,j] - sum1) /lower[j,j]
            end#if
        end#if else
    end#for
end#for

```

Running this and comparing to built in `cholesky()`:

```julia
julia> L
3×3 Matrix{Float64}:
  2.44949 0.0 0.0
  6.12372 4.1833 0.0
 22.4537 20.9165 6.1101

julia> cholesky(A).L
3×3 LowerTriangular{Float64, Matrix{Float64}}:
  2.44949 ⋅ ⋅ 
  6.12372 4.1833 ⋅ 
 22.4537 20.9165 6.1101

```

---

<div class="post-metadata">

**Author:** ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)\
**Post date:** [October 10, 2021, 10:26am UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/7 "2021-10-10T10:26:18Z")

</div>

Huge thanks for Help !

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 10, 2021, 10:32am UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/8 "2021-10-10T10:32:06Z")

</div>

Be aware that running the code like that will be slow. Make of your computations functions. See the [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/)

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.julialang.org/u/jroon)\
**Post date:** [October 10, 2021, 2:34pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/9 "2021-10-10T14:34:08Z")

</div>

> [@Jakub\_Mitura](#):
>
> Huge thanks for Help !

It actually helped me too I was about to google how to code this when I saw your question - you saved me some googling 😅

> [@lmiq](#):
>
> Be aware that running the code like that will be slow. Make of your computations functions. See the [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/)

Good point, so I made it into a function:

```julia
function chol(A)
    n = size(A)[1] # note I didn't check it is square here
    lower = zeros(n,n)

    for i in 1:n
        for j in 1:i
            sum1 = 0;
            # summation for diagonals
            if (j == i)
                for k in 1:j
                    sum1 += (lower[j,k])^2
                end#for
                lower[j,j] = sqrt(matrix[j,j] - sum1)
            else 
                # Evaluating L(i, j)
                # using L(j, j)
                for k in 1:j
                    sum1 += (lower[i,k] * lower[j,k]);
                end#for
                if(lower[j,j] > 0)
                    lower[i,j] = (matrix[i,j] - sum1) /lower[j,j]
                end#if
            end#if else
        end#for
    end#for
    lower
end

```

Timings versus built in `cholesky`:

```julia
julia> @btime cholesky($A).L;
  359.445 ns (5 allocations: 384 bytes)

julia> @btime chol($A);
  949.684 ns (23 allocations: 512 bytes)

```

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [October 10, 2021, 3:05pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/10 "2021-10-10T15:05:04Z")

</div>

Actually, there is a non-exported function `LinearAlgebra.checksquare(M)` which returns a common size if `M` is square and throws an exception otherwise.

And I think you’d want the inner loop go along columns for better performance. In that case, you’ll get an upper triangular matrix but that’s what the built-in `cholesky` gets for `Matrix` type, too. `cholesky($A).U'` must be faster than `cholesky($A).L` for that reason, as the `L` matrix is not stored in a `Cholesky` object but instead allocated and filled on demand.

For consistency, it’s also better to return `Cholesky(lower, 'L', 0)` rather than just the factor.

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.julialang.org/u/jroon)\
**Post date:** [October 10, 2021, 3:20pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/11 "2021-10-10T15:20:48Z")

</div>

Great tips thsnk @Vasily_Pisarev

---

<div class="post-metadata">

**Author:** ![Jakub\_Mitura](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_mitura/32/19496_2.png) [@Jakub\_Mitura](https://discourse.julialang.org/u/Jakub_Mitura)\
**Post date:** [October 10, 2021, 6:19pm UTC](https://discourse.julialang.org/t/cholesky-simple-example/69480/12 "2021-10-10T18:19:20Z")

</div>

thanks ! also to share it in case somebody will need it the unrolled equation for 3 by 3 covariance matrix looks like that  
unrolled 3 by 3 cholesky decomposition than forward substitution in order to calculate Mahalanobis distance according to forumula from [computational statistics - Efficient/fast Mahalanobis distance computation - Cross Validated (stackexchange.com)](https://stats.stackexchange.com/questions/147210/efficient-fast-mahalanobis-distance-computation/147222#147222?newreg=a68aa51b2f8c45daaece49163105845c)

```julia
a = sqrt(varianceX)

b = (covarianceXY)/a

c = (covarianceXZ)/a

e = sqrt(varianceY - b*b)

d = (covarianceYZ -(c * b))/e

#unrolled forward substitiution

ya= x[1]/a

yb = (x[2]-b*ya)/e

yc= (x[3]-yb*d-ya* c)/sqrt(varianceZ - c*c -d*d )

#taking square euclidean distance

return ya*ya+yb*yb+yc*yc

```

I got to this by printing from loop as you specified above and then simplifying … , so thanks !! this is huge speed up