code: def svds(
      A,
      k=6,
      ncv=None,
      tol=0,
      which="LM",
      v0=None,
      maxiter=None,
      return_singular_vectors=True,
      solver="arpack",
      random_state=None,
      options=None,
  ):
      """
      Partial singular value decomposition of a sparse matrix.
  
      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 : ndarray, sparse matrix, or LinearOperator
          Matrix to decompose of a floating point numeric dtype.
      k : int, default: 6
          Number of singular values and singular vectors to compute.
          Must satisfy ``1 <= k <= kmax``, where ``kmax=min(M, N)`` for
          ``solver='propack'`` and ``kmax=min(M, N) - 1`` otherwise.
      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'`` 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>`,
          :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>`,
          :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`), or
          :ref:`'propack' <sparse.linalg.svds-propack>` for details.
      return_singular_vectors : {True, False, "u", "vh"}
          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"``: if ``M <= N``, compute only the left singular vectors and
            return ``None`` for the right singular vectors. Otherwise, compute
            all singular vectors.
          - ``"vh"``: if ``M > N``, compute only the right singular vectors and
            return ``None`` for the left singular vectors. Otherwise, compute
            all singular vectors.
  
          If ``solver='propack'``, the option is respected regardless of the
          matrix shape.
  
      solver :  {'arpack', 'propack', 'lobpcg'}, optional
              The solver used.
              :ref:`'arpack' <sparse.linalg.svds-arpack>`,
              :ref:`'lobpcg' <sparse.linalg.svds-lobpcg>`, and
              :ref:`'propack' <sparse.linalg.svds-propack>` are supported.
              Default: `'arpack'`.
      random_state : {None, int, `numpy.random.Generator`,
                      `numpy.random.RandomState`}, optional
  
          Pseudorandom number generator state used to generate resamples.
  
          If `random_state` is ``None`` (or `np.random`), the
          `numpy.random.RandomState` singleton is used.
          If `random_state` is an int, a new ``RandomState`` instance is used,
          seeded with `random_state`.
          If `random_state` is already a ``Generator`` or ``RandomState``
          instance then that instance is used.
      options : dict, optional
          A dictionary of solver-specific options. No solver-specific options
          are currently supported; this parameter is reserved for future use.
  
      Returns
      -------
      u : ndarray, shape=(M, k)
          Unitary matrix having left singular vectors as columns.
      s : ndarray, shape=(k,)
          The singular values.
      vh : ndarray, shape=(k, N)
          Unitary matrix having right singular vectors as rows.
  
      Notes
      -----
      This is a naive implementation using ARPACK or LOBPCG as an eigensolver
      on the matrix ``A.conj().T @ A`` or ``A @ A.conj().T``, depending on
      which one is smaller size, followed by the Rayleigh-Ritz method
      as postprocessing; see
      Using the normal matrix, in Rayleigh-Ritz method, (2022, Nov. 19),
      Wikipedia, https://w.wiki/4zms.
  
      Alternatively, the PROPACK solver can be called.
  
      Choices of the input matrix `A` numeric dtype may be limited.
      Only ``solver="lobpcg"`` supports all floating point dtypes
      real: 'np.float32', 'np.float64', 'np.longdouble' and
      complex: 'np.complex64', 'np.complex128', 'np.clongdouble'.
      The ``solver="arpack"`` supports only
      'np.float32', 'np.float64', and 'np.complex128'.
  
      Examples
      --------
      Construct a matrix `A` from singular values and vectors.
  
      >>> import numpy as np
      >>> from scipy import sparse, linalg, stats
      >>> from scipy.sparse.linalg import svds, aslinearoperator, LinearOperator
  
      Construct a dense matrix `A` from singular values and vectors.
  
      >>> rng = np.random.default_rng(258265244568965474821194062361901728911)
      >>> orthogonal = stats.ortho_group.rvs(10, random_state=rng)
      >>> s = [1e-3, 1, 2, 3, 4]  # non-zero singular values
      >>> u = orthogonal[:, :5]         # left singular vectors
      >>> vT = orthogonal[:, 5:].T      # right singular vectors
      >>> A = u @ np.diag(s) @ vT
  
      With only four singular values/vectors, the SVD approximates the original
      matrix.
  
      >>> u4, s4, vT4 = svds(A, k=4)
      >>> A4 = u4 @ np.diag(s4) @ vT4
      >>> np.allclose(A4, A, atol=1e-3)
      True
  
      With all five non-zero singular values/vectors, we can reproduce
      the original matrix more accurately.
  
      >>> u5, s5, vT5 = svds(A, k=5)
      >>> A5 = u5 @ np.diag(s5) @ vT5
      >>> np.allclose(A5, A)
      True
  
      The singular values match the expected singular values.
  
      >>> np.allclose(s5, s)
      True
  
      Since the singular values are not close to each other in this example,
      every singular vector matches as expected up to a difference in sign.
  
      >>> (np.allclose(np.abs(u5), np.abs(u)) and
      ...  np.allclose(np.abs(vT5), np.abs(vT)))
      True
  
      The singular vectors are also orthogonal.
  
      >>> (np.allclose(u5.T @ u5, np.eye(5)) and
      ...  np.allclose(vT5 @ vT5.T, np.eye(5)))
      True
  
      If there are (nearly) multiple singular values, the corresponding
      individual singular vectors may be unstable, but the whole invariant
      subspace containing all such singular vectors is computed accurately
      as can be measured by angles between subspaces via 'subspace_angles'.
  
      >>> rng = np.random.default_rng(178686584221410808734965903901790843963)
      >>> s = [1, 1 + 1e-6]  # non-zero singular values
      >>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
      >>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
      >>> vT = v.T
      >>> A = u @ np.diag(s) @ vT
      >>> A = A.astype(np.float32)
      >>> u2, s2, vT2 = svds(A, k=2, random_state=rng)
      >>> np.allclose(s2, s)
      True
  
      The angles between the individual exact and computed singular vectors
      may not be so small. To check use:
  
      >>> (linalg.subspace_angles(u2[:, :1], u[:, :1]) +
      ...  linalg.subspace_angles(u2[:, 1:], u[:, 1:]))
      array([0.06562513])  # may vary
      >>> (linalg.subspace_angles(vT2[:1, :].T, vT[:1, :].T) +
      ...  linalg.subspace_angles(vT2[1:, :].T, vT[1:, :].T))
      array([0.06562507])  # may vary
  
      As opposed to the angles between the 2-dimensional invariant subspaces
      that these vectors span, which are small for rights singular vectors
  
      >>> linalg.subspace_angles(u2, u).sum() < 1e-6
      True
  
      as well as for left singular vectors.
  
      >>> linalg.subspace_angles(vT2.T, vT.T).sum() < 1e-6
      True
  
      The next example follows that of 'sklearn.decomposition.TruncatedSVD'.
  
      >>> rng = np.random.RandomState(0)
      >>> X_dense = rng.random(size=(100, 100))
      >>> X_dense[:, 2 * np.arange(50)] = 0
      >>> X = sparse.csr_matrix(X_dense)
      >>> _, singular_values, _ = svds(X, k=5, random_state=rng)
      >>> print(singular_values)
      [ 4.3293...  4.4491...  4.5420...  4.5987... 35.2410...]
  
      The function can be called without the transpose of the input matrix
      ever explicitly constructed.
  
      >>> rng = np.random.default_rng(102524723947864966825913730119128190974)
      >>> G = sparse.rand(8, 9, density=0.5, random_state=rng)
      >>> Glo = aslinearoperator(G)
      >>> _, singular_values_svds, _ = svds(Glo, k=5, random_state=rng)
      >>> _, singular_values_svd, _ = linalg.svd(G.toarray())
      >>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
      True
  
      The most memory efficient scenario is where neither
      the original matrix, nor its transpose, is explicitly constructed.
      Our example computes the smallest singular values and vectors
      of 'LinearOperator' constructed from the numpy function 'np.diff' used
      column-wise to be consistent with 'LinearOperator' operating on columns.
  
      >>> diff0 = lambda a: np.diff(a, axis=0)
  
      Let us create the matrix from 'diff0' to be used for validation only.
  
      >>> n = 5  # The dimension of the space.
      >>> M_from_diff0 = diff0(np.eye(n))
      >>> print(M_from_diff0.astype(int))
      [[-1  1  0  0  0]
       [ 0 -1  1  0  0]
       [ 0  0 -1  1  0]
       [ 0  0  0 -1  1]]
  
      The matrix 'M_from_diff0' is bi-diagonal and could be alternatively
      created directly by
  
      >>> M = - np.eye(n - 1, n, dtype=int)
      >>> np.fill_diagonal(M[:,1:], 1)
      >>> np.allclose(M, M_from_diff0)
      True
  
      Its transpose
  
      >>> print(M.T)
      [[-1  0  0  0]
       [ 1 -1  0  0]
       [ 0  1 -1  0]
       [ 0  0  1 -1]
       [ 0  0  0  1]]
  
      can be viewed as the incidence matrix; see
      Incidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU,
      of a linear graph with 5 vertices and 4 edges. The 5x5 normal matrix
      ``M.T @ M`` thus is
  
      >>> print(M.T @ M)
      [[ 1 -1  0  0  0]
       [-1  2 -1  0  0]
       [ 0 -1  2 -1  0]
       [ 0  0 -1  2 -1]
       [ 0  0  0 -1  1]]
  
      the graph Laplacian, while the actually used in 'svds' smaller size
      4x4 normal matrix ``M @ M.T``
  
      >>> print(M @ M.T)
      [[ 2 -1  0  0]
       [-1  2 -1  0]
       [ 0 -1  2 -1]
       [ 0  0 -1  2]]
  
      is the so-called edge-based Laplacian; see
      Symmetric Laplacian via the incidence matrix, in Laplacian matrix,
      (2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.
  
      The 'LinearOperator' setup needs the options 'rmatvec' and 'rmatmat'
      of multiplication by the matrix transpose ``M.T``, but we want to be
      matrix-free to save memory, so knowing how ``M.T`` looks like, we
      manually construct the following function to be
      used in ``rmatmat=diff0t``.
  
      >>> def diff0t(a):
      ...     if a.ndim == 1:
      ...         a = a[:,np.newaxis]  # Turn 1D into 2D array
      ...     d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
      ...     d[0, :] = - a[0, :]
      ...     d[1:-1, :] = a[0:-1, :] - a[1:, :]
      ...     d[-1, :] = a[-1, :]
      ...     return d
  
      We check that our function 'diff0t' for the matrix transpose is valid.
  
      >>> np.allclose(M.T, diff0t(np.eye(n-1)))
      True
  
      Now we setup our matrix-free 'LinearOperator' called 'diff0_func_aslo'
      and for validation the matrix-based 'diff0_matrix_aslo'.
  
      >>> def diff0_func_aslo_def(n):
      ...     return LinearOperator(matvec=diff0,
      ...                           matmat=diff0,
      ...                           rmatvec=diff0t,
      ...                           rmatmat=diff0t,
      ...                           shape=(n - 1, n))
      >>> diff0_func_aslo = diff0_func_aslo_def(n)
      >>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)
  
      And validate both the matrix and its transpose in 'LinearOperator'.
  
      >>> np.allclose(diff0_func_aslo(np.eye(n)),
      ...             diff0_matrix_aslo(np.eye(n)))
      True
      >>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
      ...             diff0_matrix_aslo.T(np.eye(n-1)))
      True
  
      Having the 'LinearOperator' setup validated, we run the solver.
  
      >>> n = 100
      >>> diff0_func_aslo = diff0_func_aslo_def(n)
      >>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')
  
      The singular values squared and the singular vectors are known
      explicitly; see
      Pure Dirichlet boundary conditions, in
      Eigenvalues and eigenvectors of the second derivative,
      (2022, Nov. 19), Wikipedia, https://w.wiki/5YX6,
      since 'diff' corresponds to first
      derivative, and its smaller size n-1 x n-1 normal matrix
      ``M @ M.T`` represent the discrete second derivative with the Dirichlet
      boundary conditions. We use these analytic expressions for validation.
  
      >>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
      >>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
      ...                              np.arange(1, 4)) / n)
      >>> np.allclose(s, se, atol=1e-3)
      True
      >>> print(np.allclose(np.abs(u), np.abs(ue), atol=1e-6))
      True
  
      """
      args = _iv(
          A,
          k,
          ncv,
          tol,
          which,
          v0,
          maxiter,
          return_singular_vectors,
          solver,
          random_state,
      )
      (
          A,
          k,
          ncv,
          tol,
          which,
          v0,
          maxiter,
          return_singular_vectors,
          solver,
          random_state,
      ) = args
  
      largest = which == "LM"
      n, m = A.shape
  
      if n >= m:
          X_dot = A.matvec
          X_matmat = A.matmat
          XH_dot = A.rmatvec
          XH_mat = A.rmatmat
          transpose = False
      else:
          X_dot = A.rmatvec
          X_matmat = A.rmatmat
          XH_dot = A.matvec
          XH_mat = A.matmat
          transpose = True
  
          dtype = getattr(A, "dtype", None)
          if dtype is None:
              dtype = A.dot(np.zeros([m, 1])).dtype
  
      def matvec_XH_X(x):
          return XH_dot(X_dot(x))
  
      def matmat_XH_X(x):
          return XH_mat(X_matmat(x))
  
      XH_X = LinearOperator(
          matvec=matvec_XH_X,
          dtype=A.dtype,
          matmat=matmat_XH_X,
          shape=(min(A.shape), min(A.shape)),
      )
  
      # Get a low rank approximation of the implicitly defined gramian matrix.
      # This is not a stable way to approach the problem.
      if solver == "lobpcg":
          if k == 1 and v0 is not None:
              X = np.reshape(v0, (-1, 1))
          else:
              X = random_state.standard_normal(size=(min(A.shape), k))
  
          _, eigvec = lobpcg(
              XH_X, X, tol=tol**2, maxiter=maxiter, largest=largest
          )
  
      elif solver == "propack":
          if not HAS_PROPACK:
              raise ValueError(
                  "`solver='propack'` is opt-in due "
                  "to potential issues on Windows, "
                  "it can be enabled by setting the "
                  "`SCIPY_USE_PROPACK` environment "
                  "variable before importing scipy"
              )
          jobu = return_singular_vectors in {True, "u"}
          jobv = return_singular_vectors in {True, "vh"}
          irl_mode = which == "SM"
          res = _svdp(
              A,
              k=k,
              tol=tol**2,
              which=which,
              maxiter=None,
              compute_u=jobu,
              compute_v=jobv,
              irl_mode=irl_mode,
              kmax=maxiter,
              v0=v0,
              random_state=random_state,
          )
  
          u, s, vh, _ = res  # but we'll ignore bnd, the last output
  
          # PROPACK order appears to be largest first. `svds` output order is not
          # guaranteed, according to documentation, but for ARPACK and LOBPCG
          # they actually are ordered smallest to largest, so reverse for
          # consistency.
          s = s[::-1]
          u = u[:, ::-1]
          vh = vh[::-1]
  
          u = u if jobu else None
          vh = vh if jobv else None
  
          if return_singular_vectors:
              return u, s, vh
          else:
              return s
  
      elif solver == "arpack" or solver is None:
          if v0 is None:
              v0 = random_state.standard_normal(size=(min(A.shape),))
          _, eigvec = eigsh(
              XH_X,
              k=k,
              tol=tol**2,
              maxiter=maxiter,
              ncv=ncv,
              which=which,
              v0=v0,
          )
          # arpack do not guarantee exactly orthonormal eigenvectors
          # for clustered eigenvalues, especially in complex arithmetic
          eigvec, _ = np.linalg.qr(eigvec)
  
      # the eigenvectors eigvec must be orthonomal here; see gh-16712
      Av = X_matmat(eigvec)
      if not return_singular_vectors:
          s = svd(Av, compute_uv=False, overwrite_a=True)
          return s[::-1]
  
      # compute the left singular vectors of X and update the right ones
      # accordingly
      u, s, vh = svd(Av, full_matrices=False, overwrite_a=True)
      u = u[:, ::-1]
      s = s[::-1]
      vh = vh[::-1]
  
      jobu = return_singular_vectors in {True, "u"}
      jobv = return_singular_vectors in {True, "vh"}
  
      if transpose:
          u_tmp = eigvec @ _herm(vh) if jobu else None
          vh = _herm(u) if jobv else None
          u = u_tmp
      else:
          if not jobu:
              u = None
          vh = vh @ _herm(eigvec) if jobv else None
  
      return u, s, vh
  
file_path: scipy/scipy/sparse/linalg/_eigen/_svds.py
imports: import os
  import numpy
  from .arpack import _arpack
  from .None import eigsh
  from scipy._lib._util import check_random_state
  from scipy.sparse.linalg._interface import LinearOperator, aslinearoperator
  from scipy.sparse.linalg._eigen.lobpcg import lobpcg
  from scipy.linalg import svd
label_desc: Use image node_id 3 for calling a global function with example usage: svds(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors, solver, random_state, options) and returns: u,  s,  vh, XH_dot, XH_mat, s, u,  s,  vh, s