def svd_solve(U, s, V, b, s_tol=1e-15):
    """
    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
    """
    # 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