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

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+1e-5*np.identity(n_size) # little trick I found for PS matrix

Answer 10

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

Answer 11

Numerical inaccuracy is a problem when determining whether a matrix is positive/negative -definite or semi-definite.

To illustrate this consider the following piece of code where I show that a positive-semi-definite matrix can look positive-definite or even indefinite due to numerical inaccuracy:

import numpy as np

np.random.seed(1234)


n = 5

A = np.random.uniform(-1, 1, (n, n))
B = A @ A.T  # guaranteed to be positive-definite

w, v = np.linalg.eigh(B)
print(w)
# [0.00571615 0.29977027 0.44425488 2.67023429 3.82585751]

C = v @ np.diag(w) @ v.T
w = np.linalg.eigvalsh(C)
print(w)
# [0.00571615 0.29977027 0.44425488 2.67023429 3.82585751]
# Same as original!

w[0] = 0  # set smallest two eigen-values to zero
E = v @ np.diag(w) @ v.T
w = np.linalg.eigvalsh(E)
print(w)
# [4.06224137e-16 2.99770273e-01 4.44254883e-01 2.67023429e+00 3.82585751e+00]
#  ^^^^^^^^^^^^^^
# Ooops! matrix is now positive-definite?!

w[:2] = 0  # set smallest two eigen-values to zero
E = v @ np.diag(w) @ v.T
w = np.linalg.eigvalsh(E)
print(w)
# [-5.14560681e-16  4.27158960e-17  4.44254883e-01  2.67023429e+00 3.82585751e+00]
#  ^^^^^^^^^^^^^^^
# Ooops! matrix is not definite?!

The code above consists of the following steps:

1. Construct a random matrix B which is guaranteed to be positive-definite.
2. Extract eigen-values/vectors, reconstruct matrix from these and extract eigen-values again to show that they coincide with the original eigen-values.
3. Make the smallest eigen-value zero, construct the resulting positive-semi-definite matrix, extract its eigen-values and show that it looks positive-definite.
4. Make the two smallest eigen-value zero, construct the resulting positive-semi-definite matrix, extract its eigen-values and show that it looks indefinite.

This example shows that even if the matrix is known to be positive-semi-definite we might not get any eigen-values which are exactly zero and we might even get very small negative eigen-values making the matrix look indefinite.

To overcome this issue we need to include a tolerance with which we will determine whether eigen-values are positive/negative/zero.

Here is an example which shows that by using a tolerance we can correctly identify eigen-values that should have been zero and determine whether the matrix is positive/negative definite or semi-definite.:

tol = np.finfo(float).eps

w_abs = np.abs(w)

is_zero = w_abs < w_abs.max() * tol
print(is_zero)
# [ True  True False False False]
#   ^^^^  ^^^^
# first two (smallest) eigen-values correctly identified as being zero

is_positive = w[~is_zero] > 0

s = 'the matrix is '
if is_positive.all():
    s += 'positive-'
elif not is_positive.any():
    s += 'negative-'
if is_zero.any():
    s += 'semi-'
s += 'definite'

print(s)
# the matrix is positive-semi-definite

Including a tolerance as in the example above primarily mitigates the issue of confusing a (semi-)definite matrix for an indefinite matrix, but there is still the issue that a definite matrix may well be misidentified as being semi-definite if some of the (strictly positive/negative) eigen-values are very small. I don't know of any robust solution for this and I suspect that it might just be an intrinsic limitation of finite precision arithmetic.

I have written the following simple function which determines whether a matrix is positive/negative definite or semi-definite with a tolerance:

def is_definite(x, sign=None, semi=None, tol=None):
    """
    Determine whether a matrix is positive/negative definite or semi-definite.

    Examples:

    >>> import numpy as np
    >>> np.random.seed(1234)

    >>> x = np.random.uniform(-1, 1, (5, 5)) + 1j * np.random.uniform(-1, 1, (5, 5))
    >>> is_definite(x)
    False

    >>> x = x @ x.conj().T  # x @ x.conj().T is guaranteed to be positive-definite
    >>> is_definite(x)
    True
    >>> is_definite(x, -1)
    False
    >>> is_definite(x, 1)
    True

    >>> w, v = np.linalg.eigh(x)
    >>> w[0] = 0
    >>> w
    array([0.        , 0.2272762 , 1.46465277, 4.61979679, 8.14691898])
    >>> x = v @ np.diag(w) @ v.conj().T
    >>> np.linalg.eigvalsh(x)
    array([-3.18173864e-16,  2.27276198e-01,  1.46465277e+00,  4.61979679e+00,  8.14691898e+00])
    >>> is_definite(x, 1)
    False
    >>> is_definite(x, -1)
    False
    >>> is_definite(x, 1, semi=True)
    True

    :param x: M*M-matrix or array of M*M-matrices
    :param sign:
            positive or negative number to check for positive or negative (semi-)definiteness respectively.
            Default is to return True for any definite matrix (positive or negative).
    :param semi: whether to check for semi-definiteness. Default is False
    :param tol: tolerance, default is machine precision
    :return: bool or array of bools
    """

    x = np.asarray(x)

    if x.ndim < 2:
        raise ValueError('x must be at least two-dimensional')

    x = (x + np.moveaxis(x, -1, -2).conj())/2

    if tol is None:
        tol = np.finfo(x.dtype).eps

    w = np.linalg.eigvalsh(x)
    w_abs = np.abs(w)
    is_zero = w_abs < w_abs.max(-1, keepdims=True) * tol

    if sign is None:
        pos = w > 0
        neg = w < 0
        if semi:
            pos |= is_zero
            neg |= is_zero
        else:
            pos &= ~is_zero
            neg &= ~is_zero
        return pos.all(-1) | neg.all(-1)

    ret = np.sign(w) == np.sign(sign)
    if semi:
        ret |= is_zero
    else:
        ret &= ~is_zero
    return ret.all(-1)


def is_positive_definite(x, tol=None):
    return is_definite(x, 1, False, tol)


def is_positive_semidefinite(x, tol=None):
    return is_definite(x, 1, True, tol)


def is_negative_definite(x, tol=None):
    return is_definite(x, -1, False, tol)


def is_negative_semidefinite(x, tol=None):
    return is_definite(x, -1, True, tol)


def is_indefinite(x, tol=None):
    return not is_definite(x, None, True, tol)