Find out if matrix is positive definite with numpy

I need to find out if matrix is positive definite. My matrix is numpy matrix. I was expecting to find any related method in numpy library, but no success. I appreciate any help.

You can also check if all the eigenvalues of matrix are positive, if so the matrix is positive definite:

import numpy as np

def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 0)

You could try computing Cholesky decomposition (numpy.linalg.cholesky). This will raise LinAlgError if the matrix is not positive definite.

There seems to be a small confusion in all of the answers above (at least concerning the question).

For real matrices, the tests for positive eigenvalues and positive-leading terms in np.linalg.cholesky only applies if the matrix is symmetric. So first one needs to test if the matrix is symmetric and then apply one of those methods (positive eigenvalues or Cholesky decomposition).

For example:

import numpy as np

#A nonsymmetric matrix
A = np.array([[9,7],[6,14]])

#check that all eigenvalues are positive:
np.all(np.linalg.eigvals(A) > 0)

#take a 'Cholesky' decomposition:
chol_A = np.linalg.cholesky(A)

The matrix A is not symmetric, but the eigenvalues are positive and Numpy returns a Cholesky decomposition that is wrong. You can check that:

chol_A.dot(chol_A.T)

is different than A.

You can also check that all the python functions above would test positive for 'positive-definiteness'. This could potentially be a serious problem if you were trying to use the Cholesky decomposition to compute the inverse, since:

>np.linalg.inv(A)
array([[ 0.16666667, -0.08333333],
   [-0.07142857,  0.10714286]])

>np.linalg.inv(chol_A.T).dot(np.linalg.inv(chol_A))
array([[ 0.15555556, -0.06666667],
   [-0.06666667,  0.1       ]])

are different.

In summary, I would suggest adding a line to any of the functions above to check if the matrix is symmetric, for example:

def is_pos_def(A):
    if np.array_equal(A, A.T):
        try:
            np.linalg.cholesky(A)
            return True
        except np.linalg.LinAlgError:
            return False
    else:
        return False

You may want to replace np.array_equal(A, A.T) in the function above for np.allclose(A, A.T) to avoid differences that are due to floating point errors.

To illustrate @NPE's answer with some ready-to-use code:

import numpy as np

def is_pd(K):
    try:
        np.linalg.cholesky(K)
        return 1 
    except np.linalg.linalg.LinAlgError as err:
        if 'Matrix is not positive definite' in err.message:
            return 0
        else:
            raise 

I don't know why the solution of NPE is so underrated. It's the best way to do this. I've found on Wkipedia that the complexity is cubic.

Furthermore, there it is said that it's more numerically stable than the Lu decomposition. And the Lu decomposition is more stable than the method of finding all the eigenvalues.

And, it is a very elegant solution, because it's a fact :

A matrix has a Cholesky decomposition if and only if it is symmetric positive.

So why not using maths ? Maybe some people are affraid of the raise of the exception, but it'a fact too, it's quite useful to program with exceptions.

If you specifically want symmetric (hermitian, if complex) positive SEMI-definite matrices than the below will do. If you don't care about symmetry (hermitian, if complex) remove the 'if' state that checks for it. If you want positive definite rather than positive SEMI-definite than remove the regularization line (and change the value passed to 'np.lingalg.cholesky()' from 'regularized_X' to 'X'). The below

import numpy as np

def is_hermitian_positive_semidefinite(X):
    if X.shape[0] != X.shape[1]: # must be a square matrix
        return False

    if not np.all( X - X.H == 0 ): # must be a symmetric or hermitian matrix
        return False

    try: # Cholesky decomposition fails for matrices that are NOT positive definite.

        # But since the matrix may be positive SEMI-definite due to rank deficiency
        # we must regularize.
        regularized_X = X + np.eye(X.shape[0]) * 1e-14

        np.linalg.cholesky(regularized_X)
    except np.linalg.LinAlgError:
        return False

    return True

For a real matrix $A$, we have $x^TAx=\frac{1}{2}(x^T(A+A^T)x)$, and $A+A^T$ is symmetric real matrix. So $A$ is positive definite iff $A+A^T$ is positive definite, iff all the eigenvalues of $A+A^T$ are positive.

import numpy as np

def is_pos_def(A):
    M = np.matrix(A)
    return np.all(np.linalg.eigvals(M+M.transpose()) > 0)

For Not symmetric Matrix you can use the Principal Minor Test :

def isPD(Y):
  row = X.shape [0]
  i = 0
  j = 0
  for i in range(row+1) : 
    Step = Y[:i,:j]
    j+=1
    i+=1
    det = np.linalg.det(Step)
    if det > 0 : 
        continue 
    else :
        return ("Not Positive Definite, Test Principal minor failed")

  return ("Positive Definite")

Positive Definite

numpy.linalg.cholesky(x) # just handle the error LinAlgError

Positive Semi-Definite

np.all(np.linalg.eigvals(x) >= 0)

Notice: Most of cases also np.all(np.linalg.eigvals(x) > 0) will give you if your matrix is PSD even if you see > and not only >=, I got into this problem some days ago. I figure it should be something about round off error due to the fact that we have really small eigenvalues and even cholesky decomposition might generate an error.

Note

In order to test, you might want to create some Positive semi-definite matrix and some positive definite matrices:

n_size=4
a = np.random.rand(n_size)
A_PSD = np.outer(a,a)  # the outer product of any vector generates a PSD matrix
A_PD = A_PSD+np.identity(n_size) # little trick I found for PS matrix