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