病态矩阵的对角化和不可能计算特征向量。 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的源代码中。的确，如果不需要特征向量，例程将矩阵缩放，将矩阵减少到Trigonal形式（DSYTRD），最后调用DSTERF 。最后一个例程使用QL或QR算法的Pal-Walker-Kahan变体计算对称三角形矩阵的所有特征值。相反，如果需要特征向量，DSYEVD 987654326在缩放和减少到Trigonal形状后切换到@ 987654326。 DSTEDC计算使用鸿沟和征服方法计算对称三角形矩阵的所有特征值和特征向量。

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

1. np.linalg.eigh()失败。它可能与使用分裂和征服策略有关。

2. np.linalg.eigvalsh()似乎工作。因此，它看起来像DSTERF工作。尽管如此，该常规仅提供特征值。

作为 lapack提供different algorithms to compute the eigenvalues and eigenvectors，您的C ++代码可能使用另一个，例如dsyev()。在缩放和减少Trigonal形式之后，这个例程首先调用DORGTR要生成正交矩阵，然后调用DSTEQR以获取特征向量。希望，可以通过SCIPY Low-level LAPACK functions (scipy.linalg.lapack)

调用此功能

我向代码添加了几行以调用此函数。 scipy.linalg.lapack.dsyev()计算与此Ill Tranningned矩阵类似于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])

您还可以避免减少到Trigonal表单，因为您的矩阵已经是三角形。可能需要缩放来避免数值问题。然后DORGTR和DSTEQR可以直接通过LAPACK functions for Cython。但它需要更多的工作和汇编步骤。

Answer 2

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

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

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