New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 645 - GeneralizedEigenSolver: missing computation of eigenvectors
GeneralizedEigenSolver: missing computation of eigenvectors
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Eigenvalues
3.2
All All
: Normal Feature Request
Assigned To: Nobody
: JuniorJob
: 1269 (view as bug list)
Depends on:
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2013-08-10 23:35 UTC by Gael Guennebaud
Modified: 2016-08-23 19:22 UTC (History)
7 users (show)



Attachments
Preliminary Patch (6.41 KB, patch)
2016-05-26 15:33 UTC, Tobias Wood
no flags Details | Diff
Patch 2 (7.86 KB, patch)
2016-05-26 21:35 UTC, Tobias Wood
no flags Details | Diff
Patch 3 (14.33 KB, patch)
2016-08-18 15:15 UTC, Tobias Wood
no flags Details | Diff

Description Gael Guennebaud 2013-08-10 23:35:25 UTC
GeneralizedEigenSolver is still incomplete: the computation of the eigenvectors is missing.
Comment 1 danny 2014-08-01 19:23:47 UTC
Hi - I was curious if there are any current plans to address this bug?
Comment 2 Gael Guennebaud 2014-08-02 01:56:36 UTC
AFAIK, nobody started working on it, so the job is open ;)
Comment 3 Tobias Wood 2016-05-24 15:07:54 UTC
Hello,
I also need this feature. I did try to work on it, using dtgevc as a template as mentioned here (http://stackoverflow.com/questions/29762623/generalizedeigensolver-not-giving-eigenvectors-in-eigen-library), but unfortunately my Fortran isn't strong enough to follow the implementation details, and I can't decipher the very brief algorithm in the "Further Details" section. Does anyone know where this algorithm comes from? If I can get another reference to follow I will have another try.
Comment 4 Tobias Wood 2016-05-26 13:35:11 UTC
Hi again,
I think I have managed to decipher the algorithm and have something that seems to work for my particular problem for real eigenvalues.

I am happy to try and extend this to complex as well, however I looked at the tests and currently they only calculate eigenvalues of random matrices. I think it might be more useful if I implemented the test matrix from here: https://www.nag.co.uk/numeric/fl/manual/pdf/F08/f08ykf.pdf so we have known eigenvalues/vectors to compare against?
Comment 5 Christoph Hertzberg 2016-05-26 14:06:43 UTC
Hi, thanks for working on this!

Regarding unit tests: You should generally avoid hard-coded unit tests. 
You can often produce a randomized test example with known results if you start by generating (randomized) results, then reconstruct the matrices from that and compare the result to the original (minding possible ambiguities).
Comment 6 Tobias Wood 2016-05-26 15:33:24 UTC
Created attachment 710 [details]
Preliminary Patch
Comment 7 Tobias Wood 2016-05-26 15:40:11 UTC
Comment on attachment 710 [details]
Preliminary Patch

Okay, I will try that. I've attached a preliminary patch with what I have done so far.

Unfortunately, this is failing to build for fixed size matrices. If I comment out line 56 of eigensolver_generalized_real.cpp, then the remaining tests build and pass fine (i.e. my changes do not interfere with the existing eigenvalue code). The build errors seem are about missing types and template parameters (ReturnType, CoeffReturnType, SizeAtCompileTime etc.) and stem from line 340 in my code:

m_eivec.col(i).real() = (mZT * v).normalized();

I am compiling from a clone of the repo, with a fresh hg pull; hg update. Any hints on whether this is something I'm doing wrong or something to do with the 3.3 changes is much appreciated!
Comment 8 Tobias Wood 2016-05-26 21:35:54 UTC
Created attachment 711 [details]
Patch 2

Apologies for spamming comments.

I fixed the problem in the previous patch, which was because I had something declared as VectorXd instead of VectorType which was making it fail when the specified type was Matrix4f.

I couldn't see a way to reconstruct both the A & B matrix from the eigenvalues/eigenvectors, so I opted to simply check that the values/vectors do actually solve the generalised eigenproblem. I think all 4 tests are passing.

Finally, I updated my code to use Eigen expressions instead of standard loops, which I hope is faster.

I will work on complex eigenvectors when I get a chance.
Comment 9 Gael Guennebaud 2016-05-27 12:48:34 UTC
Comment on attachment 711 [details]
Patch 2


Great start!

Some clean-up suggestions:

remove:
>+#include <iostream>


Should be preserved:
>@@ -127,10 +128,7 @@
>         m_alphas(size),
>         m_betas(size),
>         m_isInitialized(false),
>-        m_eigenvectorsOk(false),
>-        m_realQZ(size),
>-        m_matS(size, size),
>-        m_tmp(size)
>+        m_realQZ(size)
>     {}
> 
>     /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair.

Should be preserved:
>@@ -150,10 +148,7 @@
>         m_alphas(A.cols()),
>         m_betas(A.cols()),
>         m_isInitialized(false),
>-        m_eigenvectorsOk(false),
>-        m_realQZ(A.cols()),
>-        m_matS(A.rows(), A.cols()),
>-        m_tmp(A.cols())
>+        m_realQZ(A.cols())
>     {
>       compute(A, B, computeEigenvectors);
>     }


This assert should be preserved:
>-//  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");


>@@ -307,41 +303,56 @@
>+        if (computeEigenvectors) {

v should be declared outside the loop to avoid multiple malloc/free
>+          VectorType v = VectorType::Zero(A.cols());
>+          v.coeffRef(i) = 1.0;

Please remove the else!
>+          if (internal::isApprox(fabs(m_betas.coeffRef(i)), 0.0)) {
>+            // Singular eigenvalue, do nothing more
>+          } else {
>+            const Scalar a = real(m_alphas.coeffRef(i));
>+            const Scalar b = m_betas.coeffRef(i);
>+            for (Index j = i-1; j >= 0; j--) {
>+              const Index st = j+1;
>+              const Index sz = i-j; 

A.transpose().cwiseProduct(B).sum() -> A.dot(B) 
>+              v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(b*mS.block(j,st,1,sz) - a*mT.block(j,st,1,sz)).sum() / (b*mS.coeffRef(j, j) - a*mT.coeffRef(j,j));


Should be preserved:
>-  m_eigenvectorsOk = computeEigenvectors;


>diff -r 53d9a496f686 test/eigensolver_generalized_real.cpp
>--- a/test/eigensolver_generalized_real.cpp	Thu May 26 16:16:41 2016 +0200
>+++ b/test/eigensolver_generalized_real.cpp	Thu May 26 22:32:21 2016 +0100

mat.array().rowwise() * vec.transpose().array() -> mat * vec.asDiagonal()
>+  MatrixType RHS = (spdB * realEigenvectors).array().rowwise() * realEigenvalues.transpose().array();
>+  VERIFY_IS_APPROX(LHS, RHS);
> }
> 
> void test_eigensolver_generalized_real()
Comment 10 Gael Guennebaud 2016-05-27 12:50:39 UTC
Also, handling non SPD matrices is a requirement. For those we already have the GeneralizedSelfAdjointEigenSolver class which is much faster (perhaps less robust though).
Comment 11 Tobias Wood 2016-05-31 12:02:43 UTC
Thanks for the feedback Gael.

I started to work on complex eigenvectors however I have hit a major problem, because the current code appears to be giving me incorrect complex eigenvalues. I suspect this is because RealQZ isn't giving a correct decomposition.

Taking for example the matrix from Eigen's page on RealQZ (https://eigen.tuxfamily.org/dox/classEigen_1_1RealQZ.html), the following Eigen code:

  Matrix4f a; //= MatrixType::Random(rows,cols);
  a <<    0.68,   0.823,  -0.444,   -0.27,
        -0.211,  -0.605,   0.108,  0.0268,
         0.566,   -0.33, -0.0452,   0.904,
         0.597,   0.536,   0.258 ,  0.832;
  Matrix4f b; //= MatrixType::Random(rows,cols);
  b << - 0.271, -0.967, -0.687,  0.998,
         0.435, -0.514, -0.198, -0.563,
        -0.717, -0.726,  -0.74, 0.0259,
         0.214,  0.608, -0.782,  0.678;
  RealQZ<Matrix4f> qz(a, b);
  std::cout << a << std::endl << std::endl;
  std::cout << b << std::endl << std::endl;
  std::cout << qz.matrixS() << std::endl << std::endl;
  std::cout << qz.matrixT() << std::endl << std::endl;

Gives me this output:

   0.68    0.823   -0.444    -0.27
 -0.211   -0.605    0.108   0.0268
  0.566    -0.33  -0.0452    0.904
  0.597    0.536    0.258    0.832

-0.271  -0.967  -0.687   0.998
 0.435  -0.514  -0.198  -0.563
-0.717  -0.726   -0.74  0.0259
 0.214   0.608  -0.782   0.678

 0.1081    1.316   0.7939  -0.1046
 0.1319  0.04982  -0.4088  -0.1396
      0        0  -0.4452  -0.1534
     -0        0        0   -1.292

  1.571   0.1476  -0.3816  -0.7557
      0  -0.8758   -0.577  -0.2295
      0        0   -1.174  -0.2325
      0       -0        0   0.7569

The equivalent call to qz() in Matlab gives me these two matrices for S & T:

   -1.1402   -0.6674   -0.1892    0.0871
         0    0.1738    1.2628    0.8126
         0   -0.1565   -0.0653    0.4159
         0         0         0    0.4559

    0.6679   -0.1078   -0.1005    0.0189
         0    1.7314         0   -0.5382
         0         0    0.8791    0.6238
         0         0         0    1.2019

Eigen gives the following as the Eigenvalues:

(0.05029,0.266)  (0.05029,-0.266)       (0.3793,-0)        (-1.707,0)

Whereas Matlab gives:

 -1.7072 + 0.0000i   0.0131 - 0.3496i   0.0131 + 0.3496i   0.3793 + 0.0000i

The two real eigenvalues agree, but the complex ones are different. I'm stumped as to whether it is the QZ factorization or the current code in GeneralizedEigenSolver which is wrong.

Apologies if I have missed anything obvious.
Comment 12 Gael Guennebaud 2016-06-08 08:36:35 UTC
hm, indeed. Numerically speaking, the Q,Z,S,T factors satisfy all required properties, but the complex eigenvalue pair does not satisfy:

    det(A - lambda B)=0

so we do have a bug somewhere, either in QZ where a slight lost in precision break the eigenvalue extraction, or in the eigenvalue extraction itself...
Comment 13 Gael Guennebaud 2016-06-08 09:26:12 UTC
The problem is indeed in the extraction where there is a big shortcoming! I'll apply a fix soon.
Comment 14 Gael Guennebaud 2016-06-08 14:42:29 UTC
Fixed: https://bitbucket.org/eigen/eigen/commits/9098c2fcf5d8/
and backported to 3.2.


Let me known if guaranteeing that the respective 2x2 diagonal block of T of a 2x2 non triangular diagonal block of S is diagonal (instead of being triangular only) would help the eigenvector extraction. Then we'll update RealQZ accordingly.
Comment 15 Gael Guennebaud 2016-06-10 12:16:39 UTC
Finally, I've updated RealQZ to diagonalize the respective 2x2 blocks of T:

https://bitbucket.org/eigen/eigen/commits/eb8d87ca2cfc/
https://bitbucket.org/eigen/eigen/commits/3a0db566465c/

That's simpler, more stable, and in par with Matlab/Lapack
Comment 16 Christoph Hertzberg 2016-08-03 10:48:07 UTC
*** Bug 1269 has been marked as a duplicate of this bug. ***
Comment 17 Tobias Wood 2016-08-18 15:15:12 UTC
Created attachment 720 [details]
Patch 3

Hello,

Now Gael has fixed the eigenvalues, I think this new patch successfully calculates the eigenvectors, both real and complex. I added a new test file.

Gael - previously you suggested that I replace A.transpose().cwiseProduct(B).sum() with A.dot(B), which I would very much like to do. However, this is giving me a compile error:

/Users/Tobias/Code/eigen_repo/Eigen/src/Core/Dot.h:71:3: error: no member named
      'YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX' in 'Eigen::internal::static_assertion<false>'
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)

which seems very strange?
Comment 18 Gael Guennebaud 2016-08-23 16:18:35 UTC
Thank you very much for the patch!! I applied it (with some adjustments):

https://bitbucket.org/eigen/eigen/commits/c0e7a361ad5a/

and made some minor improvements:

https://bitbucket.org/eigen/eigen/commits/fea26b6579c9/
Comment 19 Tobias Wood 2016-08-23 19:10:54 UTC
I'm glad to have helped.

Out of curiosity, do you understand why replacing A.transpose().cwiseProduct(B).sum() with A.dot(B) is causing that error? It looks like it ought to work to me.
Comment 20 Gael Guennebaud 2016-08-23 19:22:23 UTC
This is because here B is a mat.block(...) expression which is not a vector at compile time. Replacing it by something like mat.row(.).segment(..) would make it work (but not for complexes as a.dot(b) <=> a.adjoint()*b, not a.transpose()*b). Anyway, I've fused them as matrix-vector products.

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