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.

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

Solution 2

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

Solution 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.

Solution 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 

Solution 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.

View more solutions

Comments

[jorgeca]

You could use np.linalg.eigvals instead, which only computes the eigenvalues. Even then, it's much slower than @NPE's approach (3x for 10x10 matrices, 40x for 1000x1000).

[Warren Weckesser]

<pedantic>It is not true in general that all positive eigenvalues implies positive definiteness, unless you know that the matrix is symmetric (real case) or Hermitian (complex case). For example, A = array([[1, -100],[0, 2]]) is not positive definite. Some might include symmetric or Hermitian as part of the definition of "positive definite", but that is not universal.</pedantic>

[jorgeca]

@WarrenWeckesser Oops, that's right, not pedantic! In fact, checking symmetry is also needed if using np.linalg.cholesky (it doesn't check it and may return a wrong result, as your example also shows). I wonder how checking whether a non symmetric matrix is positive definite can be done numerically...

[MRocklin]

This should be substantially more efficient than the eigenvalue solution.

[shinjin]

You can do np.all(x-x.T==0) to check for symmetry

[Alex Flint]

This is terribly inefficient! For matrices larger than about 6 or 7 rows/columns, use cholesky as pointed out by NPE below. The cholesky route feels less convenient (catching an exception etc) but it is much less wasteful.

[user3731622]

From the same Wikipedia page, it seems like your statement is wrong. The page says " If the matrix A is Hermitian and positive semi-definite, then it still has a decomposition of the form A = LL* if the diagonal entries of L are allowed to be zero.[3]" Thus a matrix with a Cholesky decomposition does not imply the matrix is symmetric positive definite since it could just be semi-definite. Am I interpreting this wrong? Also, it seems like you've just thrown "symmetric" across the implication. i.e. shouldn't it be every Hermitian positive-definite matrix has unique Cholesky decomposition?

[jawknee]

Just a note that in the positive semi-definite case, numerically speaking, one can also add a little identity to the matrix (thus shifting all eigenvalues a small amount e.g. a few times machine precision) then use the cholesky method as usual.

[Imperishable Night]

Nitpick: you should probably check for numpy.linalg.LinAlgError unless you import it with from numpy.linalg import LinAlgError, which isn't a thing I would want to do if I would only catch this specific exception once or twice in my code.

[H. Rev.]

If working with complex matrices, this might lead to error (namely if A is complex positive definite, hence hermitian with strictly positive eigenvalues, the cholesky trick is still correct but it will not pass the first if statement for np.allclose(A, A.T)). Using np.allclose(A, A.H) will fix this (the OP says he is working with numpy matrix, if working with ndarray, use A.conj().T instead of .H

[Guillaume]

This is the only answer properly answering the question by OP : "how to determine if a matrix is DP". All the other answers confusingly make the assumption that symmetry is needed for a matrix to be definite positive, which is not the case.

[CognizantApe]

@DeepRazi Numpy's Cholesky decomposition implementation works on complex numbers (i.e. complex np.dtype). So yes it works in that sense. But my code above originally checked if the transpose rather than the conjugate-transpose is equal to itself which makes the overall function invalid for complex numbers. I have now change the transpose to conjugate-transpose and it is now valid for complex numbers.