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?

Take a look at: Please read: make it easier to help you

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

Good point, so I made it into a function:

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> @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.

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