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

【问题标题】：scipy eigh gives negative eigenvalues for positive semidefinite matrixscipy eigh 给出半正定矩阵的负特征值
【发布时间】：2016-08-17 14:44:37
【问题描述】：

我在使用 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 的返回语句中将np.float64 替换为np.float32，我得到

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

相反，所以我猜这是某种精度问题。

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

【问题讨论】：

如果我使用 LA.eig 而不是 'LA.eigh'，我会得到不同的特征值：[ 5.15260753e+18+0.j 3.22785571e+01+0.j 4.62855397e+16+0.j]
恕我直言，您的 Hess_R 函数不会返回实际的 Hessian 矩阵。所以eigh 在您的情况下返回错误结果。
@B.M.你能进一步解释你的意思吗？返回的函数是什么？
您可能会看到不同的结果，具体取决于您的 numpy 版本所链接的 LAPACK 实现。使用 OpenBLAS v2.18 我得到 3.14000000e+02, 4.62855397e+16, 5.15260753e+18 和 scipy.linalg.eigh 和 1.06862038e+02, 4.62855397e+16, 5.15260753e+18 和 numpy.linalg.eigh 的特征值，例如 H （差异也可能部分是由于打印出浮点值时精度损失在x 作为字符串）。
仅供参考：您可以将 np.ones(d*d).reshape(d,d) 替换为 np.ones((d, d))。

标签： python numpy scipy precision eigenvalue

【解决方案1】：

我认为问题在于您已经达到了浮点精度的极限。线性代数结果的一个很好的经验法则是，对于 float32，它们大约是 10^8 的一部分，对于 float 64，大约是 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 位精度：这意味着算法可以按照预期的精度运行。

【讨论】：