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
```

<div type='discussions-op-text'>

<sup>Originally posted by **elcorto** October 26, 2022</sup>
## 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)`.

```py
>>> 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][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`][call_setter], which calls [the `theta` setter][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

```py
>>> 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.


[setter]: https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/kernels.py#L290
[getter]: https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/kernels.py#L266
[call_setter]: https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L543
</div>

COMMENT:
I went through the code once more and I'm now certain that the optimizer code path calls the equivalent of `log_marginal_likelihood(log(theta))`, via `log_marginal_likelihood(kernel_.theta)` where the `kernel_.theta` getter returns `log(theta)` as stated above.

Also I tested a GP fit and calculating the LML using the same data, same kernel etc. against gpytorch, which also confirms that in in sklearn one must use `log_marginal_likelihood(log(theta))` when feeding in `theta`.

I'll send a PR to update the docs.

Labels: Bug, module:gaussian_process