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
Summary: Fail to compute eigenvalues for a simple 3x3 companion matrix for root finding
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Eigenvalues (show other bugs)
Version: 3.3 (current stable)
Hardware: All All
: Normal Wrong Result
Assignee: Nobody
URL:
Whiteboard:
Keywords: test-needed
: 1174 (view as bug list)
Depends on:
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2018-06-15 07:16 UTC by Gary Tan
Modified: 2018-12-12 17:05 UTC (History)
6 users (show)



Attachments

Description Gary Tan 2018-06-15 07:16:59 UTC
The following

#include <Eigen/Eigenvalues>
#include <iostream>
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<MatrixXd> 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<double> lambda = es.eigenvalues()[0];
  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<complex<double> >() * 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
Comment 1 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.
Comment 2 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;`
Comment 3 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).
Comment 4 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.
Comment 5 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)
Comment 6 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.