我在scipy的eigh
函数返回正半定矩阵的负特征值时遇到一些问题。以下是MWE。hess_R
函数返回一个正半定矩阵(它是秩为1的矩阵和对角矩阵的总和,均为非负项)。
import numpy as np
from scipy import linalg as LA
def hess_R(x):
d = len(x)
H = np.ones(d*d).reshape(d,d) / (1 - np.sum(x))**2
H = H + np.diag(1 / (x**2))
return H.astype(np.float64)
x = np.array([ 9.98510710e-02 , 9.00148922e-01 , 4.41547488e-10])
H = hess_R(x)
w,v = LA.eigh(H)
print w
打印的特征值是
[ -6.74055241e-271 4.62855397e+016 5.15260753e+018]
如果我在
np.float64
的return语句中将np.float32
替换为hess_R
,我会得到[ -5.42905303e+10 4.62854925e+16 5.15260506e+18]
相反,所以我猜测这是某种精度问题。
有没有办法解决这个问题?从技术上讲,我不需要使用eigh,但是我认为这是我其他错误的根本问题(取这些矩阵的平方根,获取NaN等)。
最佳答案
我认为问题在于您已经达到浮点精度的极限。线性代数结果的一个很好的经验法则是,对于float32,它们约占10 ^ 8的一部分,对于float64而言,约占10 ^ 16的一部分。这里的特征值小于10 ^ -16。因此,返回值不能真正被信任,并且将取决于您使用的特征值实现的详细信息。
例如,这里有四个您应该可用的求解器。看一下他们的结果:
# using the 64-bit version
for impl in [np.linalg.eig, np.linalg.eigh, LA.eig, LA.eigh]:
w, v = impl(H)
print(np.sort(w))
reconstructed = np.dot(v * w, v.conj().T)
print("Allclose:", np.allclose(reconstructed, H), '\n')
输出:
[ -3.01441754e+02 4.62855397e+16 5.15260753e+18]
Allclose: True
[ 3.66099625e+02 4.62855397e+16 5.15260753e+18]
Allclose: True
[ -3.01441754e+02+0.j 4.62855397e+16+0.j 5.15260753e+18+0.j]
Allclose: True
[ 3.83999999e+02 4.62855397e+16 5.15260753e+18]
Allclose: True
注意,它们都同意较大的两个特征值,但是最小特征值的值随实现而变化。尽管如此,在所有四种情况下,输入矩阵都可以重构为最高64位精度:这意味着算法可以达到预期的精度运行。
关于python - scipy eigh为正半定矩阵给出负特征值,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/36819739/