CODEWITHSUNDEEP
c++eigeneigen3
dim_tz
dim_tz

Reputation: 1531

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

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.

Upvotes: 11

Views: 14900

Answers (3)

brice rebsamen
brice rebsamen

Reputation: 714

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

Upvotes: 2

Yixing Lao
Yixing Lao

Reputation: 1428

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.

Upvotes: 10

vsoftco
vsoftco

Reputation: 56567

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

Upvotes: 19

Related Questions