Issue with low-rank matrices

[jdbancal]

Consider the rank-one 5x5 matrix

a = 

   15.035447086947079479    3.932587856183598677   -4.848070276813470542   -8.027254936523050904   -2.865327349780228231
    3.932587856183598677    1.028585791773944732   -1.268034278346991263   -2.099564123322002035   -0.749439073848281425
   -4.848070276813470542   -1.268034278346991263    1.563224909309606855    2.588329820664053864    0.923903910371237535
   -8.027254936523050904   -2.099564123322002035    2.588329820664053864    4.285660509016328222    1.529765824738644411
   -2.865327349780228231   -0.749439073848281425    0.923903910371237535    1.529765824738644411    0.546049663433429209

A direct computation shows that its eigenvalues are

   22.458967960480388496
                       0
                       0
                       0
                       0

with first eigenvector

   0.818207159789096
   0.214005710759507
  -0.263824932422074
  -0.436831537135830
  -0.155926946446774

However, using Spectra::SymEigsSolver on Spectra::DenseSymMatProd, Spectra seems to find the eigenvalue

275.20814154032755067

The same result is found when describing the matrix in a sparse way, and calling Spectra::SymEigsSolver on Spectra::SparseSymMatProd. This was not an issue on earlier versions of Spectra.

[jdbancal]

The example above was actually when computing with 50 significant digits, so incompletely formulated. Here is an example obtained when computing with 10 significant digits

>> gem.workingPrecision(10);
>> v = gem([-0.7821975807         
   3.909983147          
  -1.699076588          
  -3.325588409          
  -1.371499025]);

>> a = v*v'

a = 

    0.6118330552           -3.058379358             1.329013596             2.601267208             1.072783220         
   -3.058379358            15.28796821             -6.643360824           -13.00299463             -5.362538075         
    1.329013596            -6.643360824             2.886861251             5.650429406             2.330281884         
    2.601267208           -13.00299463              5.650429406            11.05953826              4.561041261         
    1.072783220            -5.362538075             2.330281884             4.561041261             1.881009576         

>> eigs(a,[],1)

ans = 

  -25093.27364    

But the largest eigenvalue is actually

>> eigs(double(a),[],1)

ans =

   31.7272

[jdbancal]

And here is an example with 15 significant digits (i.e. with a precision comparable of that of double type):

>> gem.workingPrecision(15)
>> v = gem([4.21544269694397     
   2.53908041061553     
  -4.54658610267737     
  -0.996463928127838    
  -2.64328807668127]);
>> a = v*v'

a = 

   17.7699571312182        10.7033479738827       -19.1658731825582        -4.20053658859459      -11.1426294187651     
   10.7033479738827         6.44692933157151      -11.5441477084849        -2.53010203979439       -6.71152097511499    
  -19.1658731825582       -11.5441477084849        20.6714451890590         4.53050904744533       12.0179368348118     
   -4.20053658859459       -2.53010203979439        4.53050904744533        0.992940360059961       2.63394122006329    
  -11.1426294187651        -6.71152097511499       12.0179368348118         2.63394122006329        6.98697185632535    

>> eigs(a,[],1)

ans = 

   823.677714297095     

>> eigs(double(a),[],1)

ans =

   52.8682

[jdbancal]

With Spectra 0.9, the two examples above produce the correct larges eigenvalues (respectively 31.7272 and 52.8682).

[yixuan]

Thanks for reporting. This should have been fixed in https://github.com/yixuan/spectra/commit/62e43ce66c96652f0e53c9fa824e48b8a88c6bb5. I have included a test in https://github.com/yixuan/spectra/blob/master/test/Example2.cpp.

[jdbancal]

Looks good, thanks!