Python - Getting the wrong solution for Cholesky decomposition?-CodePudding

I'm trying to translate some pseudocode from matlab to a python script, but I'm having some trouble with getting the correct answer? Can someone help me identify where I messed up the translation?

The pseudocode I'm given in the one for Cholesky decomposition given in the Trefethan & Bau book is

This is done for an upper triangular matrix if i understand the description given correctly, but I think this should work for a general matrix too, no?

Anyway, I wrote the following code in python:

def is_SPD(A):
    if np.all(A == A.T):
        if np.all(la.eigvals(A) > 0):
            return True
    return False

def cholesky_decomp(A):
    #first check if matrix is symmetric and positive definite
    if is_SPD(A) == True:
        R = np.copy(A)
        for k in range (0, len(A)):
            for j in range (k 1, len(A)):
                R[j,j:] = R[j,j:] - (R[k,j:]*(R[k,j]/R[k,k]))
            R[k,k:] = R[k,k:]/sqrt(R[k,k])
        return R
    else:
        return print('Cholesky decomposition not applicable')

I'm working on a 4x4 matrix, and I've done the decomposition by the np.linalg method, and my answers are completely different.

I think this might? be because of my unfamiliarity with MATLAB and my lack of coding skills in general, but I'm not getting any correct answers at all and I can't figure out where I'm going wrong.

I'm adding a sample matrix here that I'm using, and that this shoulddd be applicable to, and should give a proper Cholesky decomposition for, but I'm getting a completely incorrect answer.

Could someone please use this to help me figure out where I'm going wrong?

A = np.array([[16, -12, -12, -16], [-12, 25, 1, -4], [-12, 1, 17, 14], [-16, -4, 14, 57]])

my code gave me:

[[  4  -3  -3  -4]
 [-12   4  -2  -4]
 [-12   1   2  -3]
 [-16  -4  14   4]]

while the numpy Cholesky function gave me:

[[ 4.  0.  0.  0.]
 [-3.  4.  0.  0.]
 [-3. -2.  2.  0.]
 [-4. -4. -3.  4.]]

CodePudding user response：

Your code correctly implements the stated algorithm, but note that the text says (added emphasis)

The input matrix A represents the superdiagonal half of the m×m Hermitian positive definite matrix to be factored.

So you need to replace the input A,

[[ 16 -12 -12 -16]
 [-12  25   1  -4]
 [-12   1  17  14]
 [-16  -4  14  57]]

by np.triu(A):

[[ 16 -12 -12 -16]
 [  0  25   1  -4]
 [  0   0  17  14]
 [  0   0   0  57]]

In addition (notation slightly changed, where ' means Hermitian transpose),

The output matrix R represents the upper-triangular factor for which A = R' R

So you get an upper-triangular R, whereas Numpy's cholesky function gives a lower-triangular result, which is just the Hermitian-transposed version of yours (see here).