对病态矩阵进行对角化，无法计算特征向量。使用numpy/scipy得到不同的结果。

Question

对病态矩阵进行对角化，无法计算特征向量。使用numpy/scipy得到不同的结果。

5

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

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]

对于那些感兴趣的人，这些数字代表1D系统的跳跃振幅。它们不是零。哈密顿量由一个稀疏矩阵给出：

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])

- Geralt

1

我必须说C++也能正确计算特征向量。你在C++中使用了哪个线性代数库？ - Warren Weckesser

2

如果您提供了一个最小化、完整和可验证的示例，那么其他人帮助您将更加容易。 - Warren Weckesser

是的，在C++中我使用LAPACK。 - Geralt

2个回答

2

有点令人惊讶的是，它们给出了不同的结果，因为我期望numpy和scipy都会调用LAPACK，但这也基本上是我的经验。

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

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

- Eelco Hoogendoorn

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- francis · Accepted Answer

简短回答：尝试使用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使用分治法计算对称三角矩阵的所有特征值和特征向量。

输入矩阵的小规范并不重要，因为在这种情况下很可能会进行缩放。由于实对称矩阵具有非常不同的特征值（从3.16224604e-63到3.16230928e-03），它是病态的。大多数线性代数过程（包括特征值计算）的准确度都受到这种特性的显著影响。提供的矩阵已经处于三角形式。
np.linalg.eigh()失败了。这可能与使用分而治之的策略有关。
np.linalg.eigvalsh()似乎有效。因此，看起来DSTERF起作用了。然而，这个例程只提供了特征值。

由于LAPACK提供了不同的算法来计算特征值和特征向量，因此你的C ++代码可能使用另一个算法，如dsyev()。在将矩阵缩放并化为三角形形式后，该例程首先调用DORGTR生成正交矩阵，然后调用DSTEQR获取特征向量。希望可以通过scipy的低级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])

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