Plotting prediction

[stbecet]

My question is about the application of a GP model, lets say SingleTaskGP. If i choose 20 random points for train_x and 20 random points for train_y I get a plot like this:

The solid line is the current prediction of the model I initialized on the 20 random points and plot with np.linspace(0.0, 9.0, 100). The blue area represents the upper and lower confidence bound. The '*' are the known observations or in my case the 20 random points.

Why is the prediction so far from the observed dots? I observed this phenomenon couple of times. Is there anything I should watch? Furthermore the upper and lower uncertainty is very high at each point. Im wondering why..

Hope you can help.

PS the code I used to generate the plot ist below:

import torch
import numpy as np
import matplotlib.pyplot as plt

from botorch.fit import fit_gpytorch_model
from gpytorch.mlls import ExactMarginalLogLikelihood
from botorch.models import SingleTaskGP


train_x = (torch.rand(20)*9.0).unsqueeze(-1)
train_y = (torch.rand(20)*8.0).unsqueeze(-1)

m = SingleTaskGP(train_x, train_y)
mll = ExactMarginalLogLikelihood(m.likelihood, m)
fit_gpytorch_model(mll);

test_x = torch.linspace(0.0, 9.0, 100)
preds = m(test_x)
lower, upper = preds.confidence_region()

plt.plot(np.linspace(0.0, 9.0, 100), preds.mean.detach().numpy())
plt.fill_between(test_x.detach().numpy(), lower.detach().numpy(), 
                 upper.detach().numpy(), alpha=0.5)
plt.plot(train_x, train_y, 'k*')
plt.grid()

[saitcakmak]

Hi @stbecet. What you're plotting above is the predictive (posterior) distribution of the true function, without the (inferred) observation noise. SingleTaskGP model assumes that the observations are with homoscedastic Gaussian noise with unknown noise level, and infers the noise level internally during model fit. When you call posterior on the trained model, either with m(test_x) as above or with m.posterior(test_x), you get the posterior on the function without this added noise. When training a model on torch.rand, the inferred noise level ends up being pretty large (you can check this in m.likelihood.noise), so you get a plot like this.
To add the observation noise to the predictions, you need to pass your posterior through the likelihood, i.e., after preds = m(test_x) call preds_w_noise = m.likelihood(preds). This should result in a plot with wider confidence region that covers your training points in most cases. If you instead used m.posterior(test_x) to get the posterior, you could do this by passing the argument m.posterior(test_x, observation_noise=True). Hope this clarifies it!

[stbecet]

Hi @saitcakmak. Thanks a lot for your reply. So the high uncertainty at the observed points results from homoskedasticity assumption of the SingleTaskGP model I guess. I don't want the uncertainty be higher to cover the observations but be lower at the observed points (and the mean posterior or prediction itself be near to these points). Should the application of the model HeteroskedasticSingleTaskGP instead of SingleTaskGP solve my problem?

Anyway I don't understand, how the model is adding the observation noise. Isn't the assumption that the observed point already includes some noise?

[stbecet]

PS: here is another plot I made with SingleTaskGP. This is the desired behavior I meant before: near the observations there is low uncertainty and the prediction of the posterior is near to these points:

In this plot I calculated the predictions with model(test_x). Wondering why the plot looks so different (and in my view more correct).

But sometimes I get with the same code I used for the previous plot something like this:

[Balandat]

The reason you're likely seeing this behavior is that the models by default use MAP fitting - they maximize the likelihood of the model parameters under the observed data. If there is very little data, this can cause issues as that likelihood may end up not having an obvious single global maximum. Since optimizing the likelihood means solving a non-convex optimization problem, the optimizer can get stuck in local minima (corresponding to the "undesirable" behavior). In the second plot above the model seems to be estimating large lengthscales, resulting in essentially predicting the mean of the observations, whereas in the first plot you end up with short lengthscales and very local behavior.

Some remedies:

1. Generally, more data will result in the likelihood to be more peaked and reduce the chance for this. But of course in practice you can't always add on more data
2. To avoid getting "stuck" in local maxima, one could re-start the optimization from a number of different randomized initial conditions on the hyperparameters and then take the best. We currently don't do this by default, since we haven't found it to be necessary in general (we may want to make this a standard option going forward though). You can use sample_all_priors from botorch.optim.utils to draw random samples of the hyperparameters from their priors as new starting conditions for the optimization and then call fit_gpytorch_model on the mll.
3. If you really are in the very low data regime, then MAP estimation may be problematic more generally, and you may want to consider fully Bayesian inference as discussed in this gpytorch tutorial. We have used this in botorch as well and have some tutorials in the works that we hope to put out soon.

[stbecet]

Thank you for that @Balandat!
I see. Understanding the complex maths behind the library is not thaaaat wrong :'D