Eigen - Check if matrix is Positive (Semi-)Definite

Question

I'm implementing a spectral clustering algorithm and I have to ensure that a matrix (laplacian) is positive semi-definite.

A check if the matrix is positive definite (PD) is enough, since the "semi-" part can be seen in the eigenvalues. The matrix is pretty big (nxn where n is in the order of some thousands) so eigenanalysis is expensive.

Is there any check in Eigen that gives a bool result in runtime?

Matlab can give a result with the chol() method by throwing an exception if a matrix is not PD. Following this idea, Eigen returns a result without complaining for LLL.llt().matrixL(), although I was expecting some warning/error. Eigen also has the method isPositive, but due to a bug it is unusable for systems with an old Eigen version.

Solution 1

You can use a Cholesky decomposition (LLT), which returns Eigen::NumericalIssue if the matrix is negative, see the documentation.

Example below:

#include <Eigen/Dense>

#include <iostream>
#include <stdexcept>

int main()
{
    Eigen::MatrixXd A(2, 2);
    A << 1, 0 , 0, -1; // non semi-positive definitie matrix
    std::cout << "The matrix A is" << std::endl << A << std::endl;
    Eigen::LLT<Eigen::MatrixXd> lltOfA(A); // compute the Cholesky decomposition of A
    if(lltOfA.info() == Eigen::NumericalIssue)
    {
        throw std::runtime_error("Possibly non semi-positive definitie matrix!");
    }    
}

Solution 2

In addition to @vsoftco 's answer, we shall also check for matrix symmetry, since the definition of PD/PSD requires symmetric matrix.

Eigen::LLT<Eigen::MatrixXd> A_llt(A);
if (!A.isApprox(A.transpose()) || A_llt.info() == Eigen::NumericalIssue) {
    throw std::runtime_error("Possibly non semi-positive definitie matrix!");
}    

This check is important, e.g. some Eigen solvers (like LTDT) requires PSD(or NSD) matrix input. In fact, there exists non-symmetric and hence non-PSD matrix A that passes the A_llt.info() != Eigen::NumericalIssue test. Consider the following example (numbers taken from Jiuzhang Suanshu, Chapter 8, Problem 1):

Eigen::Matrix3d A;
Eigen::Vector3d b;
Eigen::Vector3d x;

// A is full rank and all its eigen values >= 0
// However A is not symmetric, thus not PSD
A << 3, 2, 1, 
     2, 3, 1, 
     1, 2, 3;
b << 39, 34, 26;

// This alone doesn't check matrix symmetry, so can't guarantee PSD
Eigen::LLT<Eigen::Matrix3d> A_llt(A);
std::cout << (A_llt.info() == Eigen::NumericalIssue) 
          << std::endl;  // false, no issue detected

// ldlt solver requires PSD, wrong answer
x = A.ldlt().solve(b);
std::cout << x << std::endl;  // Wrong solution [10.625, 1.5, 4.125]
std::cout << b.isApprox(A * x) << std::endl;  // false

// ColPivHouseholderQR doesn't assume PSD, right answer
x = A.colPivHouseholderQr().solve(b);
std::cout << x << std::endl;  // Correct solution [9.25, 4.25, 2.75]
std::cout << b.isApprox(A * x) << std::endl;  // true

Notes: to be more exact, one could apply the definition of PSD by checking A is symmetric and all of A's eigenvalues >= 0. But as mentioned in the question, this could be computationally expensive.

Comments

[Taylor]

Confusing because that documentation also says "Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices, use LDLT instead for the semidefinite case." I guess it depends on how you define "rank-revealing."