Strange behaviour using Arpack EigenLibrary Wrapper for small matrices

Question

I want to solve the following generalized EVP using the ArpackWrapper of the Eigen library:

K_e is SPD. Normally K_g is indefinite and singular but for this MVP it is just indefinite. Furthermore, I'm interested in the smallest eigenvalues. For large systems 300kx300k i obtained reasonable results but for this small example the results appear strange. The first 4 resulting eigenvalues read

 0.8987
-0.720851 
 0.607632 
 0.729297

For less requested eigenvalues solverEig.info() returns Eigen::NoConvergence.

If i use Maple to compute the eigenvalues i get

 522.991427951073 
-175.66558721639944
  66.23707710214939
  -7.756864956770603
 355.6461914072188

which can be inserted in the problem statement to see their correctness.

Therefore, my question is why does it not yield the correct result? Is this a problem of arpack or of the eigen wrapper? Or most probably my wrong usage/understanding of arpack or the eigen wrapper.

Versions:

Eigen 3.3.7

Arpack https://github.com/opencollab/arpack-ng different versions result to the same behaviour

The code:

#include <iostream>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include "eigen3/unsupported/Eigen/ArpackSupport"
#include "eigen3/unsupported/Eigen/SparseExtra"
using namespace std;
int main() {

    Eigen::SparseMatrix<double> ke; 
    Eigen::SparseMatrix<double> kg; 

    Eigen::loadMarket(ke,"ke.txt"); //5x5 Matrix
    Eigen::loadMarket(kg,"kg.txt"); //5x5 Matrix


    Eigen::ArpackGeneralizedSelfAdjointEigenSolver<Eigen::SparseMatrix<double>> solverEig;

    solverEig.compute(ke,-kg,4,"SM",Eigen::ComputeEigenvectors);


    cout<<solverEig.eigenvalues()<<endl;
    return 0;
}

The files:

Ke.txt

%%MatrixMarket matrix coordinate  real general
5 5 19
1 1 3.2621670111997303820317029021680355072021484375000000000000000000e+01
2 1 1.4310835055998653686515353911090642213821411132812500000000000000e+01
3 1 -2.0000000000000000000000000000000000000000000000000000000000000000e+00
5 1 -1.4310835055998653686515353911090642213821411132812500000000000000e+01
1 2 1.4310835055998653686515353911090642213821411132812500000000000000e+01
2 2 8.7155417527999333060506614856421947479248046875000000000000000000e+01
3 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 2 -7.1554175279993268432576769555453211069107055664062500000000000000e+00
1 3 -2.0000000000000000000000000000000000000000000000000000000000000000e+00
2 3 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 3 8.2000000000000000000000000000000000000000000000000000000000000000e+01
4 3 -8.0000000000000000000000000000000000000000000000000000000000000000e+01
3 4 -8.0000000000000000000000000000000000000000000000000000000000000000e+01
4 4 8.0000000000000000000000000000000000000000000000000000000000000000e+01
5 4 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 5 -1.4310835055998653686515353911090642213821411132812500000000000000e+01
2 5 -7.1554175279993268432576769555453211069107055664062500000000000000e+00
4 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 5 8.7155417527999333060506614856421947479248046875000000000000000000e+01

Kg.txt

%%MatrixMarket matrix coordinate  real general
5 5 19
1 1 -2.1389117730512491322159007722802925854921340942382812500000000000e-01
2 1 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 1 1.8999822924102233168142106478626374155282974243164062500000000000e-01
5 1 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
2 2 -2.1389117730512491322159007722802925854921340942382812500000000000e-01
3 2 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 2 -9.3693312577685552988704387189500266686081886291503906250000000000e-02
1 3 1.8999822924102233168142106478626374155282974243164062500000000000e-01
2 3 0.0000000000000000000000000000000000000000000000000000000000000000e+00
3 3 -3.8974822924103957877406401166808791458606719970703125000000000000e-01
4 3 1.9975000000001727484821856251073768362402915954589843750000000000e-01
3 4 1.9975000000001727484821856251073768362402915954589843750000000000e-01
4 4 1.8161606213113259955527212241577217355370521545410156250000000000e-01
5 4 0.0000000000000000000000000000000000000000000000000000000000000000e+00
1 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
2 5 -9.3693312577685552988704387189500266686081886291503906250000000000e-02
4 5 0.0000000000000000000000000000000000000000000000000000000000000000e+00
5 5 4.7505937470883541351440726430155336856842041015625000000000000000e-01

Answer 1

Maybe in your case the parameter NCV of ARPACK is by default NEV (the nb of eigen vector to be computed), or the size of the matrix N, which are both small while it is recommanded that this parameter should not be less than like 20 (I'd put at least 70 or 100 really..). See those links :

https://www.caam.rice.edu/software/ARPACK/UG/node136.html (NCV, NEV and remark 4)

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html#scipy.sparse.linalg.eigs (ncv description)

Answer 2

It's odd that you get real eigenvalues, not a complex pair -- did you perhaps run eigsh for symmetric problems, which looks only at the lower halves ?

from scipy.sparse.linalg import eigs  # arpack

v0 = np.ones( K.shape[0] )  # else random

evals, V = eigs( A=K, M=M, k=3, tol=0, v0=v0 )  # k < n-1
print( "eigs evals:", evals )
eigcheck( K, evals, V, M=M )

K: 
 [[32.6217 14.3108 -2 0 -14.3108]
 [14.3108 87.1554 0 0 -7.15542]
 [-2 0 82 -80 0]
 [0 0 -80 80 0]
 [-14.3108 -7.15542 0 0 87.1554]]
M: 
 [[-0.213891 0 0.189998 0 0]
 [0 -0.213891 0 0 -0.0936933]
 [0.189998 0 -0.389748 0.19975 0]
 [0 0 0.19975 0.181616 0]
 [0 -0.0936933 0 0 0.475059]]
eigs evals: [ -1301.07 -2136.22j   -1301.07 +2136.22j   36.0965 -2477.99j ]
eigcheck : |Av - λv|_max [2e-13 2e-13 2e-13]

Could you point me to similar big problems on the web ? Thanks