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 '14 at 15:00
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. –  Nick Alger Jun 24 '14 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!) –  Brian Borchers Jun 25 '14 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. –  Abel Molina Jun 30 '14 at 19:34

1 Answer

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$ is 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_{\mbox{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 diagional 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. –  Daniel Shapero Jun 24 '14 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. –  Brian Borchers Jun 25 '14 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. –  Brian Borchers Jun 25 '14 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. –  Brian Borchers Jun 25 '14 at 1:57