Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Randomized Runge-Kutta-Nyström Methods for Unadjusted Hamiltonian and Kinetic Langevin Monte Carlo

Nawaf Bou-Rabee, Tore Selland Kleppe

TL;DR

This work tackles efficient sampling from high-dimensional targets using unadjusted Hamiltonian and kinetic Langevin kernels, where discretization introduces bias that scales with computational cost. It introduces randomized Runge-Kutta-Nyström schemes of order $5/2$ and $7/2$ that exploit the second-order structure and a triangular random variable to achieve higher $L^2$-accuracy under gradient Lipschitz constant $L$ and Hessian Lipschitz constant $L_H$. Theoretical contributions provide quantitative $L^2$-accuracy bounds for the $2.5$-order scheme (and analogous bounds for the kinetic Langevin variant), along with moment stability and local error lemmas; these underpin a rigorous strong convergence analysis. Numerical experiments demonstrate substantial reductions in bias per gradient evaluation and overall computational cost compared with Verlet and stratified Monte Carlo across well-behaved targets, highlighting the practical impact for high-dimensional MCMC.

Abstract

We introduce $5/2$- and $7/2$-order $L^2$-accurate randomized Runge-Kutta-Nyström methods, tailored for approximating Hamiltonian flows within non-reversible Markov chain Monte Carlo samplers, such as unadjusted Hamiltonian Monte Carlo and unadjusted kinetic Langevin Monte Carlo. We establish quantitative $5/2$-order $L^2$-accuracy upper bounds under gradient and Hessian Lipschitz assumptions on the potential energy function. The numerical experiments demonstrate the superior efficiency of the proposed unadjusted samplers on a variety of well-behaved, high-dimensional target distributions.

Randomized Runge-Kutta-Nyström Methods for Unadjusted Hamiltonian and Kinetic Langevin Monte Carlo

TL;DR

This work tackles efficient sampling from high-dimensional targets using unadjusted Hamiltonian and kinetic Langevin kernels, where discretization introduces bias that scales with computational cost. It introduces randomized Runge-Kutta-Nyström schemes of order and that exploit the second-order structure and a triangular random variable to achieve higher -accuracy under gradient Lipschitz constant and Hessian Lipschitz constant . Theoretical contributions provide quantitative -accuracy bounds for the -order scheme (and analogous bounds for the kinetic Langevin variant), along with moment stability and local error lemmas; these underpin a rigorous strong convergence analysis. Numerical experiments demonstrate substantial reductions in bias per gradient evaluation and overall computational cost compared with Verlet and stratified Monte Carlo across well-behaved targets, highlighting the practical impact for high-dimensional MCMC.

Abstract

We introduce - and -order -accurate randomized Runge-Kutta-Nyström methods, tailored for approximating Hamiltonian flows within non-reversible Markov chain Monte Carlo samplers, such as unadjusted Hamiltonian Monte Carlo and unadjusted kinetic Langevin Monte Carlo. We establish quantitative -order -accuracy upper bounds under gradient and Hessian Lipschitz assumptions on the potential energy function. The numerical experiments demonstrate the superior efficiency of the proposed unadjusted samplers on a variety of well-behaved, high-dimensional target distributions.
Paper Structure (11 sections, 8 theorems, 86 equations, 4 figures)

This paper contains 11 sections, 8 theorems, 86 equations, 4 figures.

Table of Contents

  1. Introduction
  2. The Randomized Algorithms
  3. Randomized Runge-Kutta-Nyström 2.5 & 3.5
  4. Application to an Unadjusted Hamiltonian Monte Carlo Chain
  5. Application to an Unadjusted Kinetic Langevin Chain
  6. Numerical Illustrations
  7. Theory
  8. L2-accuracy Theorems
  9. Lemmas used in proofs of Theorems
  10. Proofs of L2-accuracy Theorems
  11. Proofs of Lemmas

Key Result

Theorem 1

Suppose that B1-B3 hold. Let $T>0$ satisfy $L T^2 \le 1/4$ and let $h>0$ satisfy $T/h \in \mathbb{Z}$. Then for any $x, v \in \mathbb{R}^d$ Here, $(q_{t_k}(x,v), v_{t_k}(x,v))$ and $(\tilde{Q}_{t_k}(x,v), \tilde{V}_{t_k}(x,v))$ are defined by exact and eq:rRKN2pt5, respectively.

Figures (4)

  • Figure 1: 95% confidence intervals for bias of (\ref{['numeric_estimator']}) associated with uHMC sampling with different time integration schemes and selected number of time steps $T/h$. The different panels correspond to the different target distributions, and the number of simulation replica $R$ specific to each target distribution are indicated. The confidence intervals are plotted according to the number of gradient evaluations $T/h$ required per transition of the $\mathcal{X}_k$ chain on the horizontal axis (log-scale).
  • Figure 2: 95% confidence intervals for bias of (\ref{['numeric_estimator']}) associated with ukLa sampling with different time integration schemes and selected number of time steps $T/h$. The different panels correspond to the different target distributions, and the number of simulation replica $R$ specific to each target distribution are indicated. The confidence intervals are plotted according to the number of gradient evaluations $T/h$ required per transition of the $\mathcal{X}_k$ chain on the horizontal axis (log-scale).
  • Figure 3: Confidence intervals for the $n \rightarrow \infty$ limit of the posterior standard deviation of the intercept term of a logistic regression model applied to the Pima data set. All results are based on $R=1000$ simulation replica with $n=1000$, $T=\pi/2$ and 100 transitions of burn-in. The gray shaded region indicates a 95% confidence interval calculated from the combined simulations with 64 (63 in the case of 3.5 scheme) gradient evaluations per transition.
  • Figure 4: $L^2$-Accuracy Verification. Left Image: A plot of the $L^2$-error in $(x,v)$-space of the rRKN time integrator for the linear oscillator with Hamiltonian $H(x,v) = (1/2) (v^2 + x^2)$. Right Image: A plot of the $L^2$-error in $(x,v)$-space of the rRKN time integrator for a double-well system with Hamiltonian $H(x,v) = (1/2) (v^2 + (1-x^2)^2)$. Both simulations have initial condition $(1,1)$ and unit duration. The time step sizes tested are $2^{-n}$ where $n$ is given on the horizontal axis. The dashed curve is $2^{-5 n/2} = h^{5/2}$ versus $n$. Corresponding results for the ukLa chain (cf. Theorem \ref{['thm:ukLa_L2_accuracy']}) are shown (in red x-marks) with friction factor $\gamma=1/4$.

Theorems & Definitions (18)

  • Remark 1
  • Theorem 1: $L^2$-accuracy of rRKN 2.5 w.r.t. Hamiltonian Flow
  • Remark 2
  • Remark 3
  • Theorem 2: $L^2$-accuracy of rRKN 2.5 based ukLa w.r.t. Exact Splitting
  • Remark 4
  • Lemma 3: Discrete Grönwall inequality
  • Lemma 4: A priori bounds for rRKN 2.5
  • Lemma 5: Moment bounds for ukLa with rRKN 2.5
  • Lemma 6: Lipschitz continuity of exact flow
  • ...and 8 more