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

Question

我试图解决大矩阵(8000×8000)的特征值问题.由于A和B是在我的情况很少,我认为这是更好地使用eigs比eig.但问题是B单数并且MATLAB eigs无法处理它.

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

PS:这两个A和B是复杂的,非对称的

Answer 1

MATLAB利用的基础ARPACK例程允许解决奇异质量矩阵问题(B到MATLAB,M到ARPACK).APRACK文档说明:

如果 中号是单数,或者如果一个爱好附近的点的特征值σ然后用户可以选择与工作Ç =(甲 - σ 中号)^-1 中号

其中C是移位运算符.因此,通过明确指定移位参数的数值,只要移位不使C奇异sigma,eigs就可以尝试收敛到最接近的支配特征值.

例如,

>> 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)

是的,有一个解决方法.但它需要更多的工作才能稳健而正确地执行,这是norm涉及单个矩阵的问题.

我想说,"最好的方法"是,如果问题适合记忆并且可以在合理的时间内直接计算,那么就这样做.否则,上述移位方法可以付出努力来补充工作.

¹我会注意到并非所有的值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矩阵,这是不好的.稍微偏离这一点可以弥补错误.