Check for positive definiteness or positive semidefiniteness

Cholesky decomposition is a good option if you're working with positive definite (PD) matrices.

However, it throws the following error on positive semi-definite (PSD) matrix, say,

A = np.zeros((3,3)) // the all-zero matrix is a PSD matrix
np.linalg.cholesky(A)
LinAlgError: Matrix is not positive definite -
Cholesky decomposition cannot be computed

For PSD matrices, you can use scipy/numpy's eigh() to check that all eigenvalues are non-negative.

>> E,V = scipy.linalg.eigh(np.zeros((3,3)))
>> E
array([ 0.,  0.,  0.])

However, you will most probably encounter numerical stability issues. To overcome those, you can use the following function.

def isPSD(A, tol=1e-8):
  E = np.linalg.eigvalsh(A)
  return np.all(E > -tol)

Which returns True on matrices that are approximately PSD up to a given tolerance.

I assume you already know your matrix is symmetric.

A good test for positive definiteness (actually the standard one !) is to try to compute its Cholesky factorization. It succeeds iff your matrix is positive definite.

This is the most direct way, since it needs O(n^3) operations (with a small constant), and you would need at least n matrix-vector multiplications to test "directly".

Check whether the whole eigenvalues of a symmetric matrix A are non-negative is time-consuming if A is very large, while the module scipy.sparse.linalg.arpack provides a good solution since one can customize the returned eigenvalues by specifying parameters.(see Scipy.sparse.linalg.arpack for more information)

As we know if both ends of the spectrum of A are non-negative, then the rest eigenvalues must also be non-negative. So we can do like this:

from scipy.sparse.linalg import arpack
def isPSD(A, tol = 1e-8):
    vals, vecs = arpack.eigsh(A, k = 2, which = 'BE') # return the ends of spectrum of A
    return np.all(vals > -tol)

By this we only need to calculate two eigenvalues to check PSD, I think it's very useful for large A