numpy/scipy eigendecompositions的不穩定結果

我發現scipy.linalg.eig有時會給出不一致的結果。但不是每次。numpy/scipy eigendecompositions的不穩定結果

>>> import numpy as np 
>>> import scipy.linalg as lin 
>>> modmat=np.random.random((150,150)) 
>>> modmat=modmat+modmat.T # the data i am interested in is described by real symmetric matrices 
>>> d,v=lin.eig(modmat) 
>>> dx=d.copy() 
>>> vx=v.copy() 
>>> d,v=lin.eig(modmat) 
>>> np.all(d==dx) 
False 
>>> np.all(v==vx) 
False 
>>> e,w=lin.eigh(modmat) 
>>> ex=e.copy() 
>>> wx=w.copy() 
>>> e,w=lin.eigh(modmat) 
>>> np.all(e==ex) 
True 
>>> e,w=lin.eigh(modmat) 
>>> np.all(e==ex) 
False

雖然我不是最偉大的線性代數嚮導，我不明白，是特徵分解本質上受到怪異的舍入誤差，但我不明白爲什麼重複計算會導致不同的值。但我的結果和可重複性是不同的。

問題的本質究竟是什麼 - 好的，有時結果是可以接受的不同，有時它們不是。以下是一些示例：

>>> d[1] 
(9.8986888573772465+0j) 
>>> dx[1] 
(9.8986888573772092+0j)

〜3e-13的上述差異似乎並不是什麼大問題。相反，真正的問題（至少對於我現在的項目來說）是某些特徵值似乎無法在正確的符號上達成一致。

>>> np.all(np.sign(d)==np.sign(dx)) 
False 
>>> np.nonzero(np.sign(d)!=np.sign(dx)) 
(array([ 38, 39, 40, 41, 42, 45, 46, 47, 79, 80, 81, 82, 83, 
    84, 109, 112]),) 
>>> d[38] 
(-6.4011617320002525+0j) 
>>> dx[38] 
(6.1888785138080209+0j)

在MATLAB中的類似的代碼似乎沒有這個問題。

來源

2013-02-04 aestrivex

（+1）有趣的... – NPE

我試圖重現此使用NumPy的1.6.1/SciPy的0.10.1，但不能。 – NPE

我使用numpy 1.6.1和scipy 0.10.0。另外，當我沒有使用copy（）時，我無法產生這個錯誤（但它確實存在於我的大型應用程序中，其中類似於copy（）的事情正在進行）。這並不意味着太多，因爲它的意思是不一致的。 – aestrivex

特徵值分解滿足AV = V LAMBDA，這是所有什麼是保證---例如特徵值的順序不是。

回答您的問題的第二部分：

現代編譯器/線性代數庫產生/包含的代碼做不同的事情根據數據是否被在存儲器排列在（例如）16字節邊界。這會影響計算中的舍入誤差，因爲浮點運算按不同的順序完成。如果算法（這裏，LAPACK/xGEEV）在這方面不能保證數值穩定性，那麼舍入誤差的小變化就可以影響諸如特徵值排序等事情。

（如果你的代碼是這樣的事情敏感，這是不正確在不同的平臺或不同的庫版本運行例如，它會導致類似的問題！）

結果通常是準確定性 - - 例如，你得到2個可能的結果之一，這取決於數組是否恰好在內存中對齊。如果您對路線感興趣，請檢查A.__array_interface__['data'][0] % 16。

見http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf更多

來源

2013-02-04 23:26:44

除了解釋之外，在這裏包含一個解決方案是非常有用的。我不知道如何/在哪裏使用'A .__ array_interface __ ['data'] [0]％16'，基本上，鑑於這是我的問題，我該如何更改我的代碼，以便它不敏感辦法？ –

我認爲你的問題是你期待特徵值以特定的順序返回，並且它們並不總是相同的。對它們進行排序，然後你就會在路上。如果我運行代碼生成d和dx與eig我得到如下：

>>> np.max(d - dx) 
(19.275224236664116+0j)

但是......

>>> d_i = np.argsort(d) 
>>> dx_i = np.argsort(dx) 
>>> np.max(d[d_i] - dx[dx_i]) 
(1.1368683772161603e-13+0j)

來源

2013-02-04 21:00:39 Jaime

啊。這不是正確的答案，但它足夠接近，它指出我在正確的方向。我的程序中發生了什麼事情，在我的代碼中的其他地方，我正在將特徵向量索引到它們不是的位置。它確實幫助我修復了我的錯誤（非常感謝），但它並沒有解釋我提出的問題，這就是爲什麼這個計算在多次運行時提供不同的答案或不同的順序。也就是說，沒問題，特徵值不能保證以任何特定的順序，但是如果你給它與輸入相同的矩陣，他們爲什麼會改變？ – aestrivex

numpy/scipy eigendecompositions的不穩定結果

回答

相關問題