我试图解决一个大矩阵(8000×8000)的特征值问题因为A
和B
在我的例子中是稀疏的,所以我认为使用eigs
比eig
更好但问题是B
是奇异的,而MATLAB的eigs
不能处理它。
有没有解决这个问题的办法?
附笔:A
和B
都是复杂和非对称的
最佳答案
MATLAB使用的底层ARPACK例程允许解决奇异质量矩阵问题(B
到MATLAB,M
到ARPACK)APRACK文档说明:
如果M是奇异的,或者某人对点σ附近的特征值感兴趣,那么用户可以选择使用C=(a-σM)-1M
其中C是移位运算符因此,通过显式指定移位参数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.
...
指定轮班补救措施1
>> 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.
对于示例系统中的大量特征值,所选择的移位产生一个奇异的C矩阵,这是不好的稍微偏离这一点就可以纠正错误。