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.
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)
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
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:
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.
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 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")
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
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 zT · 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)
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:
Construct a random matrix B which is guaranteed to be positive-definite.
Extract eigen-values/vectors, reconstruct matrix from these and extract eigen-values again to show that they coincide with the original eigen-values.
Make the smallest eigen-value zero, construct the resulting positive-semi-definite matrix, extract its eigen-values and show that it looks positive-definite.
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)
