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参数。

请注意，“可重现的”在吓唬人的引号中，因为由于编译器的优化，浮点舍入错误并不总是可重现的。