scipy eighは正の半定理行列に対して負の固有値を与える

Question

私は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 

私はhess_Rのreturn文でnp.float32でnp.float64を交換した場合、印刷固有値が

[ -6.74055241e-271 4.62855397e+016 5.15260753e+018] 

私が代わりに

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

を取得しているので、私はこれを推測していますが、ある種のです精度の問題。

これを修正する方法はありますか？技術的に私は一口を使用する必要はありませんが、これは私の他のエラーの根本的な問題だと思う（これらの行列の平方根を取ったり、NaNなどを取得する）

Answer 1

私は問題は、浮動小数点精度の限界。線形代数結果の良い経験則は、float32の場合は10^8、float64の場合は10^16の約1つに適しているということです。最小のものから最大のものまでの比率ここでの固有値は10^-16未満です。このため、戻り値は実際には信頼できず、使用する固有値実装の詳細に依存します。

たとえば、ここには4つの異なるソルバーがあります。その結果を見てみましょう：

# 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') 

出力：彼らはすべての方が大きい2個の固有値に同意するが、実現への実装から最小の固有値の値が変化すること

[ -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ビットの精度で再構成できます。これは、アルゴリズムが使用可能な精度まで期待どおりに動作していることを意味します。