# Cholesky simple example

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

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

``````# 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

``````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?

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

2 Likes

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

3 Likes

NB: there seems to be a very old entry on this topic in Rosetta code, that a competent Julia user might want to update.

1 Like

You are right! I updated the code

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:

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

``````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> 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
``````
4 Likes

Huge thanks for Help !

1 Like

Be aware that running the code like that will be slow. Make of your computations functions. See the Performance Tips · The Julia Language

2 Likes

It actually helped me too I was about to google how to code this when I saw your question - you saved me some googling Good point, so I made it into a function:

``````function chol(A)
n = size(A) # 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> @btime cholesky(\$A).L;
359.445 ns (5 allocations: 384 bytes)

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

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.

2 Likes

Great tips thsnk @Vasily_Pisarev

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)

``````a = sqrt(varianceX)

b = (covarianceXY)/a

c = (covarianceXZ)/a

e = sqrt(varianceY - b*b)

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

#unrolled forward substitiution

ya= x/a

yb = (x-b*ya)/e

yc= (x-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