Lengthscale parameters of the posterior distribution

[madhavkrishnan94]

Hello!

I am working with a multi-input, multi-output dataset.

So far, I have been using a SingleTaskGP model to fit the data, which gives me the option to play around with the hyperparameters such as the concentration and rate of the Gamma distribution of the lengthscale prior. Is there a way to obtain the lengthscale distribution of the posterior after the model training has happened?

Meanwhile, I was reading up on Fully Bayesian GP, which learns a posterior distribution for the hyperparameters using the No-U-Turn-Sampler. The documentation also says that

our “standard” models (e.g. SingleTaskGP) learn only a single best value for each hyperparameter using MAP.

However, the SaasFullyBayesianSingleTaskGP requires the train_Y to be of dimensions n x 1. This is not particularly useful for me as I require train_Y to support multi-ouput (n x d in other words).

I'd like to clarify that I am not particular on using a Fully Bayesian GP, in my attempt to get the posterior distribution of the lengthscale. Any leads from the community on this will be greatly helpful.

Thanks.

[saitcakmak]

Hi @madhavkrishnan94. As you found out, the SingleTaskGP model has a single lengthscale rather than a distribution over lengthscales. You can access the lengthscale of the model by checking the attributes of its covariance kernel: model.covar_module.base_kernel.lengthscale.

However, the SaasFullyBayesianSingleTaskGP requires the train_Y to be of dimensions n x 1. This is not particularly useful for me as I require train_Y to support multi-ouput (n x d in other words).

We have a ModelListGP that allows you to combine multiple GP models together into a multi-output model. For the SAAS model, you can construct one model for each of your outputs and combine them using a ModelListGP to obtain a multi-output SAAS model.

Note that the SAAS model does not have a proper distribution over lengthscales either. It is an ensemble model where each model in this ensemble has a different lengthscale sampled from the posterior distribution using NUTS. By default, there are 16 models on the ensemble.

[saitcakmak]

The design of priors is generally a tricky task. You can design priors for a specific task based on domain knowledge, though it is much harder to design priors that work well in a large class of problems. If we're trying out a new prior, we would typically benchmark them against the old priors on a large number of problems to see how they perform. If you have intuition / domain knowledge on how smooth you expect the function to be, you can plot some draws from the model prior with different lenghtscale values, see how the plots look like based on how smooth you expect the function to be, and set a likelihood that assigns higher probability to the lengthscales you think are more likely.

FWIW, even a simple GP with no priors works quite well surprisingly often. Most Gamma priors will assign very small probability to very large lengthscales, which can be harmful in high dimensional search spaces. The SAAS model / prior mentioned above does the opposite of this and works quite well in high dimensions.
cc @dme65 in case you have any suggestions on designing priors

[madhavkrishnan94]

Thanks a lot for your insights, @saitcakmak

Although I did think about it earlier, I am currently not using a SAAS model in my work.

Here's another thing on the same note:

After I define the SingleTaskGP model (like I showed in my previous reply) and before I fit the model using fit_gpytorch_model(mll), I print some values like

print("\nLengthscale before: ", model.covar_module.base_kernel.lengthscale)
print("\nOutputscale before: ",model.covar_module.outputscale)
print("\nNoise before: ", model.likelihood.noise)

The output reads as

Lengthscale before: tensor([[0.6931, 0.6931, 0.6931]], device='cuda:0', dtype=torch.float64,
grad_fn=)

Outputscale before: tensor(0.6931, device='cuda:0', dtype=torch.float64,
grad_fn=)

Noise before: Parameter containing:
tensor([[2.0000],
[2.0000],
[2.0000],
[2.0000]], device='cuda:0', dtype=torch.float64, requires_grad=True)

I have two questions here:

1. The lengthscale and outputscale values that are printed- how are they calculated from their corresponding gamma priors (because that is how I initialized them inside the GP function)?

2. After fitting the model, when I print the same, are the values that get displayed calculated in the same way from a distribution of the posterior? Is there a way for me to print/plot this distribution?

Thanks!

[saitcakmak]

Hi @madhavkrishnan94. These questions are going deep into GPyTorch territory and reaching the limits of my understanding. For more detailed answers, you may want to ask in the GPyTorch repo or read through the code to see how things work.

At a high level, the values you see when you print lenghtscale/outputscale/noise are transformed values. GPyTorch stores a raw_lengthscale value and transforms this to get the lengthscale. See: https://github.com/cornellius-gp/gpytorch/blob/master/gpytorch/kernels/kernel.py#L293. The raw_lengthscale seems to be initialized to 0, and you end up with the value shown after transforming. I believe either constraint or the prior can be used to specify an initial value for these, though I am not fully sure how that works.

The values you see after fitting are similarly transformed values of the trained raw_lengthscale etc. They're singular values and not distributions (i.e., this is not fully-Bayesian inference).

[madhavkrishnan94]

Thanks for your reply.

It was quite helpful, I went through the GPyTorch forum and confirmed a few things.

The initialization always sets the raw_lengthscale and other "raw" quantities to zero. Since a log(1+exp(x)) transform is applied (to ensure positivity constraint), the default lengthscale we see before the fitting is the transformed value of 0.6931.

The other option is to set the values of the lenghtscale/outputscale manually.

Coming to the role of the lengthscale prior:

The log likelihood of the prior wrt the lengthscale is added into the log marginal likelihood during the training process.

This still begs the question. How do I arrive at the right choice of arguments to the lengthscale_prior, in my case the GammaPrior(alpha, beta), in the case of me not having too much knowledge about the problem beforehand to know what the typical lengthscales will be?

[saitcakmak]

How do I arrive at the right choice of arguments to the lengthscale_prior, in my case the GammaPrior(alpha, beta), in the case of me not having too much knowledge about the problem beforehand to know what the typical lengthscales will be?

The short answer is you don't. The whole point of the prior is to inject prior domain knowledge into the model. If you don't have any prior knowledge on what to expect from the problem, you can't design a prior for it. The best you can do is to benchmark several priors on some related problems and pick the best performing one.