This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
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
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2017-11-20 19:36 UTC by Felix Schindler
Modified: 2019-12-04 17:17 UTC (History)
4 users (show)



Attachments
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]
MWE

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();
Comment 2 Nobody 2019-12-04 17:17:51 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1488.

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