Greetings!
I've encountered a rather curious issue (heretofore, Eigen has served me perfectly without any problems!).
I am trying to solve a linear eigenvalue problem (the Schrodinger equation, precisely) :
H psi = E psi; Where H is an NxN matrix (Self-adjoint, real, &c).
I defined it as such, and then called on :

Code: Select all

SelfAdjointEigenSolver<MatrixXd> solver

, and requested it to "compute". I then located the desired eigenvalue, and extracted the eigenvector, only to find that it is completely wrong!
Here's a comparison of the result from Mathematica with that of Eigen's:
Eigen:
Mathematica:
(Obviously, Mathematica's output is the correct one; and Eigen's behaviour here is particularly strange: the entries of the vector seem to oscillate between positive and negative -- an utterly improbable outcome).
How do I resolve this problem?
Thank you for your attention,
Daniel
P.S.
H was populated identically in Mathematica & Eigen. All elements match.

In order to help you, we need a way to reproduce the issue: either the input matrix or a self-contained example.

Hi!
Thanks for your reply.
Here's the link for it -- I've had to compress it (the matrix itself is very large (as you might imagine). It's also (rather) sparse).
http://speedy.sh/nm39y/HMatrix.tar.bz2
Please let me know if there's anything else I can provide to resolve this issue.
Thankful,
Daniel

Furthermore, to ease matters, here's the code that populates the matrix:

Code: Select all

H(i, i-1) = -0.5/(dr*dr);
H(i, i) = 1.0/(dr*dr)-2/(r0+i*dr);
H(i, i+1) = H(i, i-1)


(A typical tridiagonal matrix);
"i" goes from 0 to N=2000; dr = 7.5*10^(-3); r0 = 0.001;

I get the same eigenvalues than mathlab:

Code: Select all

-15804.1  -1.97059  -0.49633 -0.220583  -0.09764 0.0545751  0.261432   0.51913  0.825403   1.17894   1.57891 .... 8877.56   8878.51   8879.42   8880.28    8881.1   8881.88   8882.61    8883.3   8883.95   8884.55   8885.12   8885.64   8886.12   8886.55   8886.95    8887.3   8887.62   8887.89   8888.13   8888.33    8888.5   8888.64


Code to reproduce:

Code: Select all

#include <iostream>
#include <fstream>
#include <Eigen/Dense>
using namespace Eigen;
int main(int argc, char **argv) {
  std::ifstream f(argv[1]);
  int n = 1000;
  MatrixXd M(n,n);
  for(int i=0;i<n;++i)
    for(int j=0;j<n;++j)
    {
      f >> M(i,j);
      f.ignore(1);
    }
 
  SelfAdjointEigenSolver<MatrixXd> eig;
 
  eig.compute(M);
  // version optimized for tridiagonal matrices (devel branch only):
  // eig.computeFromTridiagonal(M.diagonal(), M.diagonal<1>());
 
  std::cout << eig.eigenvalues().transpose() << "\n";
}


Hi, thanks for your reply.
The problem is specifically with the eigenvectors. Take a look at the eigenvector corresponding to the 2nd eigenvalue (-1.97..). If you compare Eigen's with Mathematica's (or Matlab's) you'll end with two fundamentally different vectors.
As you can see in my original post, the vector itself contains entries which simply oscillate between positive and negative -- i.e. a nonsensical result.
Can you reproduce the same result, on your end?
Beholden for your help,
Daniel

Note that if V is an eigenvector, then -V is also one.

edit:
for the record, here are the first entries of the first eigenvectors, Eigen (i.e., eig.eigenvectors().topLeftCorner(30,6)):

Code: Select all

    0.993726  -0.00122545 -0.000435866  0.000240444  0.000219208  0.000263186
    0.111145   0.00857709   0.00305096  -0.00168308  -0.00153445  -0.00184231
     0.01239     0.017876   0.00635663  -0.00350646  -0.00319671  -0.00383795
  0.00137964    0.0266563   0.00947365  -0.00522535  -0.00476354  -0.00571875
 0.000153538    0.0349283    0.0124037  -0.00684047  -0.00623551  -0.00748528
 1.70814e-05    0.0427082    0.0151508  -0.00835384  -0.00761438  -0.00913953
  1.8999e-06    0.0500142    0.0177198  -0.00976793  -0.00890233    -0.010684
 2.11286e-07    0.0568649    0.0201156   -0.0110854   -0.0101018   -0.0121216
  2.3494e-08    0.0632791    0.0223438   -0.0123091   -0.0112151    -0.013455
 2.61216e-09    0.0692752    0.0244094   -0.0134417   -0.0122449   -0.0146872
 2.90398e-10    0.0748712    0.0263179   -0.0144861   -0.0131937   -0.0158213
 3.22765e-11    0.0800848    0.0280745   -0.0154451   -0.0140639   -0.0168601
 3.58493e-12    0.0849327    0.0296844   -0.0163215   -0.0148581   -0.0178066
 3.97285e-13    0.0894316    0.0311525   -0.0171179   -0.0155787   -0.0186636
 4.37353e-14    0.0935972    0.0324838    -0.017837   -0.0162282   -0.0194342
 4.73083e-15    0.0974449    0.0336833   -0.0184815   -0.0168088   -0.0201209
  4.9135e-16     0.100989    0.0347556    -0.019054   -0.0173229   -0.0207268
 4.69477e-17     0.104245    0.0357054   -0.0195569   -0.0177729   -0.0212544
  3.8354e-18     0.107225    0.0365372   -0.0199927    -0.018161   -0.0217065
 2.35336e-19     0.109943    0.0372554   -0.0203639   -0.0184893   -0.0220858
   8.056e-21     0.112411    0.0378643   -0.0206729   -0.0187601   -0.0223948
           0     0.114642    0.0383682   -0.0209219   -0.0189754   -0.0226361
           0     0.116647     0.038771   -0.0211132   -0.0191374   -0.0228122
           0     0.118438    0.0390768   -0.0212491   -0.0192481   -0.0229255
           0     0.120025    0.0392893   -0.0213317   -0.0193094   -0.0229785
           0     0.121418    0.0394125    -0.021363   -0.0193232   -0.0229735
           0     0.122627    0.0394499   -0.0213452   -0.0192916   -0.0229128
           0     0.123663    0.0394051   -0.0212803   -0.0192163   -0.0227986
           0     0.124533    0.0392815   -0.0211703   -0.0190992   -0.0226333
           0     0.125247    0.0390826   -0.0210169    -0.018942   -0.0224189


Lapack/MatLab (i.e., V(1:30,1:6))

Code: Select all

 9.9373e-01  -1.2255e-03  -4.3587e-04   2.4044e-04   2.1921e-04   2.6319e-04
   1.1115e-01   8.5771e-03   3.0510e-03  -1.6831e-03  -1.5344e-03  -1.8423e-03
   1.2390e-02   1.7876e-02   6.3566e-03  -3.5065e-03  -3.1967e-03  -3.8380e-03
   1.3796e-03   2.6656e-02   9.4737e-03  -5.2253e-03  -4.7635e-03  -5.7187e-03
   1.5354e-04   3.4928e-02   1.2404e-02  -6.8405e-03  -6.2355e-03  -7.4853e-03
   1.7081e-05   4.2708e-02   1.5151e-02  -8.3538e-03  -7.6144e-03  -9.1395e-03
   1.8999e-06   5.0014e-02   1.7720e-02  -9.7679e-03  -8.9023e-03  -1.0684e-02
   2.1129e-07   5.6865e-02   2.0116e-02  -1.1085e-02  -1.0102e-02  -1.2122e-02
   2.3494e-08   6.3279e-02   2.2344e-02  -1.2309e-02  -1.1215e-02  -1.3455e-02
   2.6122e-09   6.9275e-02   2.4409e-02  -1.3442e-02  -1.2245e-02  -1.4687e-02
   2.9042e-10   7.4871e-02   2.6318e-02  -1.4486e-02  -1.3194e-02  -1.5821e-02
   3.2285e-11   8.0085e-02   2.8075e-02  -1.5445e-02  -1.4064e-02  -1.6860e-02
   3.5890e-12   8.4933e-02   2.9684e-02  -1.6321e-02  -1.4858e-02  -1.7807e-02
   3.9895e-13   8.9432e-02   3.1152e-02  -1.7118e-02  -1.5579e-02  -1.8664e-02
   4.4345e-14   9.3597e-02   3.2484e-02  -1.7837e-02  -1.6228e-02  -1.9434e-02
   4.9291e-15   9.7445e-02   3.3683e-02  -1.8482e-02  -1.6809e-02  -2.0121e-02
   5.4785e-16   1.0099e-01   3.4756e-02  -1.9054e-02  -1.7323e-02  -2.0727e-02
   6.0886e-17   1.0424e-01   3.5705e-02  -1.9557e-02  -1.7773e-02  -2.1254e-02
   6.7652e-18   1.0722e-01   3.6537e-02  -1.9993e-02  -1.8161e-02  -2.1707e-02
   7.5120e-19   1.0994e-01   3.7255e-02  -2.0364e-02  -1.8489e-02  -2.2086e-02
   8.3258e-20   1.1241e-01   3.7864e-02  -2.0673e-02  -1.8760e-02  -2.2395e-02
   9.1836e-21   1.1464e-01   3.8368e-02  -2.0922e-02  -1.8975e-02  -2.2636e-02
   1.0017e-21   1.1665e-01   3.8771e-02  -2.1113e-02  -1.9137e-02  -2.2812e-02
   1.0671e-22   1.1844e-01   3.9077e-02  -2.1249e-02  -1.9248e-02  -2.2926e-02
   1.0868e-23   1.2002e-01   3.9289e-02  -2.1332e-02  -1.9309e-02  -2.2978e-02
   1.0235e-24   1.2142e-01   3.9412e-02  -2.1363e-02  -1.9323e-02  -2.2973e-02
   8.4857e-26   1.2263e-01   3.9450e-02  -2.1345e-02  -1.9292e-02  -2.2913e-02
   5.7508e-27   1.2366e-01   3.9405e-02  -2.1280e-02  -1.9216e-02  -2.2799e-02
   2.8076e-28   1.2453e-01   3.9282e-02  -2.1170e-02  -1.9099e-02  -2.2633e-02
   7.3624e-30   1.2525e-01   3.9083e-02  -2.1017e-02  -1.8942e-02  -2.2419e-02


Hi Ggael!
I've finally managed to resolve the problem. I can't thank you enough for your time and help. The matter was due to a ridiculous mistake of mine, in unloading the eigenvectors into a dedicated file. So that's been cleared. Thanks again!
One more question, if I may: is there a dedicated routine for the eigenvalues/vectors of a sparse matrix? (I am aware that Armadillo/Lapack has those, not sure if Eigen's already implemented that)
Once again,
Bound and grateful for your attention and aid,
Daniel

For tridiagonal, yes, as I showed in a previous example (see above):

Code: Select all

eig.computeFromTridiagonal(M.diagonal(), M.diagonal<1>());


For general sparse matrices, we have a wrapper to ARPACK.