# Marginalization of GP regression hyperparameters with Laplace approximation

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)?

• In the meanwhile, I tried Elliptical Slice Sampling, which is still much slower than optimization. I also implemented the method above with a numerical computation of the Hessian (as a proof of concept, before implementing an analytical computation of the Hessian). For the few test cases I tried, it doesn't seem to work noticeably better than bare optimization (but the difference might emerge for some specific hard cases; so I'd need more extensive benchmarking). – lacerbi Sep 21 '15 at 7:27
• I suspect the problem may be that the posterior distribution of the hyper-parameters is not very Gaussian, in which case the Laplace approximation is unlikely to work well. Do the results of slice sampling suggest a Gaussian is a reasonable assumption? IIRC Maurizio Fillipone (ETA e.g. eurecom.fr/~filippon/Publications/ieee-tpami14.pdf) has done some excellent work on this problem. I suspect the real value of marginalisation only shows up is when there are many hyper-parameters of the covariance to tune. Great question (+1). – Dikran Marsupial Sep 21 '15 at 7:36
• Thanks. When I have more than a hundred inputs, the Gaussian approximation seems okay (not for all cases, but in many cases). I found some work by Filippone but it's for estimating the posterior over hyper-parameters when the latent variables cannot be marginalized (e.g., in classification). My case is a bit peculiar because it's technically easy for the typical GP literature; but I want it to be very fast. Btw I share the same suspicion: so far I've only tested cases with a few hyper-parameters because the numerical computation of the Hessian is demanding; I should test with more. – lacerbi Sep 21 '15 at 7:58
• If there are only a couple of hyper-parameters, in my experience, just optimizing them works reasonably well as there is little scope for overfitting the evidence. It falls apart very quickly if you try to use a composite covariance function, or one with many hyper-parameters. This is a bit of a problem for GPs/kernel machines and marginalisation is probably the best approach. IIRC Fillipone has done work on marginalisation over the hyper-parameters; I was hoping to have a go at it myself but have lacked the time/energy. Be keen to hear how your experiments turn out. – Dikran Marsupial Sep 21 '15 at 8:45
• Thanks. At a quick glance, the paper you linked is about marginalizing hyper-parameters when the marginal likelihood is not analytical (which adds another layer of complication). For that, it uses a pseudo-marginal approach, and then pretty standard MH with importance sampling. Unfortunately, this doesn't help much with my problem (my marginal likelihood is already analytical). I'll keep looking at other Filippone's papers. – lacerbi Sep 21 '15 at 15:43

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.

• Thanks very much for the links. I downloaded/forked the gpml_extensions repository a few weeks ago and been using it for a while, it's a very nice extension (I wrote a couple of additional functions and found a couple of bugs in the original repository -- was planning to upload changes but I hadn't had time, will do). – lacerbi Oct 31 '15 at 4:01
• Also, the paper sounds interesting and I'll keep it in mind. However, for the moment for this project I dropped the marginalization part; it seems that [for my problem] marginalization of hyper-parameters yields little improvement wrt optimization; and the improvement doesn't offset the much higher computational cost. – lacerbi Oct 31 '15 at 4:09

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).