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

eigs 返回的特征值是随机顺序的。如果对它们进行排序，您应该会发现每次都得到相同的结果（除非起始向量运气不佳）。

默认情况下，ARPACK 对 Krylov 进程使用随机起始向量，这解释了为什么每次调用都会得到不同的结果。如果您需要“可重现”的结果，请指定 v0 参数。

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

Answer 2

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 在您的系统中无法解决问题，则可能是版本差异或错误。