scipy eigh给出正半定矩阵的负特征值

Question

我有一些问题，scipy的eigh函数返回正半定矩阵的负特征值。以下是MWE。

hess_R函数返回一个正半定矩阵（它是一个秩一矩阵和一个对角矩阵的和，都带有非负项）。

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.float32取代np.float64我得到

[ -5.42905303e+10 4.62854925e+16 5.15260506e+18] 

代替，所以我猜测这是某种形式的精确问题。

有没有办法解决这个问题？从技术上讲，我不需要使用eigh，但我认为这是我的其他错误的基本问题（取这些矩阵的平方根，得到NaN等）

Answer 1

我认为问题在于你已经击中了浮点精度的限制。线性代数结果的一个很好的经验法则是它们对于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位精度：这意味着算法按预期运行，直至达到其可用的精度。