New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 1557 - Fail to compute eigenvalues for a simple 3x3 companion matrix for root finding
Fail to compute eigenvalues for a simple 3x3 companion matrix for root finding
 Status: RESOLVED FIXED None Eigen Unclassified Eigenvalues (show other bugs) 3.3 (current stable) All All Normal Wrong Result Nobody test-needed (view as bug list) 3.4 Show dependency tree / graph

 Reported: 2018-06-15 07:16 UTC by Gary Tan 2018-12-12 17:05 UTC (History) 6 users (show) chtz gael.guennebaud jacob.benoit.1 jitseniesen pcv22 tgntanguangning

Attachments

 Gary Tan 2018-06-15 07:16:59 UTC ```The following #include #include using namespace Eigen; using namespace std; int main() { Matrix3d A; A << 0, 0, 0, 1, 0, 0.5887907064808635127, 0, 1, 0; cout << "Here is a random 3x3 matrix, A:" << endl << A << endl << endl; EigenSolver es(A); cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl; cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl; complex lambda = es.eigenvalues(); cout << "Consider the first eigenvalue, lambda = " << lambda << endl; VectorXcd v = es.eigenvectors().col(0); cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl; cout << "... and A * v = " << endl << A.cast >() * v << endl << endl; MatrixXcd D = es.eigenvalues().asDiagonal(); MatrixXcd V = es.eigenvectors(); cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl; return 0; } produces Here is a random 3x3 matrix, A: 0 0 0 1 0 0.588791 0 1 0 The eigenvalues of A are: (0,0) (0,0) (0,0) The matrix of eigenvectors, V, is: (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) Consider the first eigenvalue, lambda = (0,0) If v is the corresponding eigenvector, then lambda * v = (0,0) (0,0) (0,0) ... and A * v = (0,0) (1,0) (0,0) Finally, V * D * V^(-1) = (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) (-nan,-nan) However, it's expected to obtain sqrt(0.5887907064808635127) = 0.76733 and its opposite. Thanks, Gary``` Gary Tan 2018-06-15 07:58:21 UTC ```Matrix3d A; A << 0, 0, 0, 1, 0, 0, 0, 1, 0; double a = double(0.5887907064808635127); A(1,2) = a; This does not work either. If the number "0.5887907..." were "shorter", then it'd work fine. However, in general that number is assigned from other variable, say a. I tried to debug, but could not figure out why. Sorry for creating this weird problem.``` Christoph Hertzberg 2018-06-15 08:06:12 UTC ```If you check `es.info()` you will notice that it did not converge -- thus the computed values will not be useful. It is strange that it does not converge, however. It may be failing, because all Eigenvalues are indeed real, e.g., you get a valid result if you set `A(0,2)=1e-30;```` Gary Tan 2018-06-15 08:19:40 UTC ```(In reply to Christoph Hertzberg from comment #2) > If you check `es.info()` you will notice that it did not converge -- thus > the computed values will not be useful. > > It is strange that it does not converge, however. It may be failing, because > all Eigenvalues are indeed real, e.g., you get a valid result if you set > `A(0,2)=1e-30;` Hi, Chris! If I use 0.588790706480863 or 0.588790706480864 then it computes fine and produces expected values. This strange problem starts with 0.5887907064808635 which is just a few epsilons away. It seems internally something wrongly handles this number. I ran the same code on macOS and got the same results. I am using g++ (Debian 6.3.0) and Apple LLVM v10.0.0 (clang-1000.10.25.5 on macOS 10.14).``` Christoph Hertzberg 2018-12-04 11:32:55 UTC ```Bisecting leads to this commit breaking your example: https://bitbucket.org/eigen/eigen/commits/c90f5b999226 There is hopefully a way to fix this without breaking Bug 933.``` Gael Guennebaud 2018-12-12 16:35:51 UTC ```Thank you both for the findings. Fixed with regression tests for both 933 and this one: https://bitbucket.org/eigen/eigen/commits/4d387242247e/ https://bitbucket.org/eigen/eigen/commits/7e7fef7bac37/ (3.3)``` Gael Guennebaud 2018-12-12 17:05:07 UTC `*** Bug 1174 has been marked as a duplicate of this bug. ***`

 Note You need to log in before you can comment on or make changes to this bug.