scipy eigh給出正半定矩陣的負特徵值

我有一些問題，scipy的eigh函數返回正半定矩陣的負特徵值。以下是MWE。scipy eigh給出正半定矩陣的負特徵值

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等）

來源

2016-04-24 angryavian

如果我使用'LA.eig'而不是'LA.eigh'，我會得到不同的特徵值：'[5.15260753e + 18 + 0.j 3.22785571e + 01 + 0.j 4.62855397e + 16 + 0.j ]' – Peaceful

恕我直言，你的'Hess_R'函數不會返回一個實際的Hessian矩陣。所以'eigh'在你的情況下返回錯誤的結果。 –

@ B.M。你能進一步解釋你的意思嗎？什麼是函數返回？ – angryavian

我認爲問題在於你已經擊中了浮點精度的限制。線性代數結果的一個很好的經驗法則是它們對於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位精度：這意味着算法按預期運行，直至達到其可用的精度。

來源

2016-04-24 19:28:28 jakevdp

scipy eigh給出正半定矩陣的負特徵值

回答

相關問題