为什么 SciPy eigsh() 在谐振子的情况下会产生错误的特征值?

计算科学 Python 稀疏矩阵 麻木的 特征值
2021-12-01 23:18:14

在计算物理学的背景下,我正在开发一些更大的代码来执行巨大稀疏矩阵的特征值计算。由于特征值在分析上是众所周知的,因此我在一个维度上针对简单谐振子测试我的例程。这样做并将我自己的例程与 SciPy 的内置求解器进行比较,我发现了下图中显示的奇怪之处。在这里您可以看到前 100 个数值计算的特征值λnum和分析特征值λana

在特征值 40 附近,数值结果开始与分析结果不同。这并不让我感到惊讶(我不会在这里讨论为什么,除非它出现在讨论中)。然而,令我惊讶的是eigsh()会产生退化的特征值(在特征值 80 附近)。为什么即使对于如此少量的特征值,eigsh() 的行为也会如此?

在此处输入图像描述

import numpy as np
from scipy.sparse.linalg import eigsh
import myFunctions as myFunc
import matplotlib.pyplot as plt

#discretize x-axis
N = 100
xmin = -10.
xmax = 10.
accuracy = 1e-5
#stepsize
h = (xmax - xmin) / (N + 1.)
#exclude first and last points since we force wave function to be zero there
x = np.linspace(-10. + h,10. - h,N)
#create potential
V = x**2

def fivePoint(N,h,V):
    C0 = (np.ones(N))*30. / (12. * h * h) + V
    C1 = (np.ones(N)) * (-16.) / (12. * h * h)
    C2 = (np.ones(N)) / (12. * h * h)
    H = sp.spdiags([C2, C1, C0, C1, C2],[-2, -1, 0, 1, 2],N,N)
    return H

H = myFunc.fivePoint(N,h,V)
eigval,eigvec = eigsh(H, k=N-1, which='SM', tol=accuracy)

#comparison analytical and numerical eigenvalues
xAxes = np.linspace(0,len(eigval)-1,len(eigval))
analyticalEigval = 2. * (xAxes + 0.5)
plt.figure()
plt.plot(xAxes,eigval, '+', label=r"$\lambda_{num}$")
plt.plot(xAxes,analyticalEigval, label=r"$\lambda_{ana}$")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.legend(loc=4)
plt.title("eigsh()-method: Comparison of $\lambda_{num}$ and $\lambda_{ana}$")
plt.show()
1个回答

在我看来,某些特征值的退化就像是Lanczos 算法崩溃的标志。Lanczos算法是比较常用的逼近Hermitian矩阵的特征值和特征向量的方法之一;这是 scipy.eigsh() 通过调用ARPACK library使用的。

在精确算术中,Lanczos 算法产生一组正交向量,但在浮点算术中,这些向量可能无法正交,甚至变为线性相关。真正令人讨厌的是,这种正交性的损失恰好发生在一个近似特征值收敛到一个真实特征值时——可以说,算法会破坏自己。结果是你会得到一些虚假的附近特征值对。对此有各种修复方法,例如使用 Gram-Schmidt 强制任何收敛的特征向量在每一步都是正交的。

尽管如此,没有一种方法是完美的,尤其是当您尝试计算矩阵的整个频谱时所以如果你想得到 50 个最小的特征值,你最好用一个有 100 个元素的向量来近似波函数,只要求eigsh()前 50 个能级,而不是使用一个有 50 个点的向量并要求所有的特征值。

如果您想了解更多信息,请查看 Yousef Saad 的大型特征值问题的数值方法