Fast Riemannian-manifold Hamiltonian Monte Carlo for hierarchical Gaussian-process models

TL;DR

This work advances efficient Bayesian inference for hierarchical GP models by combining a reduced-rank GP representation with RMHMC using a soft-absolute Hessian metric. By reorganizing computations to exploit model structure and implementing dynamically updated eigendecompositions, the approach achieves substantial speedups over general-purpose libraries and enables posterior sampling and model evidence calculation beyond Laplace approximations. The method is validated on Bayesian logistic regression and real-world NMES data, showing improved posterior exploration and the ability to compare nonlinear mean/variance GP models via BME. The study highlights the need for customisable, structure-aware libraries in GP-based inference to enable scalable, unbiased analysis of complex real-world datasets.

Abstract

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 optimising the computation order according to the model structure and dynamically programming the eigendecomposition. This improvement cannot be achieved when using an existing library based on a naive automatic differentiator. 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 finely optimises the mode of automatic differentiation depending on the model structure.

Figures (2)

• Figure 1: RMHMC for Bayesian logistic regression on simulated data. Comparison between our implementations based on the dynamic and static eigendecomposition of the Hessian and the implementation based on Hamiltorch. (A) Representative trajectories of log-posterior density for $D=1$ ($d=34$) and $D=16$ ($d=484$). (B) Mean and standard deviation of the wall time spent on computing $100$ blocks of $100$ leapfrogs are shown on the linear (main panel) and logarithmic (inset) scales for different model sizes. (C) Mean and standard deviation of the number of sweeps in the cyclic Jacobi method for a single dynamic and static eigendecomposition of the metric are shown on the logarithmic scale for different model sizes. In (B) and (C), model sizes are shown on the logarithmic scales. (Abbreviations) prop.(dyn.): the proposed method based on dynamic eigendecomposition, prop.(syevd): the proposed method based on static eigendecomposition with the divide-and-conquer algorithm, prop.(syevj): the proposed method based on static eigendecomposition with the Jacobi method.
• Figure 2: Comparison of proposed implementation of RMHMC and NUT-HMC in Hamiltorch. (A) Representative trajectories of the log-posterior density for $D=1$ ($d=34$) obtained using the proposed implementation of RMHMC based on dynamic eigendecomposition (prop.(dyn)) and NUT-HMC in Hamiltorch from the same initial condition $f=0$ and $c_g, \sigma _g, c_\ell =1$. A representative trajectory for the proposed implementation of RMHMC starting with the parameter values obtained after 2000 seconds of sampling with NUT-HMC is also shown. (B) Eigenvalues of the Hessians of the negative log-posterior density for the samples of parameter values obtained after running the proposed implementation of RMHMC and NUT-HMC in Hamiltorch for 2,000 seconds. The eigenvalues are plotted in descending order on the logarithmic scale. For negative eigenvalues, their absolute values are plotted. The absolute values of the 32nd to 34th eigenvalues for the sample from NUT-HMC were extremely small and are not shown in the plot.