Testing if a matrix is positive semi-definite

I have a list ${\cal L}$ of symmetric matrices that I need to check for positive semi-definiteness (i.e their eigenvalues are non-negative.)

The comment above implies that one could do it by computing the respective eigenvalues and checking if they are non-negative (perhaps having to take care of rounding errors.)

Computing the eigenvalues is quite expensive in my scenario but I have noticed that the library that I am using has quite a fast test for positive definiteness (that is, if the eigenvalues of a matrix are strictly positive.)

Hence the idea would be, that given a matrix $B \in {\cal L}$, one tests if $B + \epsilon I$ is positive definite. If it is not then $B$ is not positive semi-definite, otherwise one can compute the eigenvalues of $B$ to make sure it is indeed positive semi-definite.

My question now is:

Is there a more direct and efficient way of testing whether a matrix is positive semi-definite, provided that an efficient test for positive definiteness is given?

• The test you noticed in the library is likely based on the proposition that symmetric real matrix $A$ is positive definite if and only if each leading principle minor gives a positive determinant, something that could be checked by elimination without pivoting in exact arithmetic. The subtle difficulty of extending this to a semi-definite case has lured many authors into misstatement. I know the topic has been broached in a Math.SE Question, so I'll try to provide a link.
– hardmath
Jun 24, 2014 at 15:00
• math.stackexchange.com/questions/753825/…
– hardmath
Jun 24, 2014 at 15:11
• For those more knowledgable people here - would it work to shift the spectrum to be positive by adding $B + cI$ for large $c$, then find the minimum eigenvalue of the shifted system (eg with inverse iteration), then check if the smallest eigenvalue of the shifted system is smaller than the shift $c$? The shift $c$ could be, for example, the largest in magnitude eigenvalue which can be found quickly. Jun 24, 2014 at 21:51
• Yes, you can shift the eigenvalues and compute the smallest eigenvalue, but you still have the problem of setting some tolerance for what you'll accept (and of making sure that your eigenvalues are computed to at least that tolerance!) Jun 25, 2014 at 1:53
• Not sure whether this would be helpful, but note that once you know a matrix is not positive definite, to check whether it is positive semidefinite you just need to check whether its kernel is non-empty. Jun 30, 2014 at 19:34

What's your working definition of "positive semidefinite" or "positive definite"? In floating point arithmetic, you'll have to specify some kind of tolerance for this.

You could define this in terms of the computed eigenvalues of the matrix. However, you should first notice that the computed eigenvalues of a matrix scale linearly with the matrix, so that for example, the matrix I get by multiplying $A$ by a factor of one million has its eigenvalues multiplied by a million. Is $\lambda=-1.0$ a negative eigenvalue? If all of the other eigenvalues of your matrix are positive and on the order of $10^{30}$, then $\lambda=-1.0$ is effectively 0 and shouldn't be treated as a negative eigenvalue. Thus it's important to take scaling into account.

A reasonable approach is to compute the eigenvalues of your matrix, and declare that the matrix is numerically positive semidefinite if all eigenvalues are larger than $-\epsilon \left| \lambda_{\max} \right|$, where $\lambda_{\max}$ is the largest eigenvalue.

Unfortunately, computing all of the eigenvalues of a matrix is rather time consuming. Another commonly used approach is that a symmetric matrix is considered to be positive definite if the matrix has a Cholesky factorization in floating point arithmetic. Computing the Cholesky factorization is an order of magnitude faster than computing the eigenvalues. You can extend this to positive semidefiniteness by adding a small multiple of the identity to the matrix. Again, there are scaling issues. One fast approach is to do a symmetric scaling of the matrix so that the diagonal elements are 1.0 and add $\epsilon$ to the diagonal before computing the Cholesky factorization.

You should be careful with this though, because there are some problems with the approach. For example, there are circumstances where the $A$ and $B$ are postive definite in the sense that they have floating point Cholesky factorizations, but $(A+B)/2$ does not have a Cholesky factorization. Thus the set of "floating point Cholesky factorizable positive definite matrices" isn't convex!

• Could you elaborate on that last paragraph or post a link to a source? That is pretty bizarre. Jun 24, 2014 at 16:56
• A classic reference to this scaling is A. van der Slui. Condition numbers and equilibration of matrices Numerische Mathematik 14(1):14-23, 1969. It's also discussed in textbooks such as Golub and van Loan. The bit in the last paragraph is from hard won personal experience in coding line search in a semidefinite programming code- I've encountered situations where $X$ and $X+\alpha \Delta X$ have Cholesky factorizations by LAPACK, but $X+0.95 \alpha \Delta X$ does not have a Cholesky factorization according to LAPACK. These kinds of problems start to occur when you're nearly singular. Jun 25, 2014 at 1:49
• It's also not uncommon to discover that some matrices can be Cholesky factored in extended or quadruple precision but not in regular double precision or single precision floating point arithmetic. Jun 25, 2014 at 1:50
• Several of the primal-dual interior point codes for SDP (CSDP, SDPT3, SDPA) always return matrices that are positive definite and have Cholesky factorizations, while another popular solver (SeDuMi) uses an eigenvalue decomposition and will return solutions that have very small negative eigenvalues. Jun 25, 2014 at 1:57
• Seems like the username pretty much discloses the relationship between the author of the answer and author of the papers. A little bit more info on what is contained in the paper would be nice; though, anyway, it is a very interesting and relevant to the question list of papers! Mar 17, 2018 at 3:25
• Anton Menshov: The three papers that I offered above have titles underlined that are, in fact, active links; so if the reader merely clicks on them, the abstract is revealed containing the summary discussion that you seek on what is contained in the papers. So I have provided both explicit citation to the papers and a contents summary (and pin-point sources that I referenced as a prelude and the sequencing of the papers in time reveals my motivation as well). Dec 10, 2022 at 3:37

Python for @Brian Borchers' suggestion, try a Cholesky decomposition:

from sksparse.cholmod import cholesky, CholmodNotPositiveDefiniteError

def testposdef( A, beta=1e-6, pr="" ):
""" try cholesky( a scipy.sparse matrix  + beta I )

Why A + beta I, beta say 1e-6 ?
If A has tiny eigenvalues, 0 to within machine precision,
about half of these "zeros" may be negative -- tough on solvers.
Also the condition number improves to ~ rho(A) / beta.
"""
try:
solve = cholesky( A, beta=beta )  # A + beta I
if pr:
print( "%s + %g I is positive-definite" % (pr, beta ))
return solve  # x = solve( b )
except CholmodNotPositiveDefiniteError:
if pr:
print( "%s + %g I is not positive-definite" % (pr, beta ))
return False