Geek Answers Handbook

Eigen - check if the matrix is ​​positive (semi)

I am implementing a spectral clustering algorithm, and I must ensure that the matrix (Laplacian) is positively semi-definite.

We check if the matrix (PD) is sufficiently well defined, since the “half-" part can be seen in eigenvalues. The matrix is ​​quite large (nxn, where n is on the order of thousands), so our own analysis is expensive.

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

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

+8
source share
2 answers

You can use 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!"); } }
+11
source

In addition to @vsoftco's answer, we will also check the symmetry of the matrix, since a symmetric matrix is ​​required to determine PD / PSD.

 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, for example, for some proprietary solvers (such as LTDT ) a matrix input PSD (or NSD) is required. In fact, there is an asymmetric matrix A which passes A_llt.info() != Eigen::NumericalIssue . Consider the following example (figures taken from Jiuzhan 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: more precisely, the definition of PSD can be applied by checking that A symmetric and all eigenvalues ​​A> = 0. But, as mentioned in this question, this can be computationally expensive.

+1
source

More articles:


All Articles