Gaussian Process regression get_params set_params does not yield same predictor?

### Describe the bug

This is on version 1.1.1

I *think* the example below is the correct way of setting GPR params saved from a previous .fit() tuning via the optimizer. Reading through the code I don't yet understand why the parameters `self.L_` and `self.alpha_` are not being updated but maybe am missing something in my reading.

See below for the example. Everything passes except the last line. 

### Steps/Code to Reproduce

```python
import sys

import numpy as np
from scipy.linalg import cholesky
from sklearn.gaussian_process import GaussianProcessRegressor, kernels
from sklearn.gaussian_process.kernels import *

np.random.seed(0)

X = np.random.randn(100, 3)
y = np.random.randn(100, 1)

KERNEL = Matern()
# make it once, fit and get the params
gp_ = GaussianProcessRegressor(kernel=KERNEL, normalize_y=True, alpha=1e-5)
gp_.fit(X, y)
params = gp_.kernel_.get_params(deep=False)  # tried deep=True but doesn't matter

# make it again, don't fit, but set the params
gp = GaussianProcessRegressor(kernel=KERNEL, normalize_y=True, alpha=1e-5)
gp.optimizer = None  # this should make fit do very little
gp.fit(X, y)  # this creates kernele

# do it again for good measure? Maybe we need to do this several times to trigger it?
gp.kernel_.set_params(**params)
gp.fit(X, y)  # do this again to now recreate L_
gp.kernel_.set_params(**params)


def get_K_L_(self):
    K = self.kernel_(self.X_train_)
    K[np.diag_indices_from(K)] += self.alpha
    L_ = cholesky(K, lower=True)
    return K, L_


Ka, La = get_K_L_(gp_)
Kb, Lb = get_K_L_(gp)

a = gp_.__getstate__()
b = gp.__getstate__()

# here you will see alpha_ and L_ are not the same. Something is not updating them correctly?
# see https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L322
for k in a:
    if isinstance(a[k], np.ndarray):
        print(k, type(a[k]), np.all(a[k] == b[k]))


assert np.all(gp_.X_train_ == gp.X_train_)  # passes
assert gp_.alpha == gp.alpha  # passes
assert np.allclose(Ka, Kb), 'Ka Kb do not match'  # passes
assert np.allclose(La, Lb), 'La Lb do not match'  # passes
# it is bizarre that the output of the kernels matchs but not L_????
assert np.allclose(gp.L_, gp_.L_), 'L_ do not match'  # FAILS
```

### Expected Results

No exceptions. The gp and gp_ should return the same results, have same state (except for the optimizer).

### Actual Results

```
_y_train_mean <class 'numpy.ndarray'> True
_y_train_std <class 'numpy.ndarray'> True
X_train_ <class 'numpy.ndarray'> True
y_train_ <class 'numpy.ndarray'> True
L_ <class 'numpy.ndarray'> False
alpha_ <class 'numpy.ndarray'> False
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-5-71be2b0a00fc> in <module>
     53 assert np.allclose(La, Lb), 'La Lb do not match'  # passes
     54 # it is bizarre that the output of the kernels matchs but not L_????
---> 55 assert np.allclose(gp.L_, gp_.L_), 'L_ do not match'  # FAILS

AssertionError: L_ do not match
```

### Versions

```shell
System:
    python: 3.10.4 (main, Mar 31 2022, 08:41:55) [GCC 7.5.0]
executable: /home/cottrell/anaconda3/envs/3.10.4/bin/python
   machine: Linux-5.15.0-43-generic-x86_64-with-glibc2.35

Python dependencies:
      sklearn: 1.1.1
          pip: 22.1.2
   setuptools: 61.2.0
        numpy: 1.23.1
        scipy: 1.8.1
       Cython: 0.29.30
       pandas: 1.4.3
   matplotlib: 3.5.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/cottrell/anaconda3/envs/3.10.4/lib/python3.10/site-packages/numpy.libs/libopenblas64_p-r0-742d56dc.3.20.so
        version: 0.3.20
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /home/cottrell/anaconda3/envs/3.10.4/lib/python3.10/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None
    num_threads: 12

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/cottrell/anaconda3/envs/3.10.4/lib/python3.10/site-packages/scipy.libs/libopenblasp-r0-8b9e111f.3.17.so
        version: 0.3.17
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/cottrell/anaconda3/envs/3.10.4/lib/python3.10/site-packages/opencv_python_headless.libs/libopenblas-r0-f650aae0.3.3.so
        version: None
threading_layer: disabled
   architecture: Haswell
    num_threads: 1
```

COMMENT:
> I think the misunderstanding is in the lines below. Calling `gp.fit` will start from scratch and recreate the kernel from its constructor parameters. It will override what you were trying to do with `gp.kernel_.set_params`:
> 
> ```python
> gp.kernel_.set_params(**params)
> gp.fit(X, y)  # do this again to now recreate L_
> gp.kernel_.set_params(**params)
> ```
> 
> You can add print statements before and after the `fit` to convince yourself.

No, that fit with optimizer set to None does not redo things. It might re-calculate the kernel but it should not set new hyper params.

And everything should be full deterministic here.

I have also tried setting np.random.seed again.

What might be happening is some shadow random state in the other libraries? 


And in case it is useful, here is a fixed version of the above also with the fixer code applied after the assert Fails. 

```python
import sys

import numpy as np
from scipy.linalg import cholesky, cho_solve
from sklearn.gaussian_process import GaussianProcessRegressor, kernels
from sklearn.gaussian_process.kernels import *

np.random.seed(0)

X = np.random.randn(100, 3)
y = np.random.randn(100, 1)

KERNEL = Matern()
# make it once, fit and get the params
gp_ = GaussianProcessRegressor(kernel=KERNEL, normalize_y=True, alpha=1e-5)
gp_.fit(X, y)
params = gp_.kernel_.get_params(deep=False)  # tried deep=True but doesn't matter

# make it again, don't fit, but set the params
gp = GaussianProcessRegressor(kernel=KERNEL, normalize_y=True, alpha=1e-5)
np.random.seed(0)
gp.optimizer = None  # this should make fit do very little
gp.fit(X, y)  # this creates kernel

# do it again for good measure? Maybe we need to do this several times to trigger it?
gp.kernel_.set_params(**params)
gp.fit(X, y)  # do this again to now recreate L_
gp.kernel_.set_params(**params)


def get_K_L_(self):
    K = self.kernel_(self.X_train_)
    K[np.diag_indices_from(K)] += self.alpha
    L_ = cholesky(K, lower=True)
    return K, L_


Ka, La = get_K_L_(gp_)
Kb, Lb = get_K_L_(gp)

a = gp_.__getstate__()
b = gp.__getstate__()

# here you will see alpha_ and L_ are not the same. Something is not updating them correctly?
# see https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L322
for k in a:
    if isinstance(a[k], np.ndarray):
        print(k, type(a[k]), np.all(a[k] == b[k]))


assert np.all(gp_.X_train_ == gp.X_train_)  # passes
assert gp_.alpha == gp.alpha  # passes
assert np.allclose(Ka, Kb), 'Ka Kb do not match'  # passes
assert np.allclose(La, Lb), 'La Lb do not match'  # passes
# it is bizarre that the output of the kernels matchs but not L_????
failed = False
try:
    assert np.allclose(gp.L_, gp_.L_), 'L_ do not match'  # FAILS
except Exception as e:
    print("yes L_ do not match as expected")
    failed = True
assert failed, 'should have failed'


# NOTE: this fixes things
def fix_things(self):
    K = self.kernel_(self.X_train_)
    K[np.diag_indices_from(K)] += self.alpha
    L_ = cholesky(K, lower=True)
    self.L_ = L_
    self.K_inv_ = None
    self.alpha_ = cho_solve((self.L_, True), self.y_train_)

fix_things(gp)

assert np.all(gp_.X_train_ == gp.X_train_)  # passes
assert gp_.alpha == gp.alpha  # passes
assert np.allclose(Ka, Kb), 'Ka Kb do not match'  # passes
assert np.allclose(La, Lb), 'La Lb do not match'  # passes
# it is bizarre that the output of the kernels matchs but not L_????
assert np.allclose(gp.L_, gp_.L_), 'L_ do not match'  # FAILS
print('everything matches after fixing')
```

With output

```shell
_y_train_mean <class 'numpy.ndarray'> True
_y_train_std <class 'numpy.ndarray'> True
X_train_ <class 'numpy.ndarray'> True
y_train_ <class 'numpy.ndarray'> True
L_ <class 'numpy.ndarray'> False
alpha_ <class 'numpy.ndarray'> False
yes L_ do not match as expected
everything matches after fixing
```

Labels: Bug, module:gaussian_process, Needs Investigation