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 6399edf0c794..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 @@ -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) 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 diff --git a/scipy/sparse/linalg/eigen/tests/test_svds.py b/scipy/sparse/linalg/eigen/tests/test_svds.py index 322e76839362..901e029223b3 100644 --- a/scipy/sparse/linalg/eigen/tests/test_svds.py +++ b/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)