Testing if a matrix is positive semi-definite

Question

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?

Answer 1

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!

Answer 2

• Kerr, T. H., “Testing Matrices for Definiteness and Application Examples that Spawn the Need,” AIAA Journal of Guidance, Control, and Dynamics, Vol. 10, No. 5, pp. 503-506, Sept.-Oct., 1987.

• Kerr, T. H., “On Misstatements of the Test for Positive Semidefinite Matrices,” AIAA Journal of Guidance, Control, and Dynamics, Vol. 13, No. 3, pp. 571-572, May-Jun. 1990. (as occurred in Navigation & Target Tracking software in the 1970’s & 1980’s using counterexamples)

• Kerr, T. H., “Fallacies in Computational Testing of Matrix Positive Definiteness/Semidefiniteness,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 26, No. 2, pp. 415-421, Mar. 1990. [Lists fallacious algorithms that the author found to routinely exist in U.S. Navy submarine navigation and sonobuoy software in the late 1970’s and early 1980’s using counterexamples to point out the problems.]

Answer 3

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

from sksparse.cholmod import cholesky, CholmodNotPositiveDefiniteError
    # https://scikit-sparse.readthedocs.io/en/latest/cholmod.html

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