BUG: sparse.linalg: mistake in unsymm. real shift-invert ARPACK eigen…

…value selection

For real-valued unsymmetric problems, ARPACK can in some cases return
one extra eigenvalue.

The code that drops the extraneous eigenvalue did not take into account
that in shift-invert mode, the which= keyword applies to the
*transformed* eigenvalues, and dropped sometimes the wrong eigenvalue.

Moreover, fix sort order in this case for 'LM', 'LR', 'LI' modes
(largest items come first).

scipy/sparse/linalg/eigen/arpack/arpack.py

@@ -853,22 +853,27 @@ def extract(self, return_eigenvectors):
                 z = z[:, :nreturned]
             else:
                 # we got one extra eigenvalue (likely a cc pair, but which?)
-                # cut at approx precision for sorting
-                rd = np.round(d, decimals=_ndigits[self.tp])
+                if self.mode in (1, 2):
+                    rd = d
+                elif self.mode in (3, 4):
+                    rd = 1 / (d - self.sigma)
+
                 if self.which in ['LR', 'SR']:
                     ind = np.argsort(rd.real)
                 elif self.which in ['LI', 'SI']:
                     # for LI,SI ARPACK returns largest,smallest
-                    # abs(imaginary) why?
+                    # abs(imaginary) (complex pairs come together)
                     ind = np.argsort(abs(rd.imag))
                 else:
                     ind = np.argsort(abs(rd))
+
                 if self.which in ['LR', 'LM', 'LI']:
-                    d = d[ind[-k:]]
-                    z = z[:, ind[-k:]]
-                if self.which in ['SR', 'SM', 'SI']:
-                    d = d[ind[:k]]
-                    z = z[:, ind[:k]]
+                    ind = ind[-k:][::-1]
+                elif self.which in ['SR', 'SM', 'SI']:
+                    ind = ind[:k]
+
+                d = d[ind]
+                z = z[:, ind]
         else:
             # complex is so much simpler...
             d, z, ierr =\

scipy/sparse/linalg/eigen/arpack/tests/test_arpack.py

@@ -7,6 +7,7 @@
 """
 
 import threading
+import itertools
 
 import numpy as np
 
@@ -17,7 +18,7 @@
 
 from numpy import dot, conj, random
 from scipy.linalg import eig, eigh, hilbert, svd
-from scipy.sparse import csc_matrix, csr_matrix, isspmatrix, diags
+from scipy.sparse import csc_matrix, csr_matrix, isspmatrix, diags, rand
 from scipy.sparse.linalg import LinearOperator, aslinearoperator
 from scipy.sparse.linalg.eigen.arpack import eigs, eigsh, svds, \
      ArpackNoConvergence, arpack
@@ -984,3 +985,42 @@ def test_eigsh_for_k_greater():
         # Test 'A' for different types
         assert_raises(TypeError, eigsh, aslinearoperator(A), k=4)
         assert_raises(TypeError, eigsh, A_sparse, M=M_dense, k=4)
+
+
+def test_real_eigs_real_k_subset():
+    np.random.seed(1)
+
+    n = 10
+    A = rand(n, n, density=0.5)
+    A.data *= 2
+    A.data -= 1
+
+    v0 = np.ones(n)
+
+    whichs = ['LM', 'SM', 'LR', 'SR', 'LI', 'SI']
+    dtypes = [np.float32, np.float64]
+
+    for which, sigma, dtype in itertools.product(whichs, [None, 0, 5], dtypes):
+        prev_w = np.array([], dtype=dtype)
+        eps = np.finfo(dtype).eps
+        for k in range(1, 9):
+            w, z = eigs(A.astype(dtype), k=k, which=which, sigma=sigma,
+                        v0=v0.astype(dtype), tol=0)
+            assert_allclose(np.linalg.norm(A.dot(z) - z * w), 0, atol=np.sqrt(eps))
+
+            # Check that the set of eigenvalues for `k` is a subset of that for `k+1`
+            dist = abs(prev_w[:,None] - w).min(axis=1)
+            assert_allclose(dist, 0, atol=np.sqrt(eps))
+
+            prev_w = w
+
+            # Check sort order
+            if sigma is None:
+                d = w
+            else:
+                d = 1 / (w - sigma)
+
+            if which == 'LM':
+                # ARPACK is systematic for 'LM', but sort order
+                # appears not well defined for other modes
+                assert np.all(np.diff(abs(d)) <= 1e-6)