Repository: scikit-learn/scikit-learn
File: sklearn/gaussian_process/tests/test_gpc.py (Python)

Code before:
def test_lml_gradient(kernel):
    # Compare analytic and numeric gradient of log marginal likelihood.
    gpc = GaussianProcessClassifier(kernel=kernel).fit(X, y)

    length_scales = np.linspace(1, 25, 1_000)

    def evaluate_grad_at_length_scales(length_scales):
        result = np.zeros_like(length_scales)
        for i, length_scale in enumerate(length_scales):
            kernel.length_scale = length_scale
            result[i] = (
                gpc.log_marginal_likelihood(kernel.theta)
                if len(kernel.theta) == 1
                else [gpc.log_marginal_likelihood([theta]) for theta in kernel.theta]
            )
        return result

    lml_gradient = np.zeros_like(length_scales)

    for i, length_scale in enumerate(length_scales):
        kernel.length_scale = length_scale
        lml_gradient[i] = gpc.log_marginal_likelihood(kernel.theta, eval_gradient=True)[
            1
        ][0]

    lml_gradient_manual = derivative(
        evaluate_grad_at_length_scales, length_scales, maxiter=20
    ).df

    assert_almost_equal(lml_gradient, lml_gradient_manual, 3)


def test_random_starts(global_random_seed):
    # Test that an increasing number of random-starts of GP fitting only
    # increases the log marginal likelihood of the chosen theta.
    n_samples, n_features = 25, 2
    rng = np.random.RandomState(global_random_seed)
    X = rng.randn(n_samples, n_features) * 2 - 1
    y = (np.sin(X).sum(axis=1) + np.sin(3 * X).sum(axis=1)) > 0

    kernel = C(1.0, (1e-2, 1e2)) * RBF(
        length_scale=[1e-3] * n_features, length_scale_bounds=[(1e-4, 1e2)] * n_features
    )
    last_lml = -np.inf
    for n_restarts_optimizer in range(5):
        gp = GaussianProcessClassifier(
            kernel=kernel,
            n_restarts_optimizer=n_restarts_optimizer,
            random_state=global_random_seed,
        ).fit(X, y)
        lml = gp.log_marginal_likelihood(gp.kernel_.theta)

Review comment (suggestion):
I prefer to keep using the `lml_gradient_approx` name for this variable, as `scipy.differentiate.derivative` is doing an approximate (numerical) computation rather than a manually implemented (symbolically derived) computation.

Code after:
        | (gpc.kernel_.theta == gpc.kernel_.bounds[:, 1])
    )


@pytest.mark.xfail(
    raises=AssertionError,
    reason="https://github.com/scikit-learn/scikit-learn/issues/31366",
)
@pytest.mark.parametrize("kernel", non_fixed_kernels)
@pytest.mark.parametrize("length_scale", np.logspace(-3, 3, 13))
def test_lml_gradient(kernel, length_scale):
    # Clone the kernel object prior to mutating it to avoid any side effects between
    # GP tests:
    kernel = clone(kernel)
    length_scale_param_name = next(
        name for name in kernel.get_params() if name.endswith("length_scale")
    )
    kernel.set_params(**{length_scale_param_name: length_scale})

    # Compare analytic and numeric gradient of log marginal likelihood.
    gpr = GaussianProcessClassifier(kernel=kernel).fit(X, y)
    _, lml_gradient = gpr.log_marginal_likelihood(kernel.theta, eval_gradient=True)
    epsilon = 1e-9
    lml_gradient_approx = approx_fprime(
        kernel.theta.copy(),
        lambda theta: gpr.log_marginal_likelihood(theta, False),
        epsilon=epsilon,
    )
    assert_allclose(lml_gradient, lml_gradient_approx, rtol=1e-4, atol=epsilon * 100)


def test_random_starts(global_random_seed):
    # Test that an increasing number of random-starts of GP fitting only
    # increases the log marginal likelihood of the chosen theta.
    n_samples, n_features = 25, 2
    rng = np.random.RandomState(global_random_seed)
    X = rng.randn(n_samples, n_features) * 2 - 1
    y = (np.sin(X).sum(axis=1) + np.sin(3 * X).sum(axis=1)) > 0

    kernel = C(1.0, (1e-2, 1e2)) * RBF(
        length_scale=[1e-3] * n_features, length_scale_bounds=[(1e-4, 1e2)] * n_features
    )
    last_lml = -np.inf
    for n_restarts_optimizer in range(5):
        gp = GaussianProcessClassifier(
            kernel=kernel,
            n_restarts_optimizer=n_restarts_optimizer,
            random_state=global_random_seed,
        ).fit(X, y)
        lml = gp.log_marginal_likelihood(gp.kernel_.theta)
        assert lml > last_lml - np.finfo(np.float32).eps