Generalized Eigenproblem where B is not positive definite

[abhijit-c]

Hello all,

In the last steps of the implementation of my paper, it asks by to solve a generalized eigenproblem for the smallest eigenvalue/eigenvector pair. This roughly looks like

\[
D^T W(x)D u(x) = \lambda C u(x)
\]

Without digressing into too many details \( D^t W(x) D \) has a good self adjoint structure, but C is not positive semidefinite. As such I cannot use igl's eigs function or Eigen's GeneralizedSelfAdjointEigenSolver. Eigen has a generic GeneralizedEigenSolver which accepts right hand sides like mine, but it seems like they have not implemented a way to get the eigenvectors yet from the decomposition.

Any pointers on what I should do? Try to implement my own GeneralizedEigenSolver or find some other library? This is a rather important final portion of the paper. If it helps here is the matrix C:

0	0	0	0	-2
0	1	0	0	0
0	0	1	0	0
0	0	0	1	0
-2	0	0	0	0

[abhijit-c]

Oh I'm truly crazy, I should just turn the problem into \( (C^{-1})D^T W(x) D u(x) = \lambda u(x) \). \( C \) is invertible and small so there should be no issues. So then I just have an eigenvalue problem, however the LHS's symmetric structure is ruined. I only need the smallest eigenvalue / eigenvector pair, so I should inverse power method it, correct? Or maybe since the size of the problem is small (5x5) we can just directly compute all of the eigen(vectors/values) and then get the corresponding ones?

[jpanetta]

You should be able to call GeneralizedEigenSolver::eigenvectors(); this method doesn't appear in the documentation, but it exists.

It seems it's still an open problem to develop an efficient algorithm that can exploit the symmetry in your case. Also, notice that the eigenvalues/eigenvectors can be complex even though your matrices are Hermitian.

[abhijit-c]

I took a look at the src and indeed it seems to be implemented. However, I'm directly recivieng an error which states that it doesn't exist? Here is the related code/error, maybe someone can spot something I'm missing? My Eigen version 3.3.4

//Solve Generalized Eigenproblem.
Eigen::MatrixXd LHS = DT*W*D;
Eigen::GeneralizedEigenSolver<Eigen::MatrixXd> ges;
ges.compute(LHS,C);
cout << ges.eigenvalues() << endl;
cout << ges.eigenvectors() << endl;

[ 96%] Building CXX object CMakeFiles/ex5_bin.dir/src/main.cpp.o
/home/trostaft/Git/GM_Final_Project/src/main.cpp: In function ‘bool key_down(igl::opengl::glfw::Viewer&, unsigned char, int)’:
/home/trostaft/Git/GM_Final_Project/src/main.cpp:328:15: error: ‘class Eigen::GeneralizedEigenSolver<Eigen::Matrix<double, -1, -1> >’ has no member named ‘eigenvectors’; did you mean ‘m_eigenvectorsOk’?
   cout << ges.eigenvectors() << endl;
               ^~~~~~~~~~~~
               m_eigenvectorsOk
make[2]: *** [CMakeFiles/ex5_bin.dir/build.make:63: CMakeFiles/ex5_bin.dir/src/main.cpp.o] Error 1
make[1]: *** [CMakeFiles/Makefile2:70: CMakeFiles/ex5_bin.dir/all] Error 2
make: *** [Makefile:84: all] Error 2

[jpanetta]

No idea... the following compiles for me with Eigen 3.3.4:

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <iostream>

int main(int argc, const char *argv[]) {
    Eigen::MatrixXd A, B;
    A.setIdentity(3, 3);
    B.setIdentity(3, 3);
    Eigen::GeneralizedEigenSolver<Eigen::MatrixXd> ges;
    ges.compute(A, B);
    std::cout << ges.eigenvalues() << std::endl;
    std::cout << ges.eigenvectors() << std::endl;
    return 0;
}

It only computes one eigenvector correctly, but maybe it's just confused by the repeated eigenvalue.

[abhijit-c]

Very strange, writing a file separate from the final project environment and manually linking seems to work perfectly, yet using the final project makefile doesn't work. I'll investigate the build system and see what I find.

[abhijit-c]

I found the issue. It seems that the version of eigen packaged with libigl does not have the eigenvectors() function yet implemented. Looking at the source here is what's in place:

//template<typename MatrixType>
//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
//{
//  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
//  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
//  Index n = m_eivec.cols();
//  EigenvectorsType matV(n,n);
//  // TODO
//  return matV;
//}

This is mildly frustrating, I doubt naively replacing the new eigen version with the currently used one in libigl will work. Any ideas on how to do that, or maybe I'm better off finding an external library?

[bmahlbrand]

build eigen to a DLL and add an alternative version that you call for this? Not sure that I really recommend that since it's bound to have some tedious flaw.

[abhijit-c]

Oof, I really don't want to deal with that. Maybe I'll try solving the eigenvalue problem \( (C^{-1} ) D^T W(x) D u(x) = \lambda u(x) \) first, and see if that works for my applications. It's not like the generalized solver will be able to take advantage of the symmetric structure anyway since C is indefinite.

[jiangzhongshi]

Eigen is header only, a drop in replacement should work with some chance. There are components in libigl not fully compatible with 3.3.4 but maybe you won’t encounter it.

…

On Tue, May 8, 2018 at 7:07 AM Abhijit Chowdhary ***@***.***> wrote: Oof, I really don't want to deal with that. Maybe I'll try solving \( (C^{-1} ) D^T W(x) D u(x) = \lambda u(x) \) first, and see if that works for my applications. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <https://github.com/danielepanozzo/gp/issues/41#issuecomment-387367423>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AMKxFzQ29nkALsnlEcdMvySWO9HUckcMks5twXx7gaJpZM4T2SVu> .

[jpanetta]

Especially for such a small problem, I don't think there's anything wrong with transforming it into an ordinary, asymmetric eigenvalue problem as you wrote. Note that C^-1 just exchanges the first and last rows of D^T W D and multiplies them by -1/2.