为何Scipy中eigs无法复现eig的广义特征值问题计算结果？

问题背景 #

我遇到一个复方阵的广义特征值问题：

$Aq = \alpha Bq$
使用scipy.linalg.eig求解能得到正确结果，特征值$\alpha$在复平面的分布符合预期；但用Scipy稀疏求解器eigs计算时结果异常，代码如下：

from scipy.sparse.linalg import eigs
alpha, q = eigs(A, k=40, B=B, sigma=1.0+0.0J)

原本eigs在稀疏矩阵场景下应更高效，目前猜测异常与A、B的矩阵特性相关，但不确定具体原因。

本地加载矩阵及结果的代码：

import numpy as np

# 加载原始矩阵
data = np.load('YourPath/matrix.npz')
A = data['A']
B = data['B']

# 加载scipy.linalg.eig的正确结果
data = np.load('YourPath/eig_result.npz')
alpha_eig = data['alpha']
q_eig = data['q']

# 加载scipy.sparse.linalg.eigs的异常结果
data = np.load('YourPath/eigs_result.npz')
alpha_eigs = data['alpha']
q_eigs = data['q']

可能的原因及解决方案 #

• B矩阵的正定性/奇异性问题：eigs在处理广义特征值问题且指定sigma时，默认依赖B矩阵可逆且正定（Hermite正定，即B为Hermite矩阵且特征值全为正）。若B奇异、非正定或非Hermite，会导致迭代过程数值不稳定，输出异常结果。

  • 验证方式：计算B的特征值，检查是否存在零特征值（奇异），或特征值实部不全为正（非正定）；同时检查B是否满足Hermite性（np.allclose(B, B.conj().T)）。
  • 解决方向：若B可逆但非正定，可尝试改用eigs(A, k=40, M=B, which='LM')（指定M为B，求解模最大的特征值）；若B奇异，需先处理矩阵奇异性。

• sigma参数与k值不匹配：sigma=1.0+0.0J指定搜索靠近1的特征值，但如果实际靠近1的特征值数量不足40，eigs会返回数值误差较大的结果，甚至引入虚假特征值。

  • 验证方式：统计scipy.linalg.eig结果中，与1的距离小于阈值（如1e-2）的特征值数量。
  • 解决方向：调整k为实际符合条件的特征值数量，或调整sigma到特征值更集中的区域。

• 迭代收敛性不足：eigs基于Arnoldi/Lanczos迭代，若默认迭代次数或收敛阈值不满足当前矩阵的收敛要求，会导致结果不收敛。

  • 解决方向：增加迭代次数、收紧收敛阈值：

    alpha, q = eigs(A, k=40, B=B, sigma=1.0+0.0J, maxiter=1000, tol=1e-8)
    

• 矩阵稀疏度不足：若A、B实际为稠密矩阵（非零元素占比高），用稀疏求解器eigs反而会引入额外的数值误差，且效率不如稠密求解器scipy.linalg.eig。

  • 验证方式：计算矩阵的稀疏度（A.nnz / (A.shape[0]*A.shape[1])），若稀疏度高于50%，建议优先使用稠密求解器。

验证步骤 #

1. 对比特征值分布：将eig和eigs的特征值绘制在复平面，定位异常结果的位置，判断是否为偏离sigma的无效特征值。
2. 测试不同搜索策略：改用which='LR'（实部最大）、which='LI'（虚部最大）等参数，查看是否能得到与eig一致的特征值，验证求解器稳定性。
3. 检查残差：对eigs得到的特征值和特征向量，计算残差$||Aq - \alpha Bq||$，若残差远大于eig的残差，说明结果不收敛或存在数值问题。

内容的提问来源于stack exchange，提问作者Adoledin