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

4

2 Answers

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. Each of these results is just the Hermitian-transposed version of the other (see here).

1

L=numpy.linalg.cholesky(A)

A=np.dot(L,L.T.conj())

1

Your Answer

Sign up or log in

Sign up using Google Sign up using Facebook Sign up using Email and Password Submit

Post as a guest

Post Your Answer Discard

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy