Scipy中eigs函數的不一致特徵值稀疏

Question

我使用scipy.sparse.linalg模塊中的eigs函數，發現一些不一致的結果。運行兩次相同的代碼會得到不同的結果，即np.allclose的輸出爲False。任何人都可以解釋爲什麼？

from scipy.sparse.linalg import eigs 
from scipy.sparse import spdiags 
import numpy as np 


n1 = 100 
x, dx = linspace(0, 2, n1, retstep=True) 
e1 = ones(n1) 
A = 1./(dx**2)*spdiags([e1, -2*e1, e1], [-1,0,1], n1, n1) 

np.allclose(eigs(A, 90)[0], eigs(A, 90)[0]) 

在IPython中的例子可以看出here（抱歉不知道如何發佈IPython的輸出）

編輯1：

這是不排序的特徵值作爲的問題由@ Kh40tiK建議。見here。

編輯2：

嘗試不同版本SciPy的和運行發表@ Kh40tiK與其他調用腳本scipy.show_config（）後，似乎與MKL編譯SciPy的版本是一個有過錯。

隨着MKL：

2.7.6 |Anaconda 1.9.1 (64-bit)| (default, Jan 17 2014, 10:13:17) 
[GCC 4.1.2 20080704 (Red Hat 4.1.2-54)] 
('numpy:', '1.8.1') 
('scipy:', '0.13.3') 
umfpack_info: 
    NOT AVAILABLE 
lapack_opt_info: 
    libraries = ['mkl_lapack95_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core',   'iomp5', 'pthread'] 
    library_dirs = ['/home/jpsilva/anaconda/lib'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['/home/jpsilva/anaconda/include'] 
blas_opt_info: 
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread'] 
    library_dirs = ['/home/jpsilva/anaconda/lib'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['/home/jpsilva/anaconda/include'] 
openblas_info: 
    NOT AVAILABLE 
lapack_mkl_info: 
    libraries = ['mkl_lapack95_lp64', 'mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread'] 
    library_dirs = ['/home/jpsilva/anaconda/lib'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['/home/jpsilva/anaconda/include'] 
blas_mkl_info: 
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread'] 
    library_dirs = ['/home/jpsilva/anaconda/lib'] 
    define_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['/home/jpsilva/anaconda/include'] 
mkl_info: 
    libraries = ['mkl_intel_lp64', 'mkl_intel_thread', 'mkl_core', 'iomp5', 'pthread'] 
    library_dirs = ['/home/jpsilva/anaconda/lib'] 
    efine_macros = [('SCIPY_MKL_H', None)] 
    include_dirs = ['/home/jpsilva/anaconda/include'] 
False 
False 
False 
False 
False 
False 
False 
False 

沒有MKL：

2.7.5+ (default, Feb 27 2014, 19:37:08) 
[GCC 4.8.1] 
('numpy:', '1.8.1') 
('scipy:', '0.13.3') 
umfpack_info: 
    NOT AVAILABLE 
atlas_threads_info: 
    NOT AVAILABLE 
blas_opt_info: 
    libraries = ['f77blas', 'cblas', 'atlas'] 
    library_dirs = ['/usr/lib/atlas-base'] 
    define_macros = [('ATLAS_INFO', '"\\"3.10.1\\""')] 
    language = c 
    include_dirs = ['/usr/include/atlas'] 
atlas_blas_threads_info: 
    NOT AVAILABLE 
openblas_info: 
    NOT AVAILABLE 
lapack_opt_info: 
    libraries = ['lapack', 'f77blas', 'cblas', 'atlas'] 
    library_dirs = ['/usr/lib/atlas-base/atlas', '/usr/lib/atlas-base'] 
    define_macros = [('ATLAS_INFO', '"\\"3.10.1\\""')] 
    language = f77 
    include_dirs = ['/usr/include/atlas'] 
atlas_info: 
    libraries = ['lapack', 'f77blas', 'cblas', 'atlas'] 
    library_dirs = ['/usr/lib/atlas-base/atlas', '/usr/lib/atlas-base'] 
    define_macros = [('ATLAS_INFO', '"\\"3.10.1\\""')] 
    language = f77 
    include_dirs = ['/usr/include/atlas'] 
lapack_mkl_info: 
    NOT AVAILABLE 
blas_mkl_info: 
    NOT AVAILABLE 
atlas_blas_info: 
    libraries = ['f77blas', 'cblas', 'atlas'] 
    library_dirs = ['/usr/lib/atlas-base'] 
    define_macros = [('ATLAS_INFO', '"\\"3.10.1\\""')] 
    language = c 
    include_dirs = ['/usr/include/atlas'] 
mkl_info: 
    NOT AVAILABLE 
True 
False 
True 
False 
True 
False 
True 
False 

Answer 1

sp.sparse.linalg.eigs()不一定返回有序的特徵值，這意味着結果特徵值可能是隨機的順序。在調用np.allclose之前，您可能需要對特徵值進行排序。

另外，嘗試不同的容忍np.allclose，如：

np.allclose(eigs(A, 90)[0]), eigs(A,90)[0], 1e-3, 1e-5) 

希望它能幫助。

編輯

我稍微修改了劇本上的Python 3（不IPython中），sort做的事情。

#!/usr/bin/python3 
import sys 
from scipy.sparse.linalg import eigs 
from scipy.sparse import spdiags 
import numpy as np 
import scipy as sp 

n1 = 100 
x, dx = np.linspace(0, 2, n1, retstep=True) 
e1 = np.ones(n1) 
A = 1./(dx**2)*spdiags([e1, -2*e1, e1], [-1,0,1], n1, n1) 

print(sys.version) 
print('numpy:', np.version.version) 
print('scipy:', sp.version.version) 
for i in range(4): 
    print (np.allclose(np.sort(eigs(A, 90)[0]), np.sort(eigs(A, 90)[0]))) 
    print (np.allclose(eigs(A, 90)[0], eigs(A, 90)[0])) 

輸出：

3.4.0 (default, Mar 22 2014, 22:51:25) 
[GCC 4.8.2] 
numpy: 1.9.0.dev-b80ef75 
scipy: 0.15.0.dev-c2b7308 
True 
False 
True 
False 
True 
False 
True 
False 

如果sort不會做的伎倆在你的系統中，它可能是一個版本差異或錯誤。

Answer 2

eigs返回的特徵值是隨機的。如果您對它們進行排序，您應該會發現每次都會得到相同的結果（禁止使用啓動向量的運氣不佳的案例）。

默認情況下，ARPACK對Krylov進程使用隨機啓動向量，這就解釋了爲什麼每次調用都會得到不同的結果。如果您需要「可重複」結果，請指定v0參數。

請注意，「可重現的」在嚇唬人的引號中，因爲由於編譯器的優化，浮點舍入錯誤並不總是可重現的。