使用ARPACK解决eigenvalueproblem，但得到用Matlab不一致的结果(Using

我是新来ARPACK，我下载了类似下面的脚本

import time
import numpy as np
from scipy.linalg import eigh
from scipy.sparse.linalg import eigs
np.set_printoptions(suppress=True)

n=30
rstart=0
rend=n

A=np.zeros(shape=(n,n))

# first row
if rstart == 0:
    A[0, :2] = [2, -1]
    rstart += 1
# last row
if rend == n:
    A[n-1, -2:] = [-1, 2]
    rend -= 1
# other rows
for i in range(rstart, rend):
    A[i, i-1:i+2] = [-1, 2, -1]

A[0,8]=30

start_time = time.time()

evals_large, evecs_large = eigs(A, 10, sigma=3.6766133, which='LM')
print evals_large

end_time=time.time()-start_time
print(" Elapsed time: %12f seconds " % end_time)

它解决了一个非常简单的特征值问题（矩阵A没有对称的，我设置A[0,8]是30 ）。最靠近3个特征值3.6766133 （ sigma=3.6766133在设定）根据ARPACK结果

[ 3.68402411+0.j  3.82005897+0.j  3.51120293+0.j]

然后我去MATLAB和解决同特征值问题，其结果是

4.144524409923138 + 0.000000000000000i 
3.642801014184622 + 0.497479798520641i
3.642801014184622 - 0.497479798520641i
2.372392770347609 + 0.762183281789166i
2.372392770347609 - 0.762183281789166i
3.979221766266502 + 0.000000000000000i
3.918541441830947 + 0.000000000000000i
3.820058967057387 + 0.000000000000000i 
3.684024113506185 + 0.000000000000000i
3.511202932803536 + 0.000000000000000i
3.307439963195127 + 0.000000000000000i
3.080265978640102 + 0.000000000000000i
2.832849552917550 + 0.000000000000000i
2.565972630556613 + 0.000000000000000i
2.283744793210587 + 0.000000000000000i
1.996972474451519 + 0.000000000000000i
0.927737801889518 + 0.670252740725955i
0.927737801889518 - 0.670252740725955i
1.714561796881689 + 0.000000000000000i
-0.015193770830045 + 0.264703483268519i
-0.015193770830045 - 0.264703483268519i
1.438919271663752 + 0.000000000000000i
0.019951101383019 + 0.000000000000000i
0.080534338862828 + 0.000000000000000i 
0.181591307101504 + 0.000000000000000i
0.318955140475174 + 0.000000000000000i
0.488231021129767 + 0.000000000000000i
0.688030188040126 + 0.000000000000000i
1.171318650526539 + 0.000000000000000i
0.917612528393044 + 0.000000000000000i

显然，第二种模式3.642801014184622 + 0.497479798520641i更接近sigma=3.6766133 ，但ARPACK没有把它挑出来。

可能是什么问题呢？你能不能帮我想出解决办法？非常感谢。

Answer 1:

首先是关于MATLAB的功能有几件事情：

通过返回特征值eig是不排序。在[V,D] = eig(A)我们只保证的列V是在在相应的右本征向量的本征值D(i,i) 。在另一方面， svd返回递减顺序排列的奇异值。
d = eigs(A,k)返回k 最大幅度特征值。但是它是用于大型稀疏矩阵，通常是不能代替：
```
 d = eig(full(A)); d = sort(d, 'descend'); d = d(1:k); 
```
（ eigs基于ARPACK ，而eig使用LAPACK例程）。
没有自然排序复数。的约定是，所述sort功能首先通过大小排序复元素（即abs(x)上），然后通过相位角[-pi,pi]间隔（即angle(x)如果幅度是相等的。

MATLAB

考虑到这一点，考虑下面的MATLAB代码：

% create the same banded matrix you're using
n = 30;
A = spdiags(ones(n,1)*[-1,2,-1], [-1 0 1], n, n);
A(1,9) = 30;
%A = full(A);

% k eigenvalues closest to sigma
k = 10; sigma = 3.6766133;
D = eigs(A, k, sigma);

% lets check they are indeed sorted by distance to sigma
dist = abs(D-sigma);
issorted(dist)

我得到：

>> D
D =
  3.684024113506185 + 0.000000000000000i
  3.820058967057386 + 0.000000000000000i
  3.511202932803535 + 0.000000000000000i
  3.918541441830945 + 0.000000000000000i
  3.979221766266508 + 0.000000000000000i
  3.307439963195125 + 0.000000000000000i
  4.144524409923134 + 0.000000000000000i
  3.642801014184618 + 0.497479798520640i
  3.642801014184618 - 0.497479798520640i
  3.080265978640096 + 0.000000000000000i

>> dist
dist =
   0.007410813506185
   0.143445667057386
   0.165410367196465
   0.241928141830945
   0.302608466266508
   0.369173336804875
   0.467911109923134
   0.498627536953383
   0.498627536953383
   0.596347321359904

您可以尝试使用得到密集类似的结果eig ：

% closest k eigenvalues to sigma
ev = eig(full(A));
[~,idx] = sort(ev - sigma);
ev = ev(idx(1:k))

% compare against eigs
norm(D - ev)

所不同的是可以接受的小（接近机器精度）：

>> norm(ev-D)
ans =
     1.257079405021441e-14

蟒蛇

同样在Python：

import numpy as np
from scipy.sparse import spdiags
from scipy.sparse.linalg import eigs

# create banded matrix
n = 30
A = spdiags((np.ones((n,1))*[-1,2,-1]).T, [-1,0,1], n, n).todense()
A[0,8] = 30

# EIGS: k closest eigenvalues to sigma
k = 10
sigma = 3.6766133
D = eigs(A, k, sigma=sigma, which='LM', return_eigenvectors=False)
D = D[::-1]
for x in D:
    print '{:.16f}'.format(x)

# EIG
ev,_ = np.linalg.eig(A)
idx = np.argsort(np.abs(ev - sigma))
ev = ev[idx[:k]]
for x in ev:
    print '{:.16f}'.format(x)

类似的结果：

# EIGS
3.6840241135061853+0.0000000000000000j
3.8200589670573866+0.0000000000000000j
3.5112029328035343+0.0000000000000000j
3.9185414418309441+0.0000000000000000j
3.9792217662665070+0.0000000000000000j
3.3074399631951246+0.0000000000000000j
4.1445244099231351+0.0000000000000000j
3.6428010141846170+0.4974797985206380j
3.6428010141846170-0.4974797985206380j
3.0802659786400950+0.0000000000000000j

# EIG
3.6840241135061880+0.0000000000000000j
3.8200589670573906+0.0000000000000000j
3.5112029328035339+0.0000000000000000j
3.9185414418309468+0.0000000000000000j
3.9792217662665008+0.0000000000000000j
3.3074399631951201+0.0000000000000000j
4.1445244099231271+0.0000000000000000j
3.6428010141846201+0.4974797985206384j
3.6428010141846201-0.4974797985206384j
3.0802659786400906+0.0000000000000000j

Answer 2:

如果您使用的结果是与NumPy与Matlab一致eigs功能：

>> format long
>> A = diag(-ones(n-1,1),-1) + diag(2*ones(n,1)) + diag(-ones(n-1,1),+1);A(1,9)=30;
>> eigs(A,3,3.6766133)'
ans =
   3.684024113506185   3.820058967057386   3.511202932803534

至于为什么没有选择真正最接近的特征值，我觉得有收敛实矩阵的特征值的复杂和一个真正的转变的选择做。我不知道怎么ARPACK计算其迭代，但我记得被告知，有一个真正的σ真正的A不能在默认情况下收敛到一个复杂的共轭对，因为他们的绝对值比1（逆幂迭代）。由于ARPACK将产生在8 ^日和9 ^日反复复特征值（+/-顺序是随机的），我猜他们是一些定为这，我已经忘记或者从来不知道：

>> ev = eigs(A,9,3.6766133);ev(8:9)
ans =
  3.642801014184617 - 0.497479798520639i
  3.642801014184617 + 0.497479798520639i

我不知道如果他们的是比猜测的移的复部分或仅仅霎那额外的特征值，直到共轭对落入收敛球为ARPACK方法一般工作围绕这个其他。

文章来源: Using ARPACK solving eigenvalueproblem, but getting inconsistent results with Matlab