Summary:  GeneralizedEigenSolver: missing computation of eigenvectors  

Product:  Eigen  Reporter:  Gael Guennebaud <gael.guennebaud>  
Component:  Eigenvalues  Assignee:  Nobody <eigen.nobody>  
Status:  RESOLVED FIXED  
Severity:  Feature Request  CC:  chtz, daniel.perry, gael.guennebaud, jacob.benoit.1, jitseniesen, tobias, tony.williams.1990  
Priority:  Normal  Keywords:  JuniorJob  
Version:  3.2  
Hardware:  All  
OS:  All  
Whiteboard:  
Bug Depends on:  
Bug Blocks:  814  
Attachments: 

Description
Gael Guennebaud
20130810 23:35:25 UTC
Hi  I was curious if there are any current plans to address this bug? AFAIK, nobody started working on it, so the job is open ;) 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/generalizedeigensolvernotgivingeigenvectorsineigenlibrary), 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. 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? Hi, thanks for working on this! Regarding unit tests: You should generally avoid hardcoded 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). Created attachment 710 [details]
Preliminary Patch
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!
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 on attachment 711 [details] Patch 2 Great start! Some cleanup 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 = i1; j >= 0; j) { >+ const Index st = j+1; >+ const Index sz = ij; 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() Also, handling non SPD matrices is a requirement. For those we already have the GeneralizedSelfAdjointEigenSolver class which is much faster (perhaps less robust though). 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. 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... The problem is indeed in the extraction where there is a big shortcoming! I'll apply a fix soon. 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. 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 *** Bug 1269 has been marked as a duplicate of this bug. *** 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?
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/ 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. 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 matrixvector products.  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/645. 