Solving a generalised eigenvalue problem

Question

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.

Answer 1

I had the same problem and spoke with my supervisor. Here is his solution:

Use -B instead of B:

In my case, matrix B was negative definite (all eigenvalues are negative).

λ(-B)x=(-A)x

which will give the same result, but -B is positive definite and you can use scipy.sparse.linalg.eigs(). By trying this solution, an error in ARPACK arose for some matrix sizes :

scipy.sparse.linalg.eigen.arpack.arpack.ArpackError: ARPACK error -9999: Could not build an Arnoldi factorization. IPARAM(5) returns the size of the current Arnoldi factorization. The user is advised to check that enough workspace and array storage has been allocated

I dont know if that error is bound to my particular configuration and code or something else.

Other alternative for LinearOperator:

In my cases, I didnt have have to matrices A and B but only to the associate LinearOperators. I upgraded the solution of sebas for LinearOperators and it worked well, but it is slower than the first solution as it needs to invert B.

from scipy.sparse.linalg import LinearOperator

N = your_size

def A_fct(v):
    # write the function you need or just define matrix A
    return np.dot(A, v)
linop_A = LinearOperator( (N,N), matvec = A_fct)

def B_fct(v):
    # write the function you need or just define matrix B
    return np.dot(B, v)
linop_B = LinearOperator( (N, N), matvec = B_fct)

def Binv_A(v):
    retour = linop_A.matvec(v)
    sol, info = scipy.sparse.linalg.gmres(linop_B, retour)
    # or similar to gmres, it calculate inv(B) * A * v
    return sol
linop_Binv_A = LinearOperator( (N, N), matvec = Binv_A)

Answer 2

It si very weird. It should be a bug in scipy.sparse.linalg.eigs.

Looking at matlab, we see that opening the function with open eigs, at line 822-825, they say that

% eigs(A,B,k,stringSigma), stringSigma~='SM'
% A*x = lambda*B*x
% => OP = inv(B)*A and use standard inner product.
mode = 2;

I have been looking for the same piece of code inside scipy, and it seems that they are doing the same at lines 541-543 in scipy/sparse/linalg/eigen/arpack/arpack.py#L541

elif self.mode == 2:
    self.workd[xslice] = self.OPb(self.workd[xslice])
    self.workd[yslice] = self.OPa(self.workd[xslice])

where the OPb and OPa have been defined before as the product by $A$ and the inverse of $B$. So I did not found the difference, so it should be a bug.

In Matlab and also in Scipy, they seem to be solving $B^{-1}Ax = \lambda x$. Because of that I was wondering why $B$ is required to be positive definite and not just nonsingular (I think that $B$ must be positive definite for the shift-and-invert algorithm, because there it is used to define a norm). Thus, you can simply modify your code to do the same:

import numpy as np
from scipy.sparse.linalg import eigs
from scipy.sparse.linalg import LinearOperator
from scipy.sparse.linalg import factorized

A = np.diag([-5, -4, -3, -2, -1]).astype(np.float64)
B = np.diag([1, 1, -1, 1, 1]).astype(np.float64)

B_inv = factorized(B);
def mv(v):
    temp = np.dot(A, v)
    return B_inv(temp)

BA = LinearOperator(A.shape, matvec = mv) # BA.matvec(x) is np.dot(inv(B),np.dot(A,x))
vals, vecs = eigs(BA, 2, which='LR')
print vals

obtaining the right solution

/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py:258: SparseEfficiencyWarning: splu requires CSC matrix format
  warn('splu requires CSC matrix format', SparseEfficiencyWarning)
[ 3.+0.j -1.+0.j]

I got a warning because $B$ is not in csr format when performing the sparse LU decomposition (splu) by means of the call to factorized, so the factorization function perform the conversion to csr and do the sparse LU. I guess that your matrices already are in csr format, so it will not be a problem.

Here you can also play with an incomplete LU factorization, solving the system with $B$ iteratively, and lots of choices for the reordering of $B$ before the LU factorization to make it efficient.