Skip to content
Advertisement

Cholesky factorisation not lower triangular

I am building a cholesky factorisation algorthim as proposed from the book:

Linear Algebra and Optimisation for Machine Learning

They provide the following algorthim:

enter image description here

I have attempted this with python with the following algorithm:

def CF(array):
    array = np.array(array, float)
    arr = np.linalg.eigvals(array)
    for o in arr:
        if o < 0:
            return print('Cannot apply Cholesky Factorisation')
        else:
            continue
            
    L = np.zeros_like(array)
    n,_ = array.shape
    for j in range(n):
        sumkj = 0
        for k in range(j):
            sumkj += L[j,k]**2
            print(array[j,j], sumkj)
        L[j,j] = np.sqrt(array[j,j] - sumkj)
        for i in range(j+1):
            sumki = 0
            for ki in range(i):
                sumki += L[i,ki]
            L[i,j] = (array[i,j] - sumki*sumkj)/L[j,j]
    return L

When testing this on a matrix I get upper triangular as opposed to lower triangular. Where was I mistaken?

a = np.array([[4, 2, 1], [2, 6, 1], [2, 8, 2]])

CF(a)
>array([[2.        , 0.81649658, 0.70710678],
       [0.        , 2.44948974, 0.70710678],
       [0.        , 0.        , 1.41421356]])

Advertisement

Answer

You don’t even need to look at the detail of the algorithm.

Just see that in your book algorithm you have

for j=1 to d
    L(j,j) = something
    for i=j+1 to d
        L(i,j) = something

So, necessarily, all elements whose line number i is greater or equal than column number j are filled. Rest is 0. Hence a lower triangle

In your algorithm on another hand you have

for j in range(n):
    L[j,j] = something()
    for i in range(j+1):
        L[i,j] = something()

So all elements whose line number i is < to column number j, +1 (since i is in range(j+1), that is i<j+1, that is i<=j) are filled with something. Rest is 0.

Hence upper triangle

You probably wanted to

    for i in range(j+1, n)
User contributions licensed under: CC BY-SA
6 People found this is helpful
Advertisement