Merge pull request #18 from mdhaber/propack

Add PROPACK to scipy.sparse.linalg documentation, tests

doc/source/reference/sparse.linalg.svds-propack.rst

@@ -0,0 +1,8 @@
+.. _sparse.linalg.svds-propack:
+
+svds(solver='propack')
+----------------------------------------
+
+.. scipy-optimize:function:: scipy.sparse.linalg.svds
+   :impl: scipy.sparse.linalg.eigen._svds_doc._svds_propack_doc
+   :method: propack

scipy/sparse/linalg/eigen/_svds.py

@@ -66,20 +66,23 @@ def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None,
     ncv : int, optional
         When ``solver='arpack'``, this is the number of Lanczos vectors
         generated. See :ref:`'arpack' <sparse.linalg.svds-arpack>` for details.
-        When ``solver='lobpcg'``, this parameter is ignored.
+        When ``solver='lobpcg'`` or ``solver='propack'``, this parameter is
+        ignored.
     tol : float, optional
         Tolerance for singular values. Zero (default) means machine precision.
     which : {'LM', 'SM'}
         Which `k` singular values to find: either the largest magnitude ('LM')
         or smallest magnitude ('SM') singular values.
     v0 : ndarray, optional
         The starting vector for iteration; see method-specific
-        documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>` or
-        :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`) for details.
+        documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
+        :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
+        :ref:`'propack' <sparse.linalg.svds-propack>` for details.
     maxiter : int, optional
         Maximum number of iterations; see method-specific
-        documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>` or
-        :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`) for details.
+        documentation (:ref:`'arpack' <sparse.linalg.svds-arpack>`,
+        :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
+        :ref:`'propack' <sparse.linalg.svds-propack>` for details.
     return_singular_vectors : bool or str, optional
         Singular values are always computed and returned; this parameter
         controls the computation and return of singular vectors.
@@ -93,8 +96,9 @@ def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None,
 
     solver : str, optional
             The solver used.
-            :ref:`'arpack' <sparse.linalg.svds-arpack>` and
-            :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>` are supported.
+            :ref:`'arpack' <sparse.linalg.svds-arpack>`,
+            :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`, and
+            :ref:`'propack' <sparse.linalg.svds-propack>` are supported.
             Default: `'arpack'`.
     options : dict, optional
         A dictionary of solver-specific options. No solver-specific options
@@ -166,8 +170,8 @@ def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None,
     else:
         raise ValueError("which must be either 'LM' or 'SM'.")
 
-    if (not (isinstance(A, LinearOperator) or isspmatrix(A)
-             or is_pydata_spmatrix(A))):
+    if not (isinstance(A, LinearOperator) or isspmatrix(A)
+            or is_pydata_spmatrix(A)):
         A = np.asarray(A)
 
     n, m = A.shape
@@ -221,12 +225,27 @@ def matmat_XH_X(x):
         eigvals, eigvec = lobpcg(XH_X, X, tol=tol ** 2, maxiter=maxiter,
                                  largest=largest)
 
+    elif solver == 'propack':
+        from scipy.sparse.linalg import svdp
+        jobu = ((isinstance(return_singular_vectors, str)
+                 and 'u' in return_singular_vectors)
+                or (isinstance(return_singular_vectors, bool)
+                    and return_singular_vectors))
+        jobv = ((isinstance(return_singular_vectors, str)
+                 and 'v' in return_singular_vectors)
+                or (isinstance(return_singular_vectors, bool)
+                    and return_singular_vectors))
+        return svdp(A, k=k, tol=tol**2, which=which, maxiter=None,
+                    compute_u=jobu, compute_v=jobv, irl_mode=True,
+                    kmax=maxiter, v0=v0)
+
     elif solver == 'arpack' or solver is None:
         eigvals, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter,
                                 ncv=ncv, which=which, v0=v0)
 
     else:
-        raise ValueError("solver must be either 'arpack', or 'lobpcg'.")
+        raise ValueError("solver must be either 'arpack', 'lobpcg', or "
+                         "'propack'.")
 
     # Gramian matrices have real non-negative eigenvalues.
     eigvals = np.maximum(eigvals.real, 0)

scipy/sparse/linalg/eigen/_svds_doc.py

@@ -49,7 +49,9 @@ def _svds_arpack_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
 
     solver : str, optional
             This is the solver-specific documentation for ``solver='arpack'``.
-            :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>` is also supported.
+            :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>` and
+            :ref:`'propack' <sparse.linalg.svds-propack>`
+            are also supported.
 
     Returns
     -------
@@ -66,7 +68,7 @@ def _svds_arpack_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
 
     Notes
     -----
-    This is a naive implementation using ARPACK or LOBPCG as an eigensolver
+    This is a naive implementation using ARPACK as an eigensolver
     on ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on which one is more
     efficient.
 
@@ -132,7 +134,7 @@ def _svds_lobpcg_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
         Number of singular values and singular vectors to compute.
         Must satisfy ``1 <= k < min(M, N)``.
     ncv : int, optional
-        Ignored; only used for ``solver='arpack'``.
+        Ignored.
     tol : float, optional
         Tolerance for singular values. Zero (default) means machine precision.
     which : {'LM', 'SM'}
@@ -159,7 +161,9 @@ def _svds_lobpcg_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
 
     solver : str, optional
             This is the solver-specific documentation for ``solver='lobpcg'``.
-            :ref:`'arpack' <sparse.linalg.svds-arpack>` is also supported.
+            :ref:`'arpack' <sparse.linalg.svds-arpack>` and
+            :ref:`'propack' <sparse.linalg.svds-propack>`
+            are also supported.
 
     Returns
     -------
@@ -176,7 +180,7 @@ def _svds_lobpcg_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
 
     Notes
     -----
-    This is a naive implementation using ARPACK or LOBPCG as an eigensolver
+    This is a naive implementation using LOBPCG as an eigensolver
     on ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on which one is more
     efficient.
 
@@ -220,3 +224,121 @@ def _svds_lobpcg_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
     True
     """
     pass
+
+
+def _svds_propack_doc(A, k=6, ncv=None, tol=0, which='LM', v0=None,
+                     maxiter=None, return_singular_vectors=True,
+                     solver='propack'):
+    """
+    Partial singular value decomposition of a sparse matrix using PROPACK.
+
+    Compute the largest or smallest `k` singular values and corresponding
+    singular vectors of a sparse matrix `A`. The order in which the singular
+    values are returned is not guaranteed.
+
+    In the descriptions below, let ``M, N = A.shape``.
+
+    Parameters
+    ----------
+    A : sparse matrix or LinearOperator
+        Matrix to decompose. If `A` is a ``LinearOperator``
+        object, it must define both ``matvec`` and ``rmatvec`` methods.
+    k : int, default: 6
+        Number of singular values and singular vectors to compute.
+        Must satisfy ``1 <= k <= min(M, N)``.
+    ncv : int, optional
+        Ignored.
+    tol : float, optional
+        The desired relative accuracy for computed singular values.
+        Zero (default) means machine precision.
+    which : {'LM', 'SM'}
+        Which `k` singular values to find: either the largest magnitude ('LM')
+        or smallest magnitude ('SM') singular values. Note that choosing
+        ``which='SM'`` will force the ``irl`` option to be set ``True``.
+    v0 : ndarray, optional
+        Starting vector for iterations: must be of length ``A.shape[0]``.
+        If not specified, PROPACK will generate a starting vector.
+    maxiter : int, optional
+        Maximum number of iterations / maximal dimension of the Krylov
+        subspace. Default is ``5 * k``.
+    return_singular_vectors : bool or str, default: True
+        Singular values are always computed and returned; this parameter
+        controls the computation and return of singular vectors.
+
+        - ``True``: return singular vectors.
+        - ``False``: do not return singular vectors.
+        - ``"u"``: only return the left singular values, without computing the
+          right singular vectors (if ``N > M``).
+        - ``"vh"``: only return the right singular values, without computing
+          the left singular vectors (if ``N <= M``).
+
+    solver : str, optional
+            This is the solver-specific documentation for ``solver='propack'``.
+            :ref:`'arpack' <sparse.linalg.svds-arpack>` and
+            :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`
+            are also supported.
+
+    Returns
+    -------
+    u : ndarray, shape=(M, k)
+        Unitary matrix having left singular vectors as columns.
+        If `return_singular_vectors` is ``"vh"``, this variable is not
+        computed, and ``None`` is returned instead.
+    s : ndarray, shape=(k,)
+        The singular values.
+    vh : ndarray, shape=(k, N)
+        Unitary matrix having right singular vectors as rows.
+        If `return_singular_vectors` is ``"u"``, this variable is not computed,
+        and ``None`` is returned instead.
+
+    Notes
+    -----
+    This is an interface to the Fortran library PROPACK [1]_.
+
+    References
+    ----------
+
+    .. [1] Larsen, Rasmus Munk. "PROPACK-Software for large and sparse SVD
+       calculations." Available online. URL
+       http://sun. stanford. edu/rmunk/PROPACK (2004): 2008-2009.
+
+    Examples
+    --------
+    Construct a matrix ``A`` from singular values and vectors.
+
+    >>> from scipy.stats import ortho_group
+    >>> from scipy.sparse import csc_matrix, diags
+    >>> from scipy.sparse.linalg import svds
+    >>> rng = np.random.default_rng()
+    >>> orthogonal = csc_matrix(ortho_group.rvs(10, random_state=rng))
+    >>> s = [0.0001, 0.001, 3, 4, 5]  # singular values
+    >>> u = orthogonal[:, :5]         # left singular vectors
+    >>> vT = orthogonal[:, 5:].T      # right singular vectors
+    >>> A = u @ diags(s) @ vT
+
+    With only three singular values/vectors, the SVD approximates the original
+    matrix.
+
+    >>> u2, s2, vT2 = svds(A, k=3, solver='propack')
+    >>> A2 = u2 @ np.diag(s2) @ vT2
+    >>> np.allclose(A2, A.todense(), atol=1e-3)
+    True
+
+    With all five singular values/vectors, we can reproduce the original
+    matrix.
+
+    >>> u3, s3, vT3 = svds(A, k=5, solver='propack')
+    >>> A3 = u3 @ np.diag(s3) @ vT3
+    >>> np.allclose(A3, A.todense())
+    True
+
+    The singular values match the expected singular values, and the singular
+    values are as expected up to a difference in sign. Consequently, the
+    returned arrays of singular vectors must also be orthogonal.
+
+    >>> (np.allclose(s3, s) and
+    ...  np.allclose(np.abs(u3), np.abs(u.todense())) and
+    ...  np.allclose(np.abs(vT3), np.abs(vT.todense())))
+    True
+    """
+    pass

scipy/sparse/linalg/eigen/tests/test_svds.py

@@ -57,7 +57,7 @@ def test_svd_simple_real():
                   [0, 0, 1, 0]], float)
     z = csc_matrix(x)
 
-    for solver in [None, 'arpack', 'lobpcg']:
+    for solver in [None, 'arpack', 'lobpcg', 'propack']:
         for m in [x.T, x, y, z, z.T]:
             for k in range(1, min(m.shape)):
                 u, s, vh = sorted_svd(m, k)
@@ -66,7 +66,8 @@ def test_svd_simple_real():
                 m_hat = svd_estimate(u, s, vh)
                 sm_hat = svd_estimate(su, ss, svh)
 
-                assert_array_almost_equal_nulp(m_hat, sm_hat, nulp=1000)
+                assert_array_almost_equal_nulp(
+                    m_hat, sm_hat, nulp=1000 if solver != 'propack' else 1116)
 
 
 def test_svd_simple_complex():
@@ -80,7 +81,7 @@ def test_svd_simple_complex():
                   [0, 0, 1, 0]], complex)
     z = csc_matrix(x)
 
-    for solver in [None, 'arpack', 'lobpcg']:
+    for solver in [None, 'arpack', 'lobpcg', 'propack']:
         for m in [x, x.T.conjugate(), x.T, y, y.conjugate(), z, z.T]:
             for k in range(1, min(m.shape) - 1):
                 u, s, vh = sorted_svd(m, k)
@@ -89,7 +90,8 @@ def test_svd_simple_complex():
                 m_hat = svd_estimate(u, s, vh)
                 sm_hat = svd_estimate(su, ss, svh)
 
-                assert_array_almost_equal_nulp(m_hat, sm_hat, nulp=1000)
+                assert_array_almost_equal_nulp(
+                    m_hat, sm_hat, nulp=1000 if solver != 'propack' else 1575)
 
 
 def test_svd_maxiter():
@@ -116,7 +118,7 @@ def test_svd_which():
     x = hilbert(6)
     for which in ['LM', 'SM']:
         _, s, _ = sorted_svd(x, 2, which=which)
-        for solver in [None, 'arpack', 'lobpcg']:
+        for solver in [None, 'arpack', 'lobpcg', 'propack']:
             ss = svds(x, 2, which=which, return_singular_vectors=False,
                       solver=solver)
             ss.sort()
@@ -127,7 +129,7 @@ def test_svd_v0():
     # check that the v0 parameter works as expected
     x = np.array([[1, 2, 3, 4], [5, 6, 7, 8]], float)
 
-    for solver in [None, 'arpack', 'lobpcg']:
+    for solver in [None, 'arpack', 'lobpcg', 'propack']:
         u, s, vh = svds(x, 1, solver=solver)
         u2, s2, vh2 = svds(x, 1, v0=u[:, 0], solver=solver)
 
@@ -164,7 +166,7 @@ def test_svd_LM_ones_matrix():
     for n, m in (6, 5), (5, 5), (5, 6):
         for t in float, complex:
             A = np.ones((n, m), dtype=t)
-            for solver in [None, 'arpack', 'lobpcg']:
+            for solver in [None, 'arpack', 'lobpcg']:  # propack fails
                 U, s, VH = svds(A, k, solver=solver)
 
                 # Check some generic properties of svd.
@@ -182,7 +184,7 @@ def test_svd_LM_zeros_matrix():
     for n, m in (3, 4), (4, 4), (4, 3):
         for t in float, complex:
             A = np.zeros((n, m), dtype=t)
-            for solver in [None, 'arpack', 'lobpcg']:
+            for solver in [None, 'arpack', 'lobpcg', 'propack']:
                 U, s, VH = svds(A, k, solver=solver)
 
                 # Check some generic properties of svd.
@@ -198,7 +200,7 @@ def test_svd_LM_zeros_matrix_gh_3452():
     # Note that for complex dype the size of this matrix is too small for k=1.
     n, m, k = 4, 2, 1
     A = np.zeros((n, m))
-    for solver in [None, 'arpack', 'lobpcg']:
+    for solver in [None, 'arpack', 'lobpcg', 'propack']:
         U, s, VH = svds(A, k, solver=solver)
 
         # Check some generic properties of svd.
@@ -238,9 +240,12 @@ def reorder(args):
         A = np.random.RandomState(52).randn(n, m)
         L = CheckingLinearOperator(A)
 
-        v0 = np.ones(min(A.shape))
+        for solver in [None, 'arpack', 'lobpcg', 'propack']:
+            if solver == 'propack':
+                v0 = np.ones(n)
+            else:
+                v0 = np.ones(min(A.shape))
 
-        for solver in [None, 'arpack', 'lobpcg']:
             U1, s1, VH1 = reorder(svds(A, k, v0=v0, solver=solver))
             U2, s2, VH2 = reorder(svds(L, k, v0=v0, solver=solver))
 
@@ -254,9 +259,13 @@ def reorder(args):
         A = np.random.RandomState(1909).randn(n, m)
         L = CheckingLinearOperator(A)
 
-        for solver in [None, 'arpack', 'lobpcg']:
-            U1, s1, VH1 = reorder(svds(A, k, which="SM", solver=solver))
-            U2, s2, VH2 = reorder(svds(L, k, which="SM", solver=solver))
+        for solver in [None, 'arpack', 'lobpcg', 'propack']:
+            # TODO: arpack crashes when v0=v0, which="SM"
+            kwargs = {'v0': v0} if solver not in {None, 'arpack'} else {}
+            U1, s1, VH1 = reorder(svds(A, k, which="SM", solver=solver,
+                                       **kwargs))
+            U2, s2, VH2 = reorder(svds(L, k, which="SM", solver=solver,
+                                       **kwargs))
 
             assert_allclose(np.abs(U1), np.abs(U2))
             assert_allclose(s1, s2)
@@ -271,11 +280,9 @@ def reorder(args):
                 A = (rng.randn(n, m) + 1j * rng.randn(n, m)).astype(dt)
                 L = CheckingLinearOperator(A)
 
-                for solver in [None, 'arpack', 'lobpcg']:
-                    U1, s1, VH1 = reorder(svds(A, k, which="LM",
-                                               solver=solver))
-                    U2, s2, VH2 = reorder(svds(L, k, which="LM",
-                                               solver=solver))
+                for solver in [None, 'arpack', 'lobpcg', 'propack']:
+                    U1, s1, VH1 = reorder(svds(A, k, which="LM", solver=solver))
+                    U2, s2, VH2 = reorder(svds(L, k, which="LM", solver=solver))
 
                     assert_allclose(np.abs(U1), np.abs(U2), rtol=eps)
                     assert_allclose(s1, s2, rtol=eps)