Reference to pylops in LinearOperator notes and ARPACK example (scipy…

…#9605)

DOC: add reference to pylops in LinearOperator notes. 

Also added example with LinearOperator to Sparse Eigenvalue Problems with ARPACK example

doc/source/tutorial/arpack.rst

@@ -3,6 +3,8 @@ Sparse Eigenvalue Problems with ARPACK
 
 .. sectionauthor:: Jake Vanderplas <vanderplas@astro.washington.edu>
 
+.. sectionauthor:: Matteo Ravasi  <matteoravasi@gmail.com>
+
 .. currentmodule:: scipy.sparse.linalg
 
 
@@ -89,8 +91,8 @@ derived from :func:`scipy.sparse.linalg.LinearOperator`.  For this example, for
 simplicity, we'll construct a symmetric, positive-definite matrix.
 
     >>> import numpy as np
-    >>> from scipy.linalg import eigh
-    >>> from scipy.sparse.linalg import eigsh
+    >>> from scipy.linalg import eig, eigh
+    >>> from scipy.sparse.linalg import eigs, eigsh
     >>> np.set_printoptions(suppress=True)
     >>>
     >>> np.random.seed(0)
@@ -204,6 +206,105 @@ but the operation can also be specified by the user.  See the docstring of
 :func:`scipy.sparse.linalg.eigs` for details.
 
 
+Use of LinearOperator
+---------------------
+We consider now the case where you'd like to avoid creating a dense matrix
+and use :func:`scipy.sparse.linalg.LinearOperator` instead.  Our first
+linear operator applies element-wise multiplication between the input vector
+and a vector :math:`\mathbf{d}` provided by the user to the operator itself.
+This operator mimics a diagonal matrix with the elements of :math:`\mathbf{d}`
+along the main diagonal and it has the main benefit that the forward and
+adjoint operations are simple element-wise multiplications other
+than matrix-vector multiplications.  For a diagonal matrix, we expect the
+eigenvalues to be equal to the elements along the main diagonal, in this case
+:math:`\mathbf{d}`. The eigenvalues and eigenvectors obtained with ``eigsh``
+are compared  those obtained by using ``eigh`` when applied to
+the dense matrix:
+
+   >>> from scipy.sparse.linalg import LinearOperator
+   >>> class Diagonal(LinearOperator):
+   ...     def __init__(self, diag, dtype='float32'):
+   ...         self.diag = diag
+   ...         self.shape = (len(self.diag), len(self.diag))
+   ...         self.dtype = np.dtype(dtype)
+   ...     def _matvec(self, x):
+   ...         return self.diag*x
+   ...     def _rmatvec(self, x):
+   ...         return self.diag*x
+
+   >>> np.random.seed(0)
+   >>> N = 100
+   >>> d = np.random.normal(0, 1, N).astype(np.float64)
+   >>> D = np.diag(d)
+   >>> Dop = Diagonal(d, dtype=np.float64)
+
+   >>> evals_all, evecs_all = eigh(D)
+   >>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
+   >>> evals_all[-3:]
+   array([1.9507754 , 2.2408932 , 2.26975462])
+   >>> evals_large
+   array([1.9507754 , 2.2408932 , 2.26975462])
+   >>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
+   array([[-1.  0.  0.],     # may vary (signs)
+          [-0. -1.  0.],
+          [ 0.  0. -1.]]
+
+In this case we have created a quick and easy ``Diagonal`` operator.
+The external library `PyLops <https://pylops.readthedocs.io>`_ provides
+similar capabilities in the `Diagonal <https://pylops.readthedocs.io/en/
+latest/api/generated/pylops.Diagonal.html#pylops.Diagonal>`_ operator
+as well as several other operators.
+
+Finally, we consider a linear operator that mimics the application of a
+first derivative stencil. In this case the operator is equivalent to a real
+nonsymmetric matrix. Once again we compare the estimated eigenvalues
+and eigenvectors with those from a dense matrix that applies the
+same first derivative to an input signal:
+
+    >>> class FirstDerivative(LinearOperator):
+    ...     def __init__(self, N, dtype='float32'):
+    ...         self.N = N
+    ...         self.shape = (self.N, self.N)
+    ...         self.dtype = np.dtype(dtype)
+    ...     def _matvec(self, x):
+    ...         y = np.zeros(self.N, self.dtype)
+    ...         y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
+    ...         return y
+    ...     def _rmatvec(self, x):
+    ...         y = np.zeros(self.N, self.dtype)
+    ...         y[0:-2] = y[0:-2] - (0.5*x[1:-1])
+    ...         y[2:] = y[2:] + (0.5*x[1:-1])
+    ...         return y
+
+    >>> N = 21
+    >>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
+    >>> D[0] = D[-1] = 0 # take away edge effects
+    >>> Dop = FirstDerivative(N, dtype=np.float64)
+
+    >>> evals_all, evecs_all = eig(D)
+    >>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
+    >>> evals_all_imag = evals_all.imag
+    >>> isort_imag = np.argsort(np.abs(evals_all_imag))
+    >>> evals_all_imag = evals_all_imag[isort_imag]
+    >>> evals_large_imag = evals_large.imag
+    >>> isort_imag = np.argsort(np.abs(evals_large_imag))
+    >>> evals_large_imag = evals_large_imag[isort_imag]
+    >>> evals_all_imag[-4:]
+    array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
+    >>> evals_large_imag
+    array([0.95105652, -0.95105652, 0.98768834, -0.98768834])
+
+Note that the eigenvalues of this operator are all imaginary. Moreover,
+the keyword ``which='LI'`` of :func:`scipy.sparse.linalg.eigs` produces
+the eigenvalues with largest absolute imaginary part (both
+positive and negative). Again, a more advanced implementation of the
+first derivative operator is available in the
+`PyLops <https://pylops.readthedocs.io>`_ library under the name of
+`FirstDerivative <https://pylops.
+readthedocs.io/en/latest/api/generated/pylops.FirstDerivative.html>`_
+operator.
+
+
 References
 ----------
 .. [1] http://www.caam.rice.edu/software/ARPACK/

scipy/sparse/linalg/interface.py

@@ -113,6 +113,11 @@ class LinearOperator(object):
     is always a new, composite LinearOperator, that defers linear
     operations to the original operators and combines the results.
 
+    More details regarding how to subclass a LinearOperator and several
+    examples of concrete LinearOperator instances can be found in the
+    external project `PyLops <https://pylops.readthedocs.io>`_.
+
+
     Examples
     --------
     >>> import numpy as np