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.

Cannot you check if it's hermitian first, then look at the eigenvalues? Checking for hermiticity is straightforward. — vsoftco, Feb 05 '16 at 15:00
You are right about the hermitian part, but ideally I would like to avoid computing eigenvalues for a huge matrix several times, since this is my desired output so I would like this to happen just once if possible. — dim_tz, Feb 05 '16 at 15:05
Perhaps you can try Cholesky decomposition from Eigen, and that returns `NumericalIssue` if the matrix is negative, see http://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html — vsoftco, Feb 05 '16 at 15:12
Probably there is a versioning problem, because: `error: 'class Eigen::LDLT has no member named 'info'` — dim_tz, Feb 05 '16 at 17:24
Oh actually it works for LLT, although it doesn't for LDLT, thanks for the pointer! If you want you can write an answer so that I accept this, otherwise I will post a snippet as an answer later. — dim_tz, Feb 05 '16 at 17:37
Glad it helped, I wrote an answer. I'm not sure if it's the best way to go, especially since you are actually interested in the smallest eigenvalue being greater or equal to zero (and hermiticity of course), but it does the job. — vsoftco, Feb 05 '16 at 18:08

vsoftco · Accepted Answer · 2016-02-05T18:15:02.940

16

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!");
    }    
}

edited Feb 05 '16 at 18:15

answered Feb 05 '16 at 18:07

vsoftco

55,410
12
139
252

2

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." – Taylor Jun 13 '19 at 21:19

Yixing Lao · Answer 2 · 2019-02-07T10:24:24.513

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.

brice rebsamen · Answer 3 · 2022-04-11T12:42:04.043

you have to test that the matrix is symmetric (A.isApprox(A.transpose())), then create the LDLT (and not LLT because LDLT takes care of the case where one of the eigenvalues is 0, ie not strictly positive), then test for numerical issues and positiveness:

template <class MatrixT>
bool isPsd(const MatrixT& A) {
  if (!A.isApprox(A.transpose())) {
    return false;
  }
  const auto ldlt = A.template selfadjointView<Eigen::Upper>().ldlt();
  if (ldlt.info() == Eigen::NumericalIssue || !ldlt.isPositive()) {
    return false;
  }
  return true;
}

I tested this on

1 2
2 3

which has a negative eigenvalue (hence not PSD). Without the isPositive() test, isPsd() incorrectly returns true here.

and on

1 2
2 4

which has a null eigenvalue (hence PSD but not PD).

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

3 Answers3