code: def svd_solve(U, s, V, b, s_tol=1e-15):
      
      # Test shapes for efficient computations
      n = U.shape[0]
      assert(b.shape[0] == n)
      m = b.shape[1] if np.ndim(b) > 1 else 1
  
      # Auto clamp SVD based on threshold
      sclamp = np.maximum(s, s_tol)
  
      # Inversion factors
      ss = 1. / np.sqrt(sclamp)
      U2 = U * ss[np.newaxis, :]
      V2 = ss[:, np.newaxis] * V
  
      if m < n:
          # Few queries
          X = U2.dot(V2.dot(b))  # O(n^2 (2m))
      else:
          X = U2.dot(V2).dot(b)  # O(n^2 (m + n))
  
      return X
docs: Solve the system :math:`A X = b` for :math:`X`.
  
      Here :math:`A` is a positive semi-definite matrix using the singular value
      decomposition. This truncates the SVD so only dimensions corresponding to
      non-negative and sufficiently large singular values are used.
  
      Parameters
      ----------
      U: ndarray
          The :code:`U` factor of :code:`U, s, V = svd(A)` positive
          semi-definite matrix.
      s: ndarray
          The :code:`s` factor of :code:`U, s, V = svd(A)` positive
          semi-definite matrix.
      V: ndarray
          The :code:`V` factor of :code:`U, s, V = svd(A)` positive
          semi-definite matrix.
      b: ndarray
          An array or matrix
      s_tol: float
          Cutoff for small singular values. Singular values smaller than
          :code:`s_tol` are clamped to :code:`s_tol`.
  
      Returns
      -------
      X: ndarray
          The result of :math:`X = A^-1 b`
      okind: ndarray
          The indices of :code:`s` that are kept in the factorisation