GaussianProcessRegressor doesn't work with multidemensional output when normalize_y=True #18065
Description
The GaussianProcessRegressor doesn't work with multidemensional output when normalize_y=True. In this example, the code runs fine when normalize_y is False, but breaks when it is true:
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
def f(x): return(np.array([np.sin(7 * x), x ** 4]))
kernel = RBF()
gp = GaussianProcessRegressor(kernel=RBF(length_scale=15.7), n_restarts_optimizer=50,
normalize_y=True) # (works when normalize_y is False)
X = np.linspace(0, 5, 5)
gp.fit(np.atleast_2d(X).T, f(X).T)
newx = np.atleast_2d([1, 2, 3, 4]).T
gp.predict(newx, return_std=False)
gp.predict(newx, return_std=True)