病态矩阵的对角化和无法计算特征向量。numpy/scipy 的不同结果

Question

我正在处理一个包含非常小的元素的稀疏矩阵。考虑一个向量：

vec=[-1.e-76 -1.e-72 -1.e-68 -1.e-64 -1.e-60 -1.e-56 -1.e-52 -1.e-48 -1.e-44
-1.e-40 -1.e-36 -1.e-32 -1.e-28 -1.e-24 -1.e-20 -1.e-16 -1.e-12 -1.e-08
-1.e-04 -1.e-02 -1.e-04 -1.e-08 -1.e-12 -1.e-16 -1.e-20 -1.e-24 -1.e-28
-1.e-32 -1.e-36 -1.e-40 -1.e-44 -1.e-48 -1.e-52 -1.e-56 -1.e-60 -1.e-64
-1.e-68 -1.e-72 -1.e-76]

对于那些感兴趣的人来说，这些数字代表一维系统的跳跃幅度。它们不为零。哈密​​顿量由稀疏矩阵给出：

H=sps.diags([vec,vec],[-1,1],dtype='f8')

我对特征值感兴趣，但对特征向量更感兴趣。据我所知，处理对角化有两种方法： scipy.linalg和 numpy.linalg，前者更好。

 denseHam=H.toarray()

正确的特征值谱由所有这些函数给出：

import numpy as np
import scipy.linalg as la
s1= la.eigvalsh(denseHam)
s2= np.linalg.eigvalsh(denseHam)
s3= np.linalg.eigvals(denseHam) #I did not expect that!

正确的频谱是：

spectrum=[-3.16230928e-03 -3.16227766e-08 -3.16227766e-13 -3.16227766e-18
-3.16227766e-23 -3.16227766e-28 -3.16227766e-33 -3.16227766e-38
-3.16227766e-43 -3.16227766e-48 -3.16227766e-53 -3.16227766e-58
-3.16224604e-63  3.16224604e-63  3.16227766e-58  3.16227766e-53
 3.16227766e-48  3.16227766e-43  3.16227766e-38  3.16227766e-33
 3.16227766e-28  3.16227766e-23  3.16227766e-18  3.16227766e-13
 3.16227766e-08  3.16230928e-03]

然而，其他函数（也涉及特征向量的计算）失败，我无法继续，因为我需要特征向量。

我不得不说 C++ 也能够正确计算特征向量。

所以我有两个问题：

1. 为什么该函数np.linalg.eigh(denseHam)给出的频谱与 不同np.linalg.eigvalsh(denseHam)？
2. 有没有办法用python正确计算特征向量？

预先非常感谢您！

--- 更新------ 我在这里粘贴一个最小的完整示例。注意 的暴露简并性numpy.linalg.eigh：

import numpy as np
import scipy.sparse as sps

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

s1=np.linalg.eigvalsh(denseHam)
(s2,basis)=np.linalg.eigh(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem])

Answer 1

简短的回答：尝试 LAPACK 的scipy.linalg.lapack.dsyev()！

# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

函数np.linalg.eigvalsh()和np.linalg.eigh()两者都调用 LAPCK，DSYEVD如其文档中所述。然而，它们提供了不同的特征值：eigh()失败。原因很可能刻在 的源代码中DSYEVD。事实上，如果不需要特征向量，例程会缩放矩阵，将矩阵简化为三角形式 ( DSYTRD)，最后调用DSTERF。最后一个例程使用 QL 或 QR 算法的 Pal-Walker-Kahan 变体计算对称三对角矩阵的所有特征值。相反，如果需要特征向量，则在缩放和缩减为三角形之后DSYEVD切换。使用分治法计算对称三对角矩阵的所有特征值和特征向量。DSTEDCDSTEDC

1. 输入矩阵的小范数并不重要，因为在这种情况下可能会应用缩放。由于实对称矩阵具有大小差异很大的特征值（从 3.16224604e-63 到 3.16230928e-03），因此它是病态的。大多数线性代数过程（包括特征值计算）的准确性都会受到此属性的显着影响。提供的矩阵已经是三角形式。

2. np.linalg.eigh()失败了。这可能与分而治之策略的使用有关。

3. np.linalg.eigvalsh()似乎有效。因此，它看起来DSTERF有效。然而，该例程仅提供特征值。

由于LAPACK 提供了不同的算法来计算特征值和特征向量，因此您的 C++ 代码可能使用另一种算法，例如dsyev()。将矩阵缩放并简化为三角形式后，该例程首先调用DORGTR生成正交矩阵，然后调用DSTEQR获取特征向量。希望这个函数可以通过 scipy 的Low-level LAPACK 函数来调用( )scipy.linalg.lapack

我在您的代码中添加了几行来调用此函数。计算与此病态矩阵scipy.linalg.lapack.dsyev()类似的特征值。np.linalg.eigvalsh()

import numpy as np
import scipy.sparse as sps
import scipy.linalg.lapack

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

# eigenvalues only, boils down to LAPACK's DSTERF
s1=np.linalg.eigvalsh(denseHam)
# LAPACK's dsyevd, divide and conquer
(s2,basis)=np.linalg.eigh(denseHam)
# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem], s3[elem])

您还可以避免简化为三角形形式，因为您的矩阵已经是三角形的。可能需要进行缩放以避免数值问题。然后可以通过Cython 的 LAPACK 函数直接调用DORGTR和。但它需要更多的工作和编译步骤。DSTEQR