I'm trying to solve eigenproblem using Eigen library. I am using SelfAdjointEigenSolver to get eigenvalues and eigenvectors. While I seem to have quite good proof that eigenvalues are alright, eigenvectors seem to not be. Or I'm understanding something wrong.

After calculating eigenvector vct, I crate another one vct2=matrix*vct. The very idea of eigenproblem is that vct2=lambda*vct, where lambda is eigenvalue. So, I printed every coefficient of vct2 divided by vct to see if it is constant. It is not. Not even close. What is wrong? In here I create a random matrix, solve it's eigenproblem, take the first eigenvector, and print the division.

For matrix that I got, the output was: 9.89816 4.10161 -0.765191 -2.27482 3.58762 0.710447 -4.79855 -2.88263 -3.05201 -3.05561

which does not look like a string of the same number. The code to generate problem:

#include<iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace Eigen;
using namespace std;
int main(int argc, char **argv)
{
    SelfAdjointEigenSolver<MatrixXd> es;

        MatrixXd matrix=MatrixXd::Random(10,10);
        espute(matrix);
        double E0=es.eigenvalues()[0];
        VectorXd vct=es.eigenvectors().col(0);
        VectorXd vct2=matrix*vct;
        for(int i=0;i<10;i++)
        {
            cout<<vct2(i)/vct(i)<<" ";
        }
    return 0;
}

Update:

Ok. I am unable to replicate the problem with simple solution. It is most likely that problem does not come from Eigen itself but indeed from my matrix. The fact that I've seen here similar problem with solving non-self adjoint matrix suggests that maybe I should check IF my matrix is self-adjoint. It should be, but I see no other than this.

Here is code with symmetrical real matrix (so, self-adjoint too). It does not show any similar problems.

#include<iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace Eigen;
using namespace std;
int main(int argc, char** argv)
{
    SelfAdjointEigenSolver<MatrixXd> es;

    MatrixXd a = MatrixXd::Random(30, 30);
    MatrixXd matrix = a + a.transpose();
    espute(matrix);
    double E0 = es.eigenvalues()[0];
    VectorXd vct = es.eigenvectors().col(0);
    VectorXd vct2 = matrix * vct;
    for (int i = 0; i < 10; i++)
    {
        cout << vct2(i) / vct(i) << " ";
    }
    return 0;
}

I'm trying to solve eigenproblem using Eigen library. I am using SelfAdjointEigenSolver to get eigenvalues and eigenvectors. While I seem to have quite good proof that eigenvalues are alright, eigenvectors seem to not be. Or I'm understanding something wrong.

After calculating eigenvector vct, I crate another one vct2=matrix*vct. The very idea of eigenproblem is that vct2=lambda*vct, where lambda is eigenvalue. So, I printed every coefficient of vct2 divided by vct to see if it is constant. It is not. Not even close. What is wrong? In here I create a random matrix, solve it's eigenproblem, take the first eigenvector, and print the division.

For matrix that I got, the output was: 9.89816 4.10161 -0.765191 -2.27482 3.58762 0.710447 -4.79855 -2.88263 -3.05201 -3.05561

which does not look like a string of the same number. The code to generate problem:

#include<iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace Eigen;
using namespace std;
int main(int argc, char **argv)
{
    SelfAdjointEigenSolver<MatrixXd> es;

        MatrixXd matrix=MatrixXd::Random(10,10);
        espute(matrix);
        double E0=es.eigenvalues()[0];
        VectorXd vct=es.eigenvectors().col(0);
        VectorXd vct2=matrix*vct;
        for(int i=0;i<10;i++)
        {
            cout<<vct2(i)/vct(i)<<" ";
        }
    return 0;
}

Update:

Ok. I am unable to replicate the problem with simple solution. It is most likely that problem does not come from Eigen itself but indeed from my matrix. The fact that I've seen here similar problem with solving non-self adjoint matrix suggests that maybe I should check IF my matrix is self-adjoint. It should be, but I see no other than this.

Here is code with symmetrical real matrix (so, self-adjoint too). It does not show any similar problems.

#include<iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
using Eigen::SelfAdjointEigenSolver;
using namespace Eigen;
using namespace std;
int main(int argc, char** argv)
{
    SelfAdjointEigenSolver<MatrixXd> es;

    MatrixXd a = MatrixXd::Random(30, 30);
    MatrixXd matrix = a + a.transpose();
    espute(matrix);
    double E0 = es.eigenvalues()[0];
    VectorXd vct = es.eigenvectors().col(0);
    VectorXd vct2 = matrix * vct;
    for (int i = 0; i < 10; i++)
    {
        cout << vct2(i) / vct(i) << " ";
    }
    return 0;
}

• c++
• eigen
• eigenvector

• 1 Try using Eigen::EigenSolver. SelfAdjointEigenSolver is only supposed to be used for selfadjoint matrices. For real matrices, this means that your matrix should be symmetric. – Ajay Commented Nov 15, 2024 at 20:50
• Ah, yes, that's true. While the problem still exists, since I've actualy seen it in matrix that was indeed selfadjoint, but it is wrong example. Thanks, I'll update it. – martes Commented Nov 15, 2024 at 21:13
• Does the problem still exist for your selfadjoint matrix? Could you please clarify? – Ajay Commented Nov 15, 2024 at 21:30
• Sorry, I was doing that on my work computer, I don't have eigen installed where I am, I will update in an hour or so, but for the whole matrix I was checking it was still broken, even though it was self adjoint. But I can't realy send that mess of a code as a "minimal case" – martes Commented Nov 15, 2024 at 21:33
• 2 Besides the already noted misuse of SelfAdjointEigenSolver in this example, computing the element-wise ratio between two vectors is not a good way to check if vct2 = lambda*vct. Better compute the difference vct2 - lambda*vct. – chtz Commented Nov 16, 2024 at 10:33

1 Answer 1

Ok, as chtz said the problem was actually in the way I was checking if the eigenvector is right. When calculating vct2(i)/vct(i) I was sometimes getting right answer, sometimes not. When using vct2-lambda*vct it showed that vector was always correct. Poking around a bit showed, that all the places where division was giving wrong answers, had coefficient of order 1e-17, so basicaly zeroes.

I was dividing zeroes by zeroes. The vectors themself were good.