[Bug] Input normalization makes GP fail to fit data

[martinakaduc]

🐛 Bug of input normalization

I tried to fit GP on Ackley function with 200 points. The GP fit very well without input normalization, but after adding input normalization, the GP posterior surface became smoother, resulting in a convex surface.

Vanilla GP

GP + Input normalization

GP + Output standardization

GP + Input normalization + Output standardization

To reproduce

** Code snippet to reproduce **

# Re-import necessary packages and redefine required functions after environment reset
import torch
import numpy as np
import matplotlib.pyplot as plt
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_mll
from botorch.models.transforms.outcome import Standardize
from botorch.models.transforms.input import Normalize
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.kernels import MaternKernel, ScaleKernel
from gpytorch.constraints import Positive
import gpytorch

# Define the Ackley function again
def ackley(x, y):
    return -20 * np.exp(-0.2 * np.sqrt(0.5 * (x**2 + y**2))) - np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) + 20 + np.e

np.random.seed(42)
num_samples = 100

x_samples = np.random.uniform(-2, 2, num_samples)
y_samples = np.random.uniform(-2, 2, num_samples)
z_samples = ackley(x_samples, y_samples)
X_train_torch = torch.tensor(np.vstack((x_samples, y_samples)).T, dtype=torch.float32)
y_train_torch = torch.tensor(z_samples, dtype=torch.float32).unsqueeze(-1)

# model_botorch = SingleTaskGP(X_train_torch, y_train_torch)

model_botorch = SingleTaskGP(
    X_train_torch, 
    y_train_torch,
    input_transform=Normalize(d=2, bounds=torch.tensor([[-2, -2], [2, 2]])),
    outcome_transform=Standardize(1),
)
mll = ExactMarginalLogLikelihood(model_botorch.likelihood, model_botorch)
fit_gpytorch_mll(mll)

# Define the grid for predictions
x_new = np.linspace(-2, 2, 400)
y_new = np.linspace(-2, 2, 400)
X_new, Y_new = np.meshgrid(x_new, y_new)
X_new_torch = torch.tensor(np.c_[X_new.ravel(), Y_new.ravel()], dtype=torch.float32)

# Get posterior mean
with torch.no_grad():
    posterior = model_botorch.posterior(X_new_torch)
    mean = posterior.mean.detach().numpy()
    mean = mean.reshape(X_new.shape)

# Plot the GP regression results using Botorch
fig, ax = plt.subplots(figsize=(5, 4))
contour_botorch = ax.contourf(X_new, Y_new, mean, 20, cmap='viridis')
plt.colorbar(contour_botorch)
ax.scatter(x_samples, y_samples, color='red', label='Samples')
ax.set_title("GP Regression on Ackley Function (Botorch with Matern Kernel)")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.legend()
plt.show()

System information

Please complete the following information:

• BoTorch: 0.11.0
• GPyTorch: 1.11
• PyTorch: 2.2.1+cu121
• Ubuntu 20.04

[hvarfner]

Not a BoTorch dev, but from experience, Ackley is very tricky to fit.

TL;DR: Most likely not a bug, but the GP hyperparameter prior's interaction with input/output scaling.

Long reply:
The rapid variation in the function in conjunction with the hyperparameter prior generally results in an Ackley fit with a large amount of noise (the 2:nd plot).

If you don't normalize the input from[-2, 2] to [0, 1], you effectively suggest for the variation in the objective to be 4 times as rapid under the lengthscale prior. This in turn makes it more likely that the GP infers the variation in the function as opposed to interpreting it as noise.

If you don't normalize the output, a similar thing happens. There is more signal, so you are less likely to fit a high-noise GP.

My guess would be that the second plot (the "poor" fit) results in better optimization, since the regularization we get from the noise tends to aid optimization in the case of Ackley.

If you removed the prior on the hyperparameters, these examples_should all give you the same fit.