Gaussian process `log_marginal_likelihood()` `theta` parameter: pass as `log(theta)`?

Discussed in https://github.com/scikit-learn/scikit-learn/discussions/24765

I opened this as a discussion, but it is probably more suited for an issue, as it may lead to a documentation fix.

System Info

System:
    python: 3.10.7 (main, Sep  8 2022, 14:34:29) [GCC 12.2.0]
executable: /usr/bin/python3
   machine: Linux-5.19.0-1-amd64-x86_64-with-glibc2.35

Python dependencies:
      sklearn: 1.1.2
          pip: 22.2
   setuptools: 59.6.0
        numpy: 1.21.5
        scipy: 1.8.1
       Cython: 0.29.32
       pandas: 1.3.5
   matplotlib: 3.5.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /usr/lib/x86_64-linux-gnu/libgomp.so.1.0.0
        version: None
    num_threads: 16

^{Originally posted by elcorto October 26, 2022}

Preliminaries

I'm a bit confused by GaussianProcess{Regressor,Classifier}'s internal hyperparameter representation and the nature of the theta argument of the log_marginal_likelihood() method.

The doc strings of these methods say

theta : array-like of shape (n_kernel_params,) default=None
    Kernel hyperparameters for which the log-marginal likelihood is
    evaluated. If None, the precomputed log_marginal_likelihood
    of ``self.kernel_.theta`` is returned.

However, the GP's internal hyper optimizer code path seems to work with log(theta).

>>> from sklearn.gaussian_process import GaussianProcessRegressor
>>> from sklearn.gaussian_process.kernels import RBF, WhiteKernel

>>> gp=GaussianProcessRegressor(kernel=RBF()+WhiteKernel())
>>> gp.fit(rand(100,3), rand(100))

>>> gp.kernel_
RBF(length_scale=59.8) + WhiteKernel(noise_level=0.0812)

>>> log_theta=gp.kernel_.theta
>>> log_theta
array([ 4.09017207, -2.51094817])

>>> exp(log_theta)
array([59.75017175,  0.08119122])

The value of the kernel_.theta attribute is actually log(theta). And indeed, the theta getter returns log(theta), so they are not stored in log format but only returned that way.

The getter's docs say

Returns the (flattened, log-transformed) non-fixed hyperparameters.

Note that theta are typically the log-transformed values of the
kernel's hyperparameters as this representation of the search space
is more amenable for hyperparameter search, as hyperparameters like
length-scales naturally live on a log-scale.

Returns
-------
theta : ndarray of shape (n_dims,)
    The non-fixed, log-transformed hyperparameters of the kernel

OK.

Question

When calling log_marginal_likelihood(theta), the code sets kernel.theta = theta, which calls the theta setter which does essentially exp(theta). So from that I'd assume that we need to call log_marginal_likelihood(log(theta)). Indeed, continuing with above's code

>>> gp.log_marginal_likelihood_value_
-19.58130329439676

>>> gp.log_marginal_likelihood(log_theta)
-19.58130329439676

shows that we need to pass log(theta) to log_marginal_likelihood(). Also when trying to plot log_marginal_likelihood on a grid of hyperparameters, I only get what looks like correct results if I pass in log(theta). If that is true, should the documentation of the log_marginal_likelihood() methods be adapted accordingly?

Thanks.

Start by reading the Gaussian process log_marginal_likelihood() documentation and the linked sklearn/gaussian_process/_gpr.py call site, then inspect the theta getter and setter in sklearn/gaussian_process/kernels.py. Confirm the representation used by both GaussianProcessRegressor and GaussianProcessClassifier, and update the method documentation so the expected parameter form is unambiguous.