Find out if matrix is positive definite with numpy

Question

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.

Answer 1

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)

Answer 2

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

Answer 3

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.

Answer 4

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 

Answer 5

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.

Answer 6

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

Answer 7

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)

Answer 8

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")

Answer 9

May I suggest this solution, valid for non-symmetric matrices. Basically it tries to find a non-zero vector z that minimizes the result of z^T · M · z. As you see, it can either converge in a minimum (minV['success']) or reach the max number of iteration (minV['status'] == 2). If at this point the result is still positive, we could consider that the matrix is positive definite.

I guess there must be analytical methods to determine it, but I see a bit of confusion here about symmetric matrices (not a premise to be positive definite!!).

from scipy.optimize import minimize
import numpy as np

def is_pos_def(m,z0='rand', maxiter=100):
    #tells if matrix is positive definite
    if m.shape[0] != m.shape[1]:
        raise Exception("Matrix is not square") 
    elif (m==m.T).all(): #symmetry testing
        return np.all(np.linalg.eigvals(m) > 0)
    else:
        def f(z):
            z=np.array(list(z))[:,np.newaxis]
            return np.dot(np.dot(z.T, m),z)
        if z0=='rand':
            z0 = list(np.random.rand(m.shape[0]))
        #constraints for a non-zero vector solution
        cons = ({'type': 'ineq', 'fun': lambda z:  np.sum(np.abs(z))})
        minV = minimize(f, z0, method='COBYLA', options={'maxiter' : maxiter},constraints=cons);

        if minV['success'] or minV['status'] == 2:
            return minV['fun']+0 > 0 
        else:        
            return minV

This method works for both symmetric and non-symmetric, you can test with the following matrix (checked with wolfram alpha too)

m=np.array([[3, 0, 0],
            [0, 2, 0],
            [4, 3, 3]])