Eigen - Check if matrix is Positive (Semi-)Definite
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!");
}
}
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.
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).