病态矩阵的对角化和计算特征向量的不可能性。 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++ 也能够正确计算特征向量。

所以我有两个问题：

为什么函数 np.linalg.eigh(denseHam) 给出的光谱与 np.linalg.eigvalsh(denseHam) 不同？
有什么方法可以用 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

它们给出不同的结果有点令人惊讶，因为我希望 numpy 和 scipy 都调用 LAPACK，但这通常也是我的经验。

请注意，scipy 绑定提供了更多可玩的参数； numpy 可能使用不同的默认值。似乎需要进行一些实验；你的问题不只是有非常小的元素（如果它导致下溢，可以通过简单的缩放来解决），但你的问题也非常 'stiff'，特征值跨越 70 个数量级。 C++ 可能会给你特征向量，但如果它们被数字噪声污染到无用的地步，我不会感到惊讶。

这听起来像是那种用某种 transformed/preconditioned space 来解决它会更可靠的问题。文档字符串没有说明 LAPACK 函数是否可以处理 128 位浮点数；上次我试过他们没有，但如果他们现在这样做，请确保至少使用它而不是 64 位。

Answer 2

简答：试试 LAPACK 的 scipy.linalg.lapack.dsyev() !

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

函数 np.linalg.eigvalsh() and np.linalg.eigh() both call LAPCK's DSYEVD as stated in their documentations. Nevertheless, they deliver different eigenvalues: eigh() fails. The reason is likely engraved in the source code of DSYEVD. Indeed, if the eigenvectors are not needed, the routine scales the matrix, reduces the matrix to trigonal form (DSYTRD) and finally calls DSTERF. This last routine computes all eigenvalues of a symmetric tridiagonal matrix using the Pal-Walker-Kahan variant of the QL or QR algorithm. On the contrary, if eigenvectors are needed, DSYEVD switches to DSTEDC 缩放和缩减为三角形后。 DSTEDC 使用分治法计算对称三对角矩阵的所有特征值和特征向量。

输入矩阵的小范数无关紧要，因为在这种情况下可能会应用缩放。因为实对称矩阵的特征值非常不同的大小（从 3.16224604e-63 到 3.16230928e-03），它是病态的。大多数线性代数程序（包括特征值计算）的准确性受此影响很大属性。 提供的矩阵已经是三角形式。
np.linalg.eigh() 失败。这可能与使用分而治之策略有关。
np.linalg.eigvalsh() seems to work. Therefore, it looks like DSTERF 有效。然而，此例程仅提供特征值。

因为 LAPACK 提供 different algorithms to compute the eigenvalues and eigenvectors, your C++ code likely uses another one, such as dsyev(). After scaling and reducing the matrix to trigonal form, this routine first call DORGTR to generate the orthogonal matrix, then call DSTEQR to get the eigenvectors. Hopefully, this function can be called through scipy's 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])

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

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

Diagonalization of ill-conditioned matrix and imposibility to compute eigenvectors. Different results with numpy/scipy

python

numpy

linear-algebra

scipy