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中的例子可以在这里查看（抱歉，我不知道怎么发布IPython的输出）

编辑 1：

这并不是像@Kh40tiK所说的那样，排序特征值的问题。可以在这里查看。

编辑 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 在你的系统中不起作用，可能是因为版本不同或者有bug。

Answer 2

通过 eigs 得到的特征值是随机顺序的。如果你把它们排序，你会发现每次得到的特征值都是一样的（当然，如果起始向量运气不好，可能会有例外）。

ARPACK 默认使用一个随机的起始向量来进行 Krylov 过程，这就是为什么每次调用时结果会不同的原因。如果你需要“可重复”的结果，可以指定 v0 参数。

需要注意的是，“可重复”这个词是加了引号的，因为由于编译器的优化和浮点数的舍入误差，结果并不总是能完全重复。