Basic SingleTaskGP fit question

[dbernstein-symbio]

Hi Everyone, I'm getting started with botorch and I'm playing with fitting a SingleTaskGP to the noisy sinusoidal data described here:

https://botorch.org/v/0.1.0/tutorials/fit_model_with_torch_optimizer

For a small number of sample points, like 10-20, the mean and confidence intervals look ok but for the same code with a large number of points, 200, 500, the confidence interval is too small. I've attached the code and its output plot. What am I doing wrong?

Thanks in advance,
Dave

Here is the code:
import torch
from botorch.models import SingleTaskGP
from botorch.fit import fit_gpytorch_model
from gpytorch.mlls import ExactMarginalLogLikelihood
from gpytorch.constraints import GreaterThan

import math
from matplotlib import pyplot as plt

N = 500
noise = 0.3

def test():
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(f'running on {device}')

dtype = torch.float

# use regular spaced points on the interval [0, 1] (training data needs an explicit dimension)
train_X = torch.linspace(0, 1, N, dtype=dtype, device=device).unsqueeze(1)

# sample observed values and add some synthetic noise
train_Y = torch.sin(train_X * (2 * math.pi)) + noise * torch.randn_like(train_X)
#train_Y = train_Y.view(-1)

gp = SingleTaskGP(train_X, train_Y)
gp.likelihood.noise_covar.register_constraint("raw_noise", GreaterThan(1e-5))

mll = ExactMarginalLogLikelihood(likelihood=gp.likelihood, model=gp)
mll.to(train_X)

fit_gpytorch_model(mll)

# Initialize plot
f, ax = plt.subplots(1, 1, figsize=(6, 4))

# test model on 101 regular spaced points on the interval [0, 1]
test_X = torch.linspace(0, 1, 101, dtype=dtype, device=device)

# no need for gradients
with torch.no_grad():
    # compute posterior
    posterior = gp.posterior(test_X)

    # Get upper and lower confidence bounds (2 standard deviations from the mean)
    lower, upper = posterior.mvn.confidence_region()

    # Plot training points as black stars
    ax.plot(train_X.cpu().numpy(), train_Y.cpu().numpy(), 'k*')

    # Plot posterior means as blue line
    ax.plot(test_X.cpu().numpy(), posterior.mean.cpu().numpy(), 'b')

    # Shade between the lower and upper confidence bounds
    ax.fill_between(test_X.cpu().numpy(), lower.cpu().numpy(), upper.cpu().numpy(), alpha=0.5)

ax.legend(['Observed Data', 'Mean', 'Confidence'])
plt.tight_layout()

plt.show()
plt.savefig('mynewfig.png')

if name == 'main':
test()

and here is a sample output:

[dbernstein-symbio]

I don't know why the code snippet came out that way.

[saitcakmak]

Hi Dave! What you're plotting is the mean and confidence interval for the underlying sinusoidal function. It doesn't include the predicted observation noise. To get the predictions with the observation noise included, you can use posterior = gp.posterior(test_X, observation_noise=True).

[dbernstein-symbio]

Omg, thank you so much!

[dbernstein-symbio]

Wondering if this should be added to the tutorial as well?

[saitcakmak]

Yeah, good call. Updating the tutorial in #1196

[dbernstein-symbio]

Thanks. Another fyi regarding this: in order to get the SGD training loop to run I had to modify the lines including loss so they used loss.mean(). Like this

for epoch in range(NUM_EPOCHS):
# clear gradients
optimizer.zero_grad()

    # forward pass through the model to obtain the output MultivariateNormal
    output = gp(train_X)

    # Compute negative marginal log likelihood
    loss = -mll(output, train_Y)

    # back prop gradients
    loss.mean().backward()

    # print every 10 iterations
    if (epoch + 1) % 10 == 0:
        print(
            f"Epoch {epoch+1:>3}/{NUM_EPOCHS} - Loss: {loss.mean().item():>4.3f} "
            f"lengthscale: {gp.covar_module.base_kernel.lengthscale.item():>4.3f} " 
            f"noise: {gp.likelihood.noise.item():>4.3f}" 
         )
    optimizer.step()