编辑

我遇到了同样的问题，并与我的主管交谈。这是他的解决方案：

使用 -B 代替 B：

在我的例子中，矩阵 B 是负定的（所有特征值都是负的）。

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

这将给出相同的结果，但 -B 是肯定的，您可以使用 scipy.sparse.linalg.eigs()。通过尝试此解决方案，某些矩阵大小出现了 ARPACK 错误：

scipy.sparse.linalg.eigen.arpack.arpack.ArpackError：ARPACK 错误 -9999：无法构建 Arnoldi 分解。IPARAM(5) 返回当前 Arnoldi 分解的大小。建议用户检查是否已分配足够的工作空间和阵列存储

我不知道该错误是否与我的特定配置和代码或其他内容有关。

LinearOperator 的其他替代方案：

在我的情况下，我不必矩阵 A 和 B，而只需要关联的 LinearOperators。我为 LinearOperators 升级了 sebas 的解决方案，它运行良好，但它比第一个解决方案慢，因为它需要反转 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)

这很奇怪。它应该是 scipy.sparse.linalg.eigs 中的一个错误。

查看 matlab，我们看到open eigs在第 822-825 行打开函数，他们说

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

我一直在 scipy 中寻找相同的代码，似乎他们在scipy/sparse/linalg/eigen/arpack/arpack.py#L541的第 541-543 行做同样的事情

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

其中 OPb 和 OPa 之前已被定义为乘积$A$和倒数$B$. 所以我没有发现差异，所以它应该是一个错误。

在 Matlab 和 Scipy 中，他们似乎正在解决$B^{-1}Ax = \lambda x$. 正因为如此，我想知道为什么$B$需要是正定的，而不仅仅是非奇异的（我认为$B$对于移位和反转算法必须是正定的，因为它用于定义范数）。因此，你可以简单地修改你的代码来做同样的事情：

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

获得正确的解决方案

/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]

我收到警告是因为$B$在通过调用 进行稀疏 LU 分解（splu）时不是 csr 格式factorized，因此因式分解函数执行到 csr 的转换并进行稀疏 LU。我猜你的矩阵已经是 csr 格式了，所以不会有问题。

在这里，您还可以使用不完整的 LU 分解，解决系统$B$迭代地，并且有很多用于重新排序的选择$B$在 LU 分解之前使其高效。