eigensolver in eigen Library

4.5k views Asked by At

I would like to translate [vec,val] = eig(A) from MATLAB to c++ using Eigen library, but I couldn't reach to the same result! I tried eigensolver,ComplexEigenSolver and SelfAdjointEigenSolver. None of them give me the result like eig(A) in MATLAB.

Sample matrices:
Tv(:,:,223) =

    0.8648   -1.9658   -0.2785
   -1.9658    4.9142    0.8646
   -0.2785    0.8646    0.3447


Tv(:,:,224) =

    1.9735   -0.4218    1.0790
   -0.4218    3.3012    0.1855
    1.0790    0.1855    3.7751


Tv(:,:,225) =

    2.4948    1.0185    1.1633
    1.0185    1.1732   -0.4479
    1.1633   -0.4479    4.3289


Tv(:,:,226) =

    0.3321    0.0317    0.1617
    0.0317    0.0020   -0.0139
    0.1617   -0.0139    0.5834

Eigen:

MatrixXcd vec(3 * n, 3);
VectorXcd val(3);
for (int k = 0; k < n; k++){
        EigenSolver<Matrix3d> eig(Tv.block<3, 3>(3 * k, 0));
        vec.block<3, 3>(3 * k, 0) = eig.eigenvectors();
        cout <<endl << vec.block<3, 3>(3 * k, 0) << endl;
        val = eig.eigenvalues();
        cout << "val= " << endl << val << endl;

    }

//results

  (0.369152,0)   (-0.830627,0)   (-0.416876,0)
  (-0.915125,0)   (-0.403106,0) (-0.00717218,0)
  (-0.162088,0)    (0.384142,0)   (-0.908935,0)
val=
  (5.86031,0)
(0.0396418,0)
 (0.223765,0)

 (0.881678,0)  (0.204005,0)  (0.425472,0)
  (0.23084,0)  (-0.97292,0) (-0.011858,0)
(-0.411531,0) (-0.108671,0)  (0.904894,0)
val=
(1.35945,0)
(3.41031,0)
(4.27996,0)

 (0.526896,0) (-0.726801,0)  (0.440613,0)
(-0.813164,0) (-0.581899,0) (0.0125466,0)
(-0.247274,0)  (0.364902,0)  (0.897609,0)
val=
(0.377083,0)
 (2.72623,0)
 (4.89367,0)

    (0.88992,0)    (-0.43968,0)    (0.121341,0)
    (0.13406,0) (-0.00214387,0)   (-0.990971,0)
   (-0.43597,0)   (-0.898152,0)  (-0.0570358,0)
val=
   (0.257629,0)
   (0.662467,0)
(-0.00267575,0)

MATLAB:

for k=1:n
    [u,d] = eig(Tv(:,:,k))
end

%results

u =

    0.8306   -0.4169   -0.3692
    0.4031   -0.0072    0.9151
   -0.3841   -0.9089    0.1621


d =

    0.0396         0         0
         0    0.2238         0
         0         0    5.8603


u =

    0.8817    0.2040    0.4255
    0.2308   -0.9729   -0.0119
   -0.4115   -0.1087    0.9049


d =

    1.3594         0         0
         0    3.4103         0
         0         0    4.2800


u =

   -0.5269    0.7268    0.4406
    0.8132    0.5819    0.0125
    0.2473   -0.3649    0.8976


d =

    0.3771         0         0
         0    2.7262         0
         0         0    4.8937


u =

   -0.1213   -0.8899    0.4397
    0.9910   -0.1341    0.0021
    0.0570    0.4360    0.8982


d =

   -0.0027         0         0
         0    0.2576         0
         0         0    0.6625

What's your suggestion?

2

There are 2 answers

6
Ander Biguri On

Well.... the results are the same....

Result eigen:

  (0.369152,0)   (-0.830627,0)   (-0.416876,0)
  (-0.915125,0)   (-0.403106,0) (-0.00717218,0)
  (-0.162088,0)    (0.384142,0)   (-0.908935,0)
val=
  (5.86031,0)
(0.0396418,0)
 (0.223765,0)

result matlab:

u =

    0.8306   -0.4169   -0.3692
    0.4031   -0.0072    0.9151
   -0.3841   -0.9089    0.1621


d =

    0.0396         0         0
         0    0.2238         0
         0         0    5.8603

I have good news....

The vectors are THE SAME, but unordered.....

eigV1 from eigen is -eigV3 from Matlab,

eigV2 from eigen is -eigV1 from Matlab,

eigV3 from eigen is -eigV2 from Matlab,

The eigenvalues are reordered equally....

0
ggael On

I don't get your question, as looking at your results they all returns the same. Recall that the eigen-decomposition of a matrix is not completely unique:

  • eigenvalues/vectors can be arbitrarily reordered
  • if v is an eigenvector, then -v is also a valid eigenvector

Since your matrices are symmetric, you should use SelfAdjointEigenSolver to get them automatically ordered as MatLab. Then the eigenvectors will only differs from their sign, but you will have to live with that.