def METHOD_NAME() -> None:
    """
    Check that when conjugate gradient is used to evaluate the quadratic term,
    the obtained solution is:
    1. Smaller than the solution computed by Cholesky decomposition
    2. Within the error tolerance of the solution computed by Cholesky
    """
    rng: np.random.RandomState = np.random.RandomState(999)
    max_error: float = 1e-2
    noise: float = 1e-2
    train, z, _ = data(rng)
    x, y = train
    k = SquaredExponential()
    K = k(x) + noise * tf.eye(x.shape[0], dtype=default_float())
    def inv_quad_term(K: tf.Tensor, y: tf.Tensor) -> tf.Tensor:
        """
        For PSD K, compute -0.5 * y.T K^{-1} y via Cholesky decomposition
        """
        L = tf.linalg.cholesky(K)
        Linvy = tf.linalg.triangular_solve(L, y)
        return -0.5 * tf.reduce_sum(tf.square(Linvy))
    cholesky_quad_term = inv_quad_term(K, y)
    cglb = CGLB(
        train,
        kernel=k,
        inducing_variable=z,
        noise_variance=noise,
        cg_tolerance=max_error,
        max_cg_iters=100,
        restart_cg_iters=10,
    )
    common = cglb._common_calculation()
    cglb_quad_term = cglb.quad_term(common)
    assert cglb_quad_term <= cholesky_quad_term
    assert np.abs(cglb_quad_term - cholesky_quad_term) <= max_error