def svds(A, k=<NUM_LIT:6>, ncv=None, tol=<NUM_LIT:0>, which='<STR_LIT>', v0=None,<EOL>maxiter=None, return_singular_vectors=True):
if not (isinstance(A, LinearOperator) or isspmatrix(A)):<EOL><INDENT>A = np.asarray(A)<EOL><DEDENT>n, m = A.shape<EOL>if isinstance(A, LinearOperator):<EOL><INDENT>if n > m:<EOL><INDENT>X_dot = A.matvec<EOL>X_matmat = A.matmat<EOL>XH_dot = A.rmatvec<EOL><DEDENT>else:<EOL><INDENT>X_dot = A.rmatvec<EOL>XH_dot = A.matvec<EOL>dtype = getattr(A, '<STR_LIT>', None)<EOL>if dtype is None:<EOL><INDENT>dtype = A.dot(np.zeros([m,<NUM_LIT:1>])).dtype<EOL><DEDENT>def X_matmat(V):<EOL><INDENT>out = np.empty((V.shape[<NUM_LIT:1>], m), dtype=dtype)<EOL>for i, col in enumerate(V.T):<EOL><INDENT>out[i, :] = A.rmatvec(col.reshape(-<NUM_LIT:1>, <NUM_LIT:1>)).T<EOL><DEDENT>return out.T<EOL><DEDENT><DEDENT><DEDENT>else:<EOL><INDENT>if n > m:<EOL><INDENT>X_dot = X_matmat = A.dot<EOL>XH_dot = _herm(A).dot<EOL><DEDENT>else:<EOL><INDENT>XH_dot = A.dot<EOL>X_dot = X_matmat = _herm(A).dot<EOL><DEDENT><DEDENT>def matvec_XH_X(x):<EOL><INDENT>return XH_dot(X_dot(x))<EOL><DEDENT>XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype,<EOL>shape=(min(A.shape), min(A.shape)))<EOL>eigvals, eigvec = eigsh(XH_X, k=k, tol=tol ** <NUM_LIT:2>, maxiter=maxiter,<EOL>ncv=ncv, which=which, v0=v0)<EOL>if which == '<STR_LIT>':<EOL><INDENT>eigvals = np.maximum(eigvals.real, <NUM_LIT:0>)<EOL>t = eigvec.dtype.char.lower()<EOL>factor = {'<STR_LIT:f>': <NUM_LIT>, '<STR_LIT:d>': <NUM_LIT>}<EOL>cond = factor[t] * np.finfo(t).eps<EOL>cutoff = cond * np.max(eigvals)<EOL>above_cutoff = (eigvals > cutoff)<EOL>nlarge = above_cutoff.sum()<EOL>nsmall = k - nlarge<EOL>slarge = np.sqrt(eigvals[above_cutoff])<EOL>s = np.zeros_like(eigvals)<EOL>s[:nlarge] = slarge<EOL>if not return_singular_vectors:<EOL><INDENT>return s<EOL><DEDENT>if n > m:<EOL><INDENT>vlarge = eigvec[:, above_cutoff]<EOL>ularge = X_matmat(vlarge) / slarge<EOL>vhlarge = _herm(vlarge)<EOL><DEDENT>else:<EOL><INDENT>ularge = eigvec[:, above_cutoff]<EOL>vhlarge = _herm(X_matmat(ularge) / slarge)<EOL><DEDENT>u = _augmented_orthonormal_cols(ularge, nsmall)<EOL>vh = _augmented_orthonormal_rows(vhlarge, nsmall)<EOL><DEDENT>else:<EOL><INDENT>s = np.sqrt(eigvals)<EOL>if not return_singular_vectors:<EOL><INDENT>return s<EOL><DEDENT>if n > m:<EOL><INDENT>v = eigvec<EOL>u = X_matmat(v) / s<EOL>vh = _herm(v)<EOL><DEDENT>else:<EOL><INDENT>u = eigvec<EOL>vh = _herm(X_matmat(u) / s)<EOL><DEDENT><DEDENT>return u, s, vh<EOL>
Compute the largest k singular values/vectors for a sparse matrix.

    Parameters
    ----------
    A : {sparse matrix, LinearOperator}
        Array to compute the SVD on, of shape (M, N)
    k : int, optional
        Number of singular values and vectors to compute.
    ncv : int, optional
        The number of Lanczos vectors generated
        ncv must be greater than k+1 and smaller than n;
        it is recommended that ncv > 2*k
        Default: ``min(n, 2*k + 1)``
    tol : float, optional
        Tolerance for singular values. Zero (default) means machine precision.
    which : str, ['LM' | 'SM'], optional
        Which `k` singular values to find:

            - 'LM' : largest singular values
            - 'SM' : smallest singular values

        .. versionadded:: 0.12.0
    v0 : ndarray, optional
        Starting vector for iteration, of length min(A.shape). Should be an
        (approximate) right singular vector if N > M and a right singular vector
        otherwise.
        Default: random

        .. versionadded:: 0.12.0
    maxiter: int, optional
        Maximum number of iterations.

        .. versionadded:: 0.12.0
    return_singular_vectors : bool, optional
        Return singular vectors (True) in addition to singular values

        .. versionadded:: 0.12.0

    Returns
    -------
    u : ndarray, shape=(M, k)
        Unitary matrix having left singular vectors as columns.
    s : ndarray, shape=(k,)
        The singular values.
    vt : ndarray, shape=(k, N)
        Unitary matrix having right singular vectors as rows.

    Notes
    -----
    This is a naive implementation using ARPACK as an eigensolver
    on A.H * A or A * A.H, depending on which one is more efficient.