python - 病态矩阵的对角化和计算特征向量的不可能性。 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切换到 DSTEDC缩放和减少到三角形后。 DSTEDC 使用分治法计算对称三对角矩阵的所有特征值和特征向量。

1. 输入矩阵的小范数无关紧要，因为在这种情况下可能会应用缩放。由于实对称矩阵具有非常不同幅度的特征值(从 3.16224604e- 63 到 3.16230928e-03), 它是病态的。大多数线性代数程序(包括特征值计算)的准确性都受此属性的显着影响。提供的矩阵已经是三角形式。

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

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

LAPACK 提供 different algorithms to compute the eigenvalues and eigenvectors ，您的 C++ 代码可能使用另一个代码，例如 dsyev() . 将矩阵缩放并简化为三角形式后，此例程首先调用DORGTR 生成正交矩阵，然后调用DSTEQR 获取特征向量。希望这个函数可以通过 scipy 的 Low-level LAPACK functions (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])

您还可以避免还原为三角形，因为您的矩阵已经是三角形了。可能需要缩放以避免数值问题。然后通过LAPACK functions for Cython可以直接调用DORGTR和DSTEQR .但它需要更多的工作和一个编译步骤。