Cholesky decomposition positive semidefinite matrix

Question

I have read that the Cholesky decomposition of a matrix in Numpy/Scipy only works if it's positive definite. Indeed, the following doesn't work, as the matrix is positive semi-definite

np.linalg.cholesky([[1, 0], [0, 0]])
numpy.linalg.linalg.LinAlgError: Matrix is not positive definite

However, I am using a symmetric positive semi-definite matrix, for which Numpy does the decomposition without error:

np.linalg.cholesky([[2, 6], [6, 18]])
array([[1.41421356e+00, 0.00000000e+00],
      [4.24264069e+00, 5.64928468e-08]])

What is happening? Both matrices in my tests are positive semi-definite, but from what I read, I expected the Cholesky decomposition in Numpy/Scipy to only work with positive definite matrices and give a LinAlgError, otherwise.

Martin Ueding Martin Ueding · Accepted Answer · 2019-03-14T08:50:54

You can have a look whether the decomposition has worked. By multiplying your matrix L with its transpose, you should get your initial matrix back:

>>> import numpy as np
>>> L = np.array([[1.41421356e+00, 0.00000000e+00],
...       [4.24264069e+00, 5.64928468e-08]])
>>> L.dot(L.T)
array([[ 1.99999999,  5.99999999],
       [ 5.99999999, 18.00000002]])

The eigenvalues of your matrix are 0 and 20, so it is positive semi-definite.

I would conclude that this algorithm might work, but it is not guaranteed that the algorithm gives sensible results. Similarly, you can use the conjugate gradient to solve a system of equations even though the system is not a hermitian matrix. Although the conjugate gradient sometimes works, it might not work as well.

Perhaps the checking criterion used in NumPy does not have to be too strict because it might be able to cope with certain semi-definite matrices.

Cholesky decomposition positive semidefinite matrix

1 Answers