Hi all,

I am getting a bit of contradicting results so I was wondering if someone can shed some light.
I have a symmetric matrix, that should be positive definite.

I usually use L*L^T to solve a normal equations matrix. In this specific case I am looking at, I get no errors.
I checked using L*D*L^T function isPositive (I understand that the check is for semi-definite) and it returns true.
I check the eigenvalues using SelfAdjointEigenSolver, there is 1 negative eigenvalue (-1.4384e+10)
I check the eigenvalues using matlab, all are > 0.
I check the eigenvalues using EigenSolver all are > 0.

What am I missing here ?

Which Eigen version are you using? Are you using float or double? complex or real?

Could you share the matrix or a self-contained exemple so that we can reproduce the issue?

Hi ggael,

I am using 3.2.9 but also tried 3.3-beta2.
The values are real doubles.

There is no problem in sharing the matrix so here is a reproduction sample and the matrix in comma delimited format.

Reproduction code
https://app.box.com/s/pbpcff3baz80abafje7k6kmy1byk3sly
Matrix
https://app.box.com/s/7isc1kh1kcp3rmsbr4jfmuy8pcyyisl4

Thank you, I can reproduce. There seems to be some overflow or the likes as scaling the matrix with:

mat = mat/mat.norm();

fixes the inconsistency. But this is weird as I'm pretty sure that SelfAdjointEigenSolver performs such a scaling under the hood... I'll investigate...

ah, wrong alarm, actually this -1.4384e+10 eigenvalue is numerically equal to 0 because the largest eigenvalue is much larger. Your test should be:

bool adj_es_ispos = !((adj_es.eigenvalues().array() < -std::numeric_limits<double>::epsilon()*adj_es.eigenvalues()(mat.cols()-1)).any());

For the record, here are the eigenvalues computed by SelfAdjointEigenSolver:

Code: Select all

-5.16561e+08    0.0768628      52.7249      102.535      226.034       425.64      471.342      494.873      540.352      561.638 ... 6.88748e+11 7.04878e+11 7.15946e+11 7.24494e+11  7.2881e+11 3.80591e+14 4.19592e+14  5.7554e+14 3.40776e+19 3.95004e+25


and the one computed by the equivalent routine in Lapack:

Code: Select all

-3.21481e+10 -6.72629e+09      1.01093      97.1779      106.476      248.211      426.369      471.102      494.625       544.21 ... 6.88747e+11 7.04948e+11 7.15957e+11  7.2453e+11 7.28839e+11 3.80591e+14 4.19592e+14  5.7554e+14 3.40776e+19 3.95004e+25


Thanks ggael,

That makes sense,
Below are the eigen values that I get from matlab (R2012a) which are all positive

Code: Select all

0.110583862411932
55.3799486706705
102.701133941637
229.403566917025
413.229196031421
471.339657806331
494.325394238103
549.302100205325
560.944047280822
....
724500242927.785
728810656908.836
380590611806053
419591813628402
575539720666378
3.40776085157530e+19
3.95004019281995e+25

Yes, matlab calls a generic eigen solver, more like Eigen::EigenSolver instead of a one specialized for symmetric problems. Without more precision there is not much we can do in SelfAdjointEigenSolver because with double, there is no way to assess the accuracy of eiegnvalues smaller than 1e-16 * largest eigenvalues. Actually, if you call it with "long double", SelfAdjointEigenSolver also returns positive eigenvalues.