Method to check for positive definite matrices

Question

I think it's already been asked, but I still can't figure out a way to do it computationally.

I had to check for positive definiteness of an $n \times n$ matrix $A$. I know that for any nonzero column vector $x$, we have $x^{t}A x > 0$.

I want to implement it in C++ to use it with the Cholesky method, in which I need:

1. symmetric matrix; and
2. positive definite matrix.

How do I approach this situation? If I go for the definition then it can be computationally inefficient as I have to check for each vector $x$. Is there any other method that I can use to test for positive definiteness?

Answer 1

The standard way of approaching this is really to attempt a Cholesky factorization and check to see if such a factorization exists. This is both fast and realiable. Here it is in MATLAB notation:

A = zeros(3); % some matrix

[~,p] = chol(A)

If the input matrix is not positive definite, then $p$ will be a positive integer, e.g. $p=1$ and MATLAB will not generate an error. $p$ could also be positive for some positive semi-definite matrices. Of course 'chol' only considers the upper (or lower) triangular part of the matrix.

Another way is to check the following condition:

p == 0 && rank(A) == size(A,1)

this involves @Bill Greene's suggestion. Likewise, you might use the Eigen decomposition:

all(eig(A) > eps)

That's it to the best of my knowledge. Note that MATLAB uses a Lapack library internally (such as MKL) and you should be able to get similar functionality out of any Lapack library in C++.

Answer 2

Mostly, I'm leaving this answer here as a cautionary tale to not use a Choleski factorization. Most of the time, this is a fine answer. However, very specifically, it can and will fail, so if this is in a process where people can be hurt, please do not use it.

First, a Choleski factorization exists for matrices that have zero eigenvalues. This is a property that exists without round-off error. To see this, we can look at a $2\times 2$ matrix. Here, we have \begin{align*} L_{11} =& \sqrt{A_{11}}\\ L_{21} =& \frac{A_{21}}{\sqrt{A_{11}}}\\ L_{22} =& \sqrt{A_{22}-A_{21}^2/A_{11}} \end{align*} Then, consider the matrix $$ \begin{bmatrix}2 & 2\\2 & 2\end{bmatrix} $$ The Choleski factorization is $$ \begin{bmatrix}\sqrt{2} & 0\\\sqrt{2} & 0\end{bmatrix} $$ Clearly, this matrix is rank-1 and has a zero eigenvalue. Further, MATLAB is perfectly willing to compute it

>> A=ones(2)*2; chol(A)    

ans =

    1.4142    1.4142
         0    0.0000

Long story short, Choleski works fine both theoretically and practically for some matrices at have zero eigenvalues.

Next, the question becomes: With round off error, can we Choleski factor a matrix with negative eigenvalues? The answer is also yes. Here's a code to find one for you:

for i=1:100
    try
        A11=rand*1e5;
        A22=rand;
        A21=sqrt(A11*A22)+1e-14;
        A=[A11 A21;A21 A22];
        U=chol(A);
        [v d]=eig(A);
        if d(1) < 0
            break;
        end 
    catch myerr
        i
    end
end

I just ran it and found:

   7.89575133795518e+04   2.47689038064322e+02
   2.47689038064322e+02   7.76998374838857e-01

Note, the matrix is design precisely to be indefinite with a negative eigenvalue. With round off, sometimes even eig gives the wrong answer.

Where does this leave us? Unfortunately, not in a good spot. Most of the time, trying a Choleski factorization is the way to go. It's fast and is mostly reliable. However, as I said at the top, this can and will fail. I had no idea until it broke one of my codes. Now, I have safety margins. Is there a more reliable way? Sometimes. If you use something like an implicitly restarted Arnoldi method, such as with ARPACK, and you ask for the smallest eigenvalue with an error tolerance on the computation. If the smallest eigenvalue minus this error tolerance is greater than zero, then you're good to go. If not, there can be trouble, which requires additional investigation. Beyond that, I don't know of any perfectly reliable way to check.

However, let me state one more time because it's important:

With round off, a Choleski factorization can be taken of a matrix with negative eigenvalues.

---

Edit 1

There was a follow-up question on how to accomplish this in C++. Unfortunately, I don't know a fully straightforward way, but I can offer a few suggestions. Note, each of these methods involves calculating the smallest eigenvalue of the symmetric matrix to a specified tolerance. If the smallest eigenvalue minus the tolerance is greater than or equal to zero, then we know we're positive definite.

Frankly, your best bet is to use ARPACK. This is a Fortran code, but if you compile it with something like gfortran and write your own headers, it will work just fine. Using Fortran is a little tricky in C++. You'll have to turn off C++ name mangling with extern "C" {} and determine the Fortran name mangling scheme which can be FOO, FOO_, foo, or foo_. To determine what your compiler is doing, use the utility nm on Linux/Unix. Alternatively, both autotools and CMake have routines to determine the Fortran name mangling scheme. Anyway, ARPACK is good, really good. Just ask it for the smallest algebraic eigenvalue and specify some kind of tolerance.

Second, use LAPACK and the routine syevr. Make sure to specify that you only want one eigenvalue, set RANGE to I and set IL and IU. Like ARPACK, LAPACK is in Fortran, so make sure to correctly deal with name mangling. There will be fewer routines that you need to write a header for, namely only one whereas ARPACK will require more, so it's simpler to use. That said, for this case, ARPACK will work better.

Third, the paper that ARPACK is based on is called "Implicit Application of Polynomial Filters in a k-Step Arnoldi Method" by Dan Sorensen (link). It's a good paper and worth a read. That said, really what you're doing is using the Arnoldi method to reduce your matrix to a tridiagonal form and then looking at the eigenvalues of this form. The Arnoldi method really is just Gram-Schmidt applied to the Krylov subspace. Basically, we get a relationship where $AV = VH + \textrm{stuff}$. The structure of $H$ is $$ H=\begin{bmatrix} \langle \cdot,\cdot \rangle & \langle \cdot,\cdot \rangle & \dots & \langle \cdot,\cdot \rangle\\ \|\cdot\| & \langle \cdot,\cdot \rangle & \dots & \langle \cdot,\cdot \rangle\\ 0 & \|\cdot\| & \dots & \langle \cdot,\cdot \rangle\\ \vdots & 0 & \ddots & \vdots\\ 0 & 0 & \dots & \langle \cdot,\cdot \rangle \end{bmatrix} $$ Here, the $\langle \cdot,\cdot \rangle$ are the Gram-Schmidt coefficients and the $\|\cdot\|$ are the normalization of the Krylov vectors after normalization. For a symmetric matrix, $H$ should be tridiagonal, but it really never will be due to round-off error. Anyway, the eigenvalues of $H$ are called Ritz values and they approximate the eigenvalues of $A$. If I recall correctly, the Ritz values will approach the eigenvalues of $A$ from the inside and find the outermost eigenvalues first. How does this help you? You can just write your own eigenvalue routine. Basically, use the method from Sorensen's paper, but you can really skip the implicit QR iterations if you only want one eigenvalue. Literally, what you'll be doing is coding the Arnoldi method and then finding the eigenvalues of $H$. Restarting will be helpful and the paper explains that, but, again, you can probably skip the implicit QR iterations. Once the outermost eigenvalues converge to a certain tolerance, pick the smallest one and that's your smallest eigenvalue to the tolerance.

Anyway, yeah, this is kind of a pain and I'm sorry about that. Really, ARPACK is the right tool even though there will be some pain using it.

Answer 3

If your idea is to use this info for the Cholesky decomposition I assume that you are working with symmetric matrix.

In the case of $A \in M(\mathbb{R}^n)$ and a, i.e. a real matrix, we can add to @Tolga_Birdal's answer the consideration that all the eigenvalues are real. With this you can check only the smallest eigenvalue, if it is $\lambda_{min} > 0$, not all.