numpy/scipy eigendecompositions的不穩定結果

Question

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

>>> 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中的類似的代碼似乎沒有這個問題。

Answer 1

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

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

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

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

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

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

Answer 2

我認爲你的問題是你期待特徵值以特定的順序返回，並且它們並不總是相同的。對它們進行排序，然後你就會在路上。如果我運行代碼生成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)