From bf4a183f50279274bda5ec477231a91ce215a538 Mon Sep 17 00:00:00 2001 From: Matt Haberland Date: Wed, 23 Jun 2021 22:45:03 -0700 Subject: [PATCH] DOC: sparse.linalg: add PROPACK documentation to svds --- .../reference/sparse.linalg.svds-propack.rst | 8 ++ scipy/sparse/linalg/eigen/_svds.py | 24 ++-- scipy/sparse/linalg/eigen/_svds_doc.py | 132 +++++++++++++++++- 3 files changed, 149 insertions(+), 15 deletions(-) create mode 100644 doc/source/reference/sparse.linalg.svds-propack.rst diff --git a/doc/source/reference/sparse.linalg.svds-propack.rst b/doc/source/reference/sparse.linalg.svds-propack.rst new file mode 100644 index 000000000000..abd213911940 --- /dev/null +++ b/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 diff --git a/scipy/sparse/linalg/eigen/_svds.py b/scipy/sparse/linalg/eigen/_svds.py index 85d53744b09c..17988b9d53a2 100644 --- a/scipy/sparse/linalg/eigen/_svds.py +++ b/scipy/sparse/linalg/eigen/_svds.py @@ -66,7 +66,8 @@ 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' ` 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'} @@ -74,12 +75,14 @@ def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None, or smallest magnitude ('SM') singular values. v0 : ndarray, optional The starting vector for iteration; see method-specific - documentation (:ref:`'arpack' ` or - :ref:`'lobpcg' `) for details. + documentation (:ref:`'arpack' `, + :ref:`'lobpcg' `), or + :ref:`'propack' ` for details. maxiter : int, optional Maximum number of iterations; see method-specific - documentation (:ref:`'arpack' ` or - :ref:`'lobpcg' `) for details. + documentation (:ref:`'arpack' `, + :ref:`'lobpcg' `), or + :ref:`'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' ` and - :ref:`'lobpcg' ` are supported. + :ref:`'arpack' `, + :ref:`'lobpcg' `, and + :ref:`'propack' ` are supported. Default: `'arpack'`. options : dict, optional A dictionary of solver-specific options. No solver-specific options @@ -231,9 +235,9 @@ def matmat_XH_X(x): 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=maxiter, - compute_u=jobu, compute_v=jobv, irl_mode=True, kmax=ncv, - v0=v0) + 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, diff --git a/scipy/sparse/linalg/eigen/_svds_doc.py b/scipy/sparse/linalg/eigen/_svds_doc.py index 3b97e76d43f2..fa7a8c8266c4 100644 --- a/scipy/sparse/linalg/eigen/_svds_doc.py +++ b/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' ` is also supported. + :ref:`'lobpcg' ` and + :ref:`'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' ` is also supported. + :ref:`'arpack' ` and + :ref:`'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' ` and + :ref:`'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