Solve the system :math:`A X = b` for :math:`X` where :math:`A` is a
    positive semi-definite matrix.

    This first tries cholesky, and if numerically unstable with solve using a
    truncated SVD (see solve_posdef_svd).

    The log-determinant of :math:`A` is also returned since it requires minimal
    overhead.

    Parameters
    ----------
    A: ndarray
        A positive semi-definite matrix.
    b: ndarray
        An array or matrix

    Returns
    -------
    X: ndarray
        The result of :math:`X = A^-1 b`
    logdet: float
        The log-determinant of :math:`A`

def solve_posdef(A, b):
    """
    Solve the system :math:`A X = b` for :math:`X` where :math:`A` is a
    positive semi-definite matrix.

    This first tries cholesky, and if numerically unstable with solve using a
    truncated SVD (see solve_posdef_svd).

    The log-determinant of :math:`A` is also returned since it requires minimal
    overhead.

    Parameters
    ----------
    A: ndarray
        A positive semi-definite matrix.
    b: ndarray
        An array or matrix

    Returns
    -------
    X: ndarray
        The result of :math:`X = A^-1 b`
    logdet: float
        The log-determinant of :math:`A`
    """
    # Try cholesky for speed
    try:
        lower = False
        L = cholesky(A, lower=lower)
        if any(L.diagonal() < CHOLTHRESH):
            raise LinAlgError("Unstable cholesky factor detected")
        X = cho_solve((L, lower), b)
        logdet = cho_log_det(L)

    # Failed cholesky, use svd to do the inverse
    except LinAlgError:

        U, s, V = svd(A)
        X = svd_solve(U, s, V, b)
        logdet = svd_log_det(s)

    return X, logdet