scipy eigh为正半定矩阵给出负特征值 [英] scipy eigh gives negative eigenvalues for positive semidefinite matrix
问题描述
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]
如果在hess_R的return语句中将np.float64替换为np.float32，我会得到
[ -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位精度:这意味着算法可以达到预期的精度运行.
I am having some issues with scipy's eigh function returning negative eigenvalues for positive semidefinite matrices. Below is a MWE.
The hess_R function returns a positive semidefinite matrix (it is the sum of a rank one matrix and a diagonal matrix, both with nonnegative entries).
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
The eigenvalues printed are
[ -6.74055241e-271 4.62855397e+016 5.15260753e+018]
If I replace np.float64 with np.float32 in the return statement of hess_R I get
[ -5.42905303e+10 4.62854925e+16 5.15260506e+18]
instead, so I am guessing this is some sort of precision issue.
Is there a way to fix this? Technically I do not need to use eigh, but I think this is the underlying problem with my other errors (taking square roots of these matrices, getting NaNs, etc.)
I think the issue is that you've hit the limits of floating-point precision. A good rule-of-thumb for linear algebra results is that they're good to about one part in 10^8 for float32, and about one part in 10^16 for float 64. It appears that the ratio of your smallest to largest eigenvalue here is less than 10^-16. Because of this, the returned value cannot really be trusted and will depend on the details of the eigenvalue implementation you use.
For example, here are four different solvers you should have available; take a look at their results:
# 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')
Output:
[ -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
Notice they all agree on the larger two eigenvalues, but that the value of the smallest eigenvalue changes from implementation to implementation. Still, in all four cases the input matrix can be reconstructed up to 64-bit precision: this means the algorithms are operating as expected up to the precision available to them.
这篇关于scipy eigh为正半定矩阵给出负特征值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
