python - semidefinite - Find out if matrix is positive definite with numpy




random positive semidefinite matrix numpy (4)

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.


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)

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.


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 

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)






scipy