Marginalization of GP regression hyperparameters with Laplace approximation

Question

I am using Gaussian Processes (GP) for regression (via the gpml package for MATLAB).

So far, I was optimizing the hyper-parameters by maximizing the log likelihood, but I would like to try a more Bayesian approach by (approximately) marginalizing over the hyper-parameters.

The hard part is that I would also like the training to remain relatively fast. At the moment the optimization is very fast since the likelihood is Gaussian, so everything is analytical (likelihood and gradient). Therefore, maximizing the likelihood with pseudo-Newton methods requires only a handful of function evaluations.

I used MCMC a lot for other projects (e.g., slice sampling) and in my experience it tends to be either slow (i.e., getting independent samples might require a lot of time) and/or it requires a lot of careful tuning and checking.

As a middle-ground method (in-between optimization and proper marginalization), I was thinking of using a Laplace approximation combined with importance sampling. The algorithm would be:

• maximize the log likelihood just as I am doing now;
• compute the Hessian analytically at the maximum (as far as I know gpml doesn't do it, but it shouldn't be too hard to write it down);
• sample parameters from the multivariate normal distribution centered on the mode and with covariance matrix equal to the inverse Hessian;
• give weights to the samples via importance sampling.

On paper, this method should work pretty quickly (the hyper-parameter space is relatively small so I don't see major problems in computing/inverting the Hessian). Of course if the posterior is not Gaussian it's not as good as a proper marginalization, but the importance sampling bit should account for some deviations from normality.

Do you foresee any problem with this? Alternatively, is there a sampling method which is (almost) as fast as optimization (and with no tuning required)?

Answer 1

You might be interested in the gpml_extensions repository here:

https://github.com/rmgarnett/gpml_extensions/

There is code for computing the Hessian of the log likelihood wrt $\theta$ for both exact inference and the Laplace approximation for approximate GP inference. There is also convenience code for using these to find the Laplace approximation to the hyperparameter posterior $p(\theta \mid \mathbf{X}, \mathbf{y})$ (theta_posterior_laplace.m).

Finally, this paper from UAI 2014 suggests a fast analytical approximation (called the MGP) to the posterior predictive distribution

$p(f^\ast \mid \mathbf{x}^\ast, \mathbf{X}, \mathbf{y}) = \int p(f^\ast \mid \mathbf{x}^\ast, \mathbf{X}, \mathbf{y}, \theta) p(\theta \mid \mathbf{X}, \mathbf{y}) \, \mathrm{d}\theta$,

under an arbitrary Gaussian approximation to the $\theta$ posterior:

$p(\theta \mid \mathbf{X}, \mathbf{y}) \approx \mathcal{N}(\theta; \hat{\theta}, \mathbf{\Sigma})$.

The Laplace approximation would be one way to derive such an approximation.

There is an implementation of the MGP built atop gpml/gpml_extensions available from the same user, but I don't have the reputation to post that link.

Answer 2

The best reference I could find online so far, and a very fitting one, is Ville Pietiläinen's MSc thesis: Approximations for Integration Over the Hyperparameters in Gaussian Processes (2010).

Pietiläinen compares a point estimate approach (so-called MAP-II, since the latent variables are marginalized analytically given a Gaussian likelihood) to three marginalization methods:

1. Grid search (with a "smart" grid)
2. Central composite design (CCD)
3. Quasi-random importance sampling with Student-$t$ proposal distribution

Interestingly, for the case studies considered in the thesis, marginalization of the hyperparameters seems to provide some benefits with respect to optimization only for small datasets (e.g., no more than 50-100 training points). More in general, as Pietiläinen argues, the advantage of marginalization should emerge for test points in regions in which the input density is low (hence prediction uncertainty is high).

The three different marginalization methods perform somewhat similarly, although the comparison in the thesis is not exhaustive. CCD is appealing since it requires many less function evaluations than the other methods.

Another interesting reference is Philip Boyle's PhD thesis: Gaussian Processes for Regression and Optimisation (2007). In particular, chapter 8 is focussed on marginalization over hyperparameters.

In conclusion, I will eventually try and implement the CCD method, although at this point is not a priority as I am not expecting a major gain over optimization. I can probably better spend my time by playing with other factors that have a larger impact on the quality of the prediction (e.g., the choice of covariance function).