MATLAB - Eigs(A,B) 当 B 为单数时

Question

我正在尝试解决大型矩阵 (8000×8000) 的特征值问题。由于 A 和 B 在我的例子中是稀疏的，我认为使用 eigs 比 eig 更好。但问题是B是单数的，MATLAB的eigs无法处理。

是否有任何变通方法可以解决此问题？

附言:A 和B 都是复数和非对称的

Answer 1

MATLAB 利用的底层 ARPACK 例程允许解决奇异质量矩阵问题(B 到 MATLAB，M 到 ARPACK)。 APRACK 文档指出:

If M is singular or if one is interested in eigenvalues near a point σ then a user may choose to work with C = (A - σ M)^-1M

其中 C 是移位运算符。因此，通过为 shift 参数明确指定一个数值 sigma , eigs只要移位不使 C 奇异，就可以尝试收敛到最近的主导特征值。

例如，

>> n = 1000;
>> A = sprandn(n,n,1/n) + speye(n,n); % ensure invertibility of A
>> B = sprandn(n,n,1/n); B(1:n-1,1:n-1) = B(1:n-1,1:n-1) + speye(n-1); % mildly pathological B
>> condest(B)
ans =
   Inf

A 不需要是可逆的，但它会让我们稍后检查我们的答案。在没有指定类次的情况下调用 eigs 会引发错误:

>> eigs(A,B,1)
Error using eigs/checkInputs/LUfactorB (line 954)
B is singular.
...

指定类次可解决此问题¹

>> eigs(A,B,1,0)
ans =
    0.1277

但这只是在 0 附近发现的主要特征值。让我们检查一些其他点

>> sort(arrayfun(@(sigma)eigs(A,B,1,sigma),linspace(-10,10,10).'),'descend')
ans =
   4.1307 + 0.2335i
   4.1307 + 0.2335i
   4.1307 + 0.2335i
   3.3349 + 0.0000i
   1.1267 + 0.0000i
   0.1277 + 0.0000i
   0.1277 + 0.0000i
   0.1277 + 0.0000i
   0.1277 + 0.0000i
   0.1277 + 0.0000i

如果看起来像 4.1307 + 0.2335i 可能是系统的主导特征值。让我们围绕该点寻找更多特征值:

>> eigs(A,B,5,4.13)
ans =
   4.1307 - 0.2335i
   4.1307 + 0.2335i
   3.3349 + 0.0000i
   2.8805 + 0.0000i
   2.6613 + 0.0000i

看起来不错。但是有更大的有限特征值吗？幸运的是，由于 A 是可逆的，我们可以通过取 eig(B/A) 的倒数直接检查特征值:

>> lam = sort(1./eig(full(B/A)),'descend')
lam =
      Inf + 0.0000i
   4.1307 + 0.2335i
   4.1307 - 0.2335i
   3.4829 + 1.6481i
   3.4829 - 1.6481i
   3.3349 + 0.0000i
   2.4162 + 2.1442i
   2.4162 - 2.1442i
   2.8805 + 0.0000i
   2.2371 + 1.7137i
   2.2371 - 1.7137i
   ...

首先，我们看到烦人的无限特征值给出了所有问题。其次，我们可以看到 eigs 确实找到了最大的有限特征值，但没有找到复平面中幅度稍小的特征值，因为纯实数特征值更接近于转移点:

>> [lam,abs(lam-4.13)]
ans =
      Inf + 0.0000i      Inf + 0.0000i
   4.1307 + 0.2335i   0.2335 + 0.0000i % found by eigs(A,B,5,4.13)
   4.1307 - 0.2335i   0.2335 + 0.0000i % found by eigs(A,B,5,4.13)
   3.4829 + 1.6481i   1.7705 + 0.0000i
   3.4829 - 1.6481i   1.7705 + 0.0000i
   3.3349 + 0.0000i   0.7951 + 0.0000i % found by eigs(A,B,5,4.13)
   2.4162 + 2.1442i   2.7450 + 0.0000i
   2.4162 - 2.1442i   2.7450 + 0.0000i
   2.8805 + 0.0000i   1.2495 + 0.0000i % found by eigs(A,B,5,4.13)
   (8 more complex eigenvalues)
   2.6613 + 0.0000i   1.4687 + 0.0000i % found by eigs(A,B,5,4.13)

所以是的，有一个解决方法。但它需要更多的工作才能稳健和正确地执行，这是涉及奇异矩阵的问题的标准。

“最好的方法”，我想说的是，如果问题适合内存并且可以在合理的时间内直接计算，那么就这样做。否则，上述轮类方法可以补充工作量。

---

¹ 我会注意到并非所有 sigma 值都有效。对于上面给出的示例，1 失败:

>> eigs(A,B,5,1)
Error using eigs/checkInputs/LUfactorAminusSigmaB (line 994)
The shifted operator is singular. The shift is an eigenvalue.
 Try to use some other shift please.

对于示例系统中的大量1 特征值，所选择的移位会创建奇异的C 矩阵，这并不好。稍微偏离这一点可以纠正错误。