2
$\begingroup$

I am trying to solve a generalized eigenvalue problem using Arpack, right now the code is using LAPACK but that's too slow, we only need a few eigenvalues and the matrices are sparse so using Arpack should be the way to go.

Before I start working with the original code I decided to test a simple case using scipy wrapper for Arpack (eigs) but the results that I am getting are wrong and change every time the code runs.

Minimum working example:

    import numpy as np
    from scipy.linalg import eig
    from scipy.sparse.linalg import eigs
    n = 8
    A = np.diag(np.arange(1,n+1,1.0))
    B = np.eye(n) # We want symmetric but a non-diagonal B. eigs gives correct answer for B=np.eye(n)
    B[0][n-1] = 2
    B[n-1][0] = 2


    evals,_ = eigs(A,k=3,M=B,which='LM') 
    print("The eigenvalues obtained by eigs (uses Arpack)")
    print(evals)

    print("Correct eigenvalues using eig (uses Lapack):")
    evals_l,_ = eig(A,b=B)
    print(evals_l)
```
$\endgroup$

2 Answers 2

8
$\begingroup$

The matrix B (M in the documentation) needs to positive definite according to the documentation: "If sigma is None, M is positive definite", this is in addition to the first requirement "M must represent a real, symmetric matrix if A is real" which your B follows. The eigenvalues of your current matrix B are -1, 1 and 6. So matrix B is not positive definite so it would stand to reason that eigs wouldn't work.

In comparison, the documentation for eig doesn't place any such limitation on the matrix B.

$\endgroup$
1
  • 1
    $\begingroup$ As a quick comment to my answer, if you change the 2's to 1's B becomes positive semi-definite and eigs appears to works, can't guarantee that in general though. $\endgroup$ Sep 30, 2020 at 10:44
1
$\begingroup$

Not really an answer to your question, but the LAPACK implementation of QZ (which solves the GEP) is known to be pretty slow. My PR (https://github.com/Reference-LAPACK/lapack/pull/421) fixes that. You could build my fork from source and use that.

Also watch out when using MKL, they use blocking factor 1 for the HT reduction for some reason, really hampering performance. Depending on the entries, my PR achieves speedups of up to 250, so it might be worth a shot.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.