Fast Riemannian Manifold Hamiltonian Monte Carlo for Hierarchical Gaussian Process Models

Hierarchical Bayesian models based on Gaussian processes are considered useful for describing complex nonlinear statistical dependencies among variables in real-world data. However, effective Monte Carlo algorithms for inference with these models have not yet been established, except for several simple cases. In this study, we show that, compared with the slow inference achieved with existing program libraries, the performance of Riemannian manifold Hamiltonian Monte Carlo (RMHMC) can be drastically improved by applying the chain rule for the differentiation of the Hamiltonian in the optimal order determined by the model structure, and by dynamically programming the eigendecomposition of the Riemannian metric with the recursive update of the eigenvectors at the previous move. This improvement cannot be achieved when using a naive automatic differentiator included in commonly used libraries. We numerically demonstrate that RMHMC effectively samples from the posterior, allowing the calculation of model evidence, in a Bayesian logistic regression on simulated data and in the estimation of propensity functions for the American national medical expenditure data using several Bayesian multiple-kernel models. These results lay a foundation for implementing effective Monte Carlo algorithms for analysing real-world data with Gaussian processes, and highlight the need to develop a customisable library set that allows users to incorporate dynamically programmed objects and to finely optimise the mode of automatic differentiation depending on the model structure.