I have a generalized eigenvalue problem in the standard form
$\lambda \mathbf{B} \mathbf{x} = \mathbf{A} \mathbf{x} $,
resulting from a finite difference discretization of a coupled system of two linear stability equations, so the system is large $(10^5 $ x $10^5)$ and sparse. I know that $\mathbf{B}$ is indefinite, and not symmetric. Since I am interested in the stability of the system, I need to find the first, say ten, eigenvalues with the largest real part.
What it the preferred way to tackle such kind of problem?
I want to solve it using scipy linear algebra package for sparse matrices, but it appears that the eigs function only works when $\mathbf{B}$ is positive definite.
Any help is appreciated.
Edit
I post a minimal example code that generate my issue. The matrices that I use in the example are not the same I have in my application, but they illustrate the problem. A working Matlab snippet, producing the correct answer is:
A = diag([-5, -4, -3, -2, -1]);
B = diag([1, 1, -1, 1, 1]);
eigs(A, B, 2, 'LR')
which produces the correct answer:
>> msolution
ans =
3.0000
-1.0000
An equivalent python version to this problem is:
import numpy as np
from scipy.sparse.linalg import eigs
A = np.diag([-5, -4, -3, -2, -1]).astype(np.float64)
B = np.diag([1, 1, -1, 1, 1]).astype(np.float64)
vals, vecs = eigs(A, 2, B, which='LR')
print vals
producing the clearly wrong result
[ 83.66243085+163.44457559j 83.66243085-163.44457559j]
Now, this might look like a bug in Scipy, but the docs explicitly ask for a positive definite matrix $\mathbf{B}$ and the one I have here is not, as the one I have in my application.