Inconsistent eigenvalues from eigs function in Scipy sparse

Question

I am using eigs function from scipy.sparse.linalg module and found some inconsistent results. Running twice the same code gives different results, ie, the output of np.allclose is False. Can anyone explain why is that?

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])

The example in IPython can be seen here (Sorry don't know how to post IPython output)

Edit 1:

It is not a matter of sorting the eigenvalues as suggested by @Kh40tiK. See here.

Edit 2:

After trying different versions of Scipy and running the script posted by @Kh40tiK with an additional call to scipy.show_config(), it seems that the SciPy version compiled with MKL is the one at fault.

With 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

Without 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

score 4 · Answer 1 · answered Apr 02 '14 at 14:05

The eigenvalues returned by eigs are in random order. If you sort them, you should find that you get the same ones each time (barring cases with bad luck with starting vector).

ARPACK by default uses a random starting vector for the Krylov process, which explains why you get different results for each call. If you need "reproducible" results, specify the v0 argument.

Note that "reproducible" is in scare quotes, because due to compiler optimizations, floating point rounding error is not always reproducible.

Kh40tiK · Answer 2 · 2014-04-02T23:54:22.210

3

sp.sparse.linalg.eigs() doesn't necessarily return ordered eigenvalues, which means the result eigenvalues could be in random order. You might want to sort the eigenvalues before calling np.allclose.

Also, try different tolerance for np.allclose, eg:

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

Hope it helps.

EDIT

I slightly modified the script on Python 3 ( not IPython ), sort do matters.

#!/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]))

Output:

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

If sort doesn't do the trick in your system, it could be a version difference or a bug.

edited Apr 02 '14 at 23:54

answered Apr 02 '14 at 09:25

Kh40tiK

2,276
19
29

Well, for `np.eig()` I would accept it but as `scipy.sparse.linalg.eigs` calculates only a few sorted eigenvalues (the first 6 and by largest magnitude, by default) I would expect it to return ordered ones. – poeticcapybara Apr 02 '14 at 09:37
And sort didn't do the trick... Check [here](http://nbviewer.ipython.org/gist/PoeticCapybara/9931042) – poeticcapybara Apr 02 '14 at 09:42
1-you're asking `eigs` for `90` eigenvalues and not `6`. 2-the order of the returned eigenvalues is not mentioned in the documentation of `eig` so you cannot assume anything about it. – gg349 Apr 02 '14 at 09:47
I ran some tests for Scipy 0.13.3 (current version available on pip and conda) and the Scipy version built with MKL breaks down, while the standard one doesn't. @flebool 1-That is irrelevant 2-Completely agree. – poeticcapybara Apr 02 '14 at 13:57
1

Why does your answer begin with a comment about `np.eig`? The question is about `scipy.sparse.linalg.eigs`. Your comment applies to `scipy.sparse.linalg.eigs`, and sorting the eigenvalues worked for me, so the reference to `np.eig` is confusing. – Warren Weckesser Apr 02 '14 at 14:26

Inconsistent eigenvalues from eigs function in Scipy sparse

2 Answers2