6
$\begingroup$

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?

$\endgroup$
  • 7
    $\begingroup$ I assume you would like to check for a positive definite matrix before attempting a Cholesky decomposition? As far as I know, this is not possible. You simply have to attempt a Cholesky factorization and abandon it if you encounter a zero or negative pivot. This is the approach the MATLAB backslash operator takes for square, symmetric matrices. $\endgroup$ – Bill Greene Feb 16 '17 at 15:55
  • $\begingroup$ Ok,if as a new question, i were to check a matrix is positive definite , then i need to check for positive definite and i am searching a way to code it efficiently! $\endgroup$ – BAYMAX Feb 16 '17 at 16:01
  • 2
    $\begingroup$ You can't check each vector $x$, since they are infinite. You can check that your matrix is symmetric, and then, theoretically, compute all eigenvalues or the determinant of all the minors of the matrix. $\endgroup$ – nicoguaro Feb 16 '17 at 16:16
  • $\begingroup$ Yes, thanks, you are right of course. I deleted my comment about SVD. $\endgroup$ – Bill Greene Feb 16 '17 at 18:46
  • $\begingroup$ Guys , i want to do it in C++ , in MATLAB due to various commands i can do it but in C++ how i can proceed ? $\endgroup$ – BAYMAX Mar 1 '17 at 13:42
9
$\begingroup$

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++.

$\endgroup$
  • 1
    $\begingroup$ Indeed, a matrix is positive definite iff it has a Cholesky decomposition. This is what you should use in practice (on a computer) to check for positive definiteness. $\endgroup$ – GoHokies Feb 16 '17 at 22:11
  • $\begingroup$ Yes the first two yellow lines do exactly that. $\endgroup$ – Tolga Birdal Feb 16 '17 at 22:21
  • 2
    $\begingroup$ This is not entirely correct. First, the Choleski factorization is actually defined on indefinite matrices. For example, in MATLAB, A=ones(2)*2;chol(A) computes even though A has a zero eigenvalue. Note, this is not a round off error issue as the factorization exists for some matrices with zero eigenvalues. However, with round off, this can become an issue. I'm having trouble generating a cute matrix at the moment, but sometimes with roundoff we can Choleski a matrix with small negative eigenvalues. Most of the time, Choleski is the way to go. If lives are at stake, do not use it. $\endgroup$ – wyer33 Feb 20 '17 at 5:13
  • 1
    $\begingroup$ There are many Lapack solutions (Armadillo, Netlib, Flens...). Maybe you could already and easily start with Eigen (eigen.tuxfamily.org/dox/classEigen_1_1LDLT.html#details). If that doesn't work, you could move to rather raw Lapack solutions. $\endgroup$ – Tolga Birdal Mar 1 '17 at 15:52
  • 1
    $\begingroup$ @BAYMAX, try to use Cholesky from Eigen and the method that I explained. If Cholesky exists, you are (most probably) good to go. $\endgroup$ – Tolga Birdal Mar 2 '17 at 19:52
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Please give some insight on how to do in C++... $\endgroup$ – BAYMAX Mar 1 '17 at 13:43
  • $\begingroup$ @BAYMAX It turns out it's sort of a pain, but I added some info. $\endgroup$ – wyer33 Mar 2 '17 at 16:30
  • $\begingroup$ Thank you for the information , it will take some time to digest , i will do it .. thanks $\endgroup$ – BAYMAX Mar 2 '17 at 16:45
  • $\begingroup$ Note that it is not necessarily the case that that matrix you passed to chol has a negative eigenvalue. In addition to the errors in eig that you mentioned, it is possible that the computational error in sqrt or in the multiplication makes it positive definite even in exact arithmetic. $\endgroup$ – Federico Poloni Mar 2 '17 at 20:09
  • 3
    $\begingroup$ Isn't this just describing a general failure of a backward-stable algorithm (Cholesky factorization is backward-stable) on an ill-conditioned problem (determining whether $\lambda\geq0$ is ill-conditioned—all comparisons are ill-conditioned at the boundary)? For your matrix, whether Cholesky factorization succeeds/fails can be changed by adding a minute random perturbation of size $O(\epsilon_{\mathrm{mach}})$. $\endgroup$ – Kirill Mar 2 '17 at 23:13
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ and what would be a good (that is, efficient) method to find this smallest (potentially negative) eigenvalue? $\endgroup$ – GoHokies Feb 16 '17 at 21:38
  • $\begingroup$ mathoverflow.net/questions/24287/… $\endgroup$ – Tolga Birdal Feb 16 '17 at 21:50
  • $\begingroup$ @Tolga Birdal all of those methods compare very unfavorably to the Cholesky decomposition in terms of computational efficiency and/or numerical stability. $\endgroup$ – GoHokies Feb 16 '17 at 21:57
  • $\begingroup$ I agree. chol is my preference - look at my answer. But that doesn't mean the alternatives don't exist. $\endgroup$ – Tolga Birdal Feb 16 '17 at 22:20
  • $\begingroup$ Well the idea of chol is good, and correct because as @GoHokies noted it is a characterization of SPD matrix. The check of eigenvalues, or the smallest eigenvalue, can be an alternative, i.e. the idea express in TolgaBirtal use the function of matlab chol that it returns the value $p$. Now I am not sure that this value $p$ is returned in a generic function of a different language/library, similar for the sparse case when you want the incomplete Cholesky (not present in the question). So the eigenvalue check is more agnostic respect a particular function used. $\endgroup$ – Mauro Vanzetto Feb 17 '17 at 9:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.