New user self-registration is disabled due to spam. Please email eigen-core-team @ if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 1488 - eigensolver gives complex eigenvectors for hyperbolic matrix
Summary: eigensolver gives complex eigenvectors for hyperbolic matrix
Status: NEW
Alias: None
Product: Eigen
Classification: Unclassified
Component: Eigenvalues (show other bugs)
Version: 3.3 (current stable)
Hardware: x86 - 64-bit Linux
: Normal Wrong Result
Assignee: Nobody
Depends on:
Reported: 2017-11-20 19:36 UTC by Felix Schindler
Modified: 2017-11-21 22:08 UTC (History)
4 users (show)

MWE (4.25 KB, text/x-c++src)
2017-11-20 19:36 UTC, Felix Schindler
no flags Details

Description Felix Schindler 2017-11-20 19:36:59 UTC
Created attachment 805 [details]

Dear all,

I am not sure wether this is a bug or me using the library in a wrong way.

I have a real 4x4 matrix arising in a Finite Volume discretization of the Euler equations in 2 space dimensions.
We know from theory that the matrix has only real eigenvalues and eigenvectors.
Using the `Eigen::EigenSolver`, I get (nearly) real eigenvalues, but complex eigenvectors.
Using other libraries (such as lapack), I get real eigenvalues and eigenvectors.

I am running Arch Linux (gcc 7.2.0) and eigen 3.3.4. The output of the attached program is:
This is the matrix in question:
           0            0            1            0
-1.13722e-20 -4.14193e-18  -0.00274562            0
 1.50769e-06   0.00109825  -6.6271e-18          0.4
 5.78686e-18 -4.54887e-21      1.39714 -5.79871e-18

These are (up to suitable reordering) analytically known eigenvalues:
-4.14193e-18 -4.14193e-18     0.747566    -0.747566

This is (one choice) of analytically known (right) eigenvectors:
           1            0      2.66116      2.66116
 -0.00274562      3.97879  -0.00730654  -0.00730654
-4.14193e-18            0      1.98939     -1.98939
 3.76922e-06   -0.0109242      3.71802      3.71802

And indeed, the max error of the decomposition is: 4.44089e-16

These are eigenvalues computed by lapack:
    0.747566    -0.747566 -4.38053e-18 -4.13868e-18

These are eigenvectors computed by lapack:
   -0.533695   0.00146532    -0.398972    -0.745647
    0.533695  -0.00146532    -0.398972     0.745647
    0.999927    0.0120743 -4.40182e-18 -3.69203e-05
          -1  0.000164793  4.13837e-18  3.31676e-06

And indeed, the max error of the decomposition is: 1.35669e-13

These are eigenvalues computed by eigen:
               (0.747566,0)               (-0.747566,0)  (-4.21066e-18,1.01069e-19) (-4.21066e-18,-1.01069e-19)

These are eigenvectors computed by eigen:
               (0.533695,0)               (-0.533695,0)         (-0.82212,0.569311)        (-0.82212,-0.569311)
            (-0.00146532,0)              (0.00146532,0)    (1.83753e-05,0.00215749)   (1.83753e-05,-0.00215749)
               (0.398972,0)                (0.398972,0)  (-6.00757e-17,2.98675e-17) (-6.00757e-17,-2.98675e-17)
               (0.745647,0)               (-0.745647,0)   (3.04829e-06,-8.0695e-06)    (3.04829e-06,8.0695e-06)

And indeed, the max error of the decomposition is: 4.44089e-16
But this decomposition contains eigenvectors with nontrivial imaginary part!
Comment 1 Gael Guennebaud 2017-11-21 22:08:34 UTC
This is indeed a tricky problem. Both Eigen and Lapack follows the same steps: reduction to Hessenberg, then Schur factorization, and finally eigenvector extraction. In order to figure out how Lapack managed to keep real eigenvalues, I asked Lapack for the hessenberg reduction and then asked Lapack for the eigenvalues, and I got complex ones. Same if I asked for the schur decomposition then for the eigenvalues. So looking at dgees implementation, I saw that it starts by permuting rows/columns to make it as close as possible to upper triangular form before starting with the hessenberg/schur stuffs. And indeed, if I first permute the first and last rows/columns (the eigenvalues/vectors are not changed), then Eigen also provides real eigenvalues/vectors:

  Eigen::PermutationMatrix<Eigen::Dynamic> P(4);
  P.indices() << 3, 1, 2, 0;
  matrix = (P * matrix) * P.inverse();

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