scipy.linalg.sparse.eigsh returns 半正定矩阵的负特征值和非一致特征值

Question

我尝试使用scipy.linalg.sparse.eigsh（我们称之为方法1：M1）来计算实对称半定矩阵的拉普拉斯矩阵的最小特征值W.

作为基准，我运行针对 scipy.linalg.eigh 的计算（方法 2：M2）。似乎M1returns与M2的特征值不同，而且那些特征值似乎是错误的。

代码如下：

# Initialization
N=10
np.random.seed(0) # for reproducibility 
Nvp=4

# Create a real symmetric semi-definite matrix W, precision float64
X=np.random.random((N,N))
W=np.dot(X, X.T) 
W=np.array(W, dtype=np.float64) 

# Compute its Laplacian matrix A, precision float64
d=np.array([sum(W[:][j]) for j in range(N)], dtype=np.float64)
D=np.diag(d)
A=D-W
for i in range(N):
    for j in range(N):
        A[i][j]/=np.sqrt(d[i], dtype=np.float64)*np.sqrt(d[j], dtype=np.float64)
A=np.array(A, dtype=np.float64)

让我们检查一下 A 的格式是否正确：

>>> A.dtype
dtype('float64')
>>> np.allclose(A, A.T)
True

现在让我们运行进行一些测试：

## 1) Compute A's smallest eigenvalues by 2 different means
wA2, vA2 = la.eigh(A)
wA1, vA1 = sparse.eigsh(A, k=Nvp, sigma=0, which='LM')

## 2) Compute W's smallest eigenvalues by 2 different means
wW2, vW2 = la.eigh(W)
wW1, vW1 = sparse.eigsh(W, k=Nvp, sigma=0, which='LM')

# Output computed eigenvalues
print(wA2[:Nvp])
print(wA1[:Nvp])
print(wW2[:Nvp])
print(wW1[:Nvp])

这是输出：

>>>[-1.88737914e-15  9.03999970e-01  9.23513143e-01  9.52678423e-01]

[-4.93242313e-01 -8.14381749e-17  9.22235466e-01  9.44848237e-01]

[0.00575077 0.04770667 0.08565863 0.16319798]

[0.00575077 0.04770667 0.08565863 0.16319798]

第一个输出显示 M1 计算了 A 的负特征值，这在数学上是不可能的。此外，如果一个检查其他人的计算，让我们说第三个：

>>> np.dot(A, vA2[:,2])-wA2[2]*vA2[:,2]
array([-0.01183104, -0.25059123,  0.47372558, -0.31358618, -0.2968036 ,
   -0.04832199,  0.40973325, -0.01369472,  0.33267859, -0.27122678])

它甚至不接近于零。我必须添加每次不同的计算特征值。对我来说，这是由于初始化向量。我会说 scipy.linalg.sparse.eigsh 的迭代次数不足以接近真实结果，但设置 maxiter=1000000 不会影响 st运行ge 结果.关于负特征值，不幸的是我一无所知。

我是运行宁 :

Python 3.7.0（默认，2018 年 6 月 28 日，13:15:42）

[GCC 7.2.0] :: Anaconda, Inc. linux

NumPy 和 SciPy 是针对 Intel MKL 构建的。

有谁能赐教吗？预先感谢您的宝贵时间。

Answer 1

您的矩阵 A 是奇异矩阵，即 0 是 A 的特征值。 shift-invert 方法中 sigma 的值本身不能是特征值。查看 ARPACK 文档的 "Shift and Invert Spectral Transformation Mode" 部分。

相同的公式显示在SciPy tutorial中。 shift-invert 方法的公式是

inv(A - sigma*M) @ M @ x = nu*x

(Argh，我真希望 Whosebug 有 LaTeX 标记！）当 sigma 为 0 时，A - sigma*M 就是 A，并且inv(A) 不存在。

scipy.linalg.sparse.eigsh returns 半正定矩阵的负特征值和非一致特征值

scipy.linalg.sparse.eigsh returns negative and non-consistent eigenvalues for positive semi-definite matrix

python

numpy

scipy

eigenvalue

python-3.x