Nawaf Bou-Rabee

Nawaf Bou-Rabee, Jerrold E. Marsden

In this paper structure-preserving time-integrators for rigid body-type mechanical systems are derived from a discrete Hamilton-Pontryagin variational principle. From this principle one can derive a novel class of variational partitioned Runge-Kutta methods on Lie groups. Included among these integrators are generalizations of symplectic Euler and Störmer-Verlet integrators from flat spaces to Lie groups. Because of their variational design, these integrators preserve a discrete momentum map (in the presence of symmetry) and a symplectic form. In a companion paper, we perform a numerical analysis of these methods and report on numerical experiments on the rigid body and chaotic dynamics of an underwater vehicle. The numerics reveal that these variational integrators possess structure-preserving properties that methods designed to preserve momentum (using the coadjoint action of the Lie group) and energy (for example, by projection) lack.

Nawaf Bou-Rabee, Martin Hairer, Eric Vanden-Eijnden

The Metropolis-Adjusted Langevin Algorithm (MALA), originally introduced to sample exactly the invariant measure of certain stochastic differential equations (SDE) on infinitely long time intervals, can also be used to approximate pathwise the solution of these SDEs on finite time intervals. However, when applied to an SDE with a nonglobally Lipschitz drift coefficient, the algorithm may not have a spectral gap even when the SDE does. This paper reconciles MALA's lack of a spectral gap with its ergodicity to the invariant measure of the SDE and finite time accuracy. In particular, the paper shows that its convergence to equilibrium happens at exponential rate up to terms exponentially small in time-stepsize. This quantification relies on MALA's ability to exactly preserve the SDE's invariant measure and accurately represent the SDE's transition probability on finite time intervals.

Nawaf Bou-Rabee, Houman Owhadi

This paper presents a Lie-Trotter splitting for inertial Langevin equations (Geometric Langevin Algorithm) and analyzes its long-time statistical properties. The splitting is defined as a composition of a variational integrator with an Ornstein-Uhlenbeck flow. Assuming the exact solution and the splitting are geometrically ergodic, the paper proves the discrete invariant measure of the splitting approximates the invariant measure of inertial Langevin to within the accuracy of the variational integrator in representing the Hamiltonian. In particular, if the variational integrator admits no energy error, then the method samples the invariant measure of inertial Langevin without error. Numerical validation is provided using explicit variational integrators with first, second, and fourth order accuracy.

Nawaf Bou-Rabee, Jesús María Sanz-Serna

This paper surveys in detail the relations between numerical integration and the Hamiltonian (or hybrid) Monte Carlo method (HMC). Since the computational cost of HMC mainly lies in the numerical integrations, these should be performed as efficiently as possible. However, HMC requires methods that have the geometric properties of being volume-preserving and reversible, and this limits the number of integrators that may be used. On the other hand, these geometric properties have important quantitative implications on the integration error, which in turn have an impact on the acceptance rate of the proposal. While at present the velocity Verlet algorithm is the method of choice for good reasons, we argue that Verlet can be improved upon. We also discuss in detail the behavior of HMC as the dimensionality of the target distribution increases.

Nawaf Bou-Rabee, Eric Vanden-Eijnden

This paper introduces time-continuous numerical schemes to simulate stochastic differential equations (SDEs) arising in mathematical finance, population dynamics, chemical kinetics, epidemiology, biophysics, and polymeric fluids. These schemes are obtained by spatially discretizing the Kolmogorov equation associated with the SDE in such a way that the resulting semi-discrete equation generates a Markov jump process that can be realized exactly using a Monte Carlo method. In this construction the spatial increment of the approximation can be bounded uniformly in space, which guarantees that the schemes are numerically stable for both finite and long time simulation of SDEs. By directly analyzing the generator of the approximation, we prove that the approximation has a sharp stochastic Lyapunov function when applied to an SDE with a drift field that is locally Lipschitz continuous and weakly dissipative. We use this stochastic Lyapunov function to extend a local semimartingale representation of the approximation. This extension permits to analyze the complexity of the approximation. Using the theory of semigroups of linear operators on Banach spaces, we show that the approximation is (weakly) accurate in representing finite and infinite-time statistics, with an order of accuracy identical to that of its generator. The proofs are carried out in the context of both fixed and variable spatial step sizes. Theoretical and numerical studies confirm these statements, and provide evidence that these schemes have several advantages over standard methods based on time-discretization. In particular, they are accurate, eliminate nonphysical moves in simulating SDEs with boundaries (or confined domains), prevent exploding trajectories from occurring when simulating stiff SDEs, and solve first exit problems without time-interpolation errors.

Nawaf Bou-Rabee, Houman Owhadi

Stochastic variational integrators for constrained, stochastic mechanical systems are developed in this paper. The main results of the paper are twofold: an equivalence is established between a stochastic Hamilton-Pontryagin (HP) principle in generalized coordinates and constrained coordinates via Lagrange multipliers, and variational partitioned Runge-Kutta (VPRK) integrators are extended to this class of systems. Among these integrators are first and second-order strongly convergent RATTLE-type integrators. We prove order of accuracy of the methods provided. The paper also reviews the deterministic treatment of VPRK integrators from the HP viewpoint.

Nawaf Bou-Rabee, Eric Vanden-Eijnden

This paper proposes a simple strategy to simulate stochastic differential equations (SDE) arising in constant temperature molecular dynamics. The main idea is to patch an explicit integrator with Metropolis accept or reject steps. The resulting `Metropolized integrator' preserves the SDE's equilibrium distribution and is pathwise accurate on finite time intervals. As a corollary the integrator can be used to estimate finite-time dynamical properties along an infinitely long solution. The paper explains how to implement the patch (even in the presence of multiple-time-stepsizes and holonomic constraints), how it scales with system size, and how much overhead it requires. We test the integrator on a Lennard-Jones cluster of particles and `dumbbells' at constant temperature.

Nawaf Bou-Rabee, Giulia Ortolan, Alessandro Saccon

This paper presents a counterexample to the conjecture that the semi-explicit Lie-Newmark algorithm is variational. As a consequence the Lie-Newmark method is not well-suited for long-time simulation of rigid body-type mechanical systems. The counterexample consists of a single rigid body in a static potential field, and can serve as a test of the variational nature of other rigid-body integrators.

Massimo Solinas, Agnes Valenti, Nawaf Bou-Rabee et al.

Scientific computing has long relied on double precision (64-bit floating point) arithmetic to guarantee accuracy in simulations of real-world phenomena. However, the growing availability of hardware accelerators such as Graphics Processing Units (GPUs) has made low-precision formats attractive due to their superior performance, reduced memory footprint, and improved energy efficiency. In this work, we investigate the role of mixed-precision arithmetic in neural-network based Variational Monte Carlo (VMC), a widely used method for solving computationally otherwise intractable quantum many-body systems. We first derive general analytical bounds on the error introduced by reduced precision on Metropolis-Hastings MCMC, and then empirically validate these bounds on the use-case of VMC. We demonstrate that significant portions of the algorithm, in particular, sampling the quantum state, can be executed in half precision without loss of accuracy. More broadly, this work provides a theoretical framework to assess the applicability of mixed-precision arithmetic in machine-learning approaches that rely on MCMC sampling. In the context of VMC, we additionally demonstrate the practical effectiveness of mixed-precision strategies, enabling more scalable and energy-efficient simulations of quantum many-body systems.

Nawaf Bou-Rabee, Siddharth Mitra, Andre Wibisono

Hamiltonian Monte Carlo (HMC) algorithms are among the most widely used sampling methods in high dimensional settings, yet their convergence properties are poorly understood in divergences that quantify relative density mismatch, such as Kullback-Leibler (KL) and Rényi divergences. These divergences naturally govern acceptance probabilities and warm-start requirements for Metropolis-adjusted Markov chains. In this work, we develop a framework for upgrading Wasserstein convergence guarantees for unadjusted Hamiltonian Monte Carlo (uHMC) to guarantees in tail-sensitive KL and Rényi divergences. Our approach is based on one-shot couplings, which we use to establish a regularization property of the uHMC transition kernel. This regularization allows Wasserstein-2 mixing-time and asymptotic bias bounds to be lifted to KL divergence, and analogous Orlicz-Wasserstein bounds to be lifted to Rényi divergence, paralleling earlier work of Bou-Rabee and Eberle (2023) that upgrade Wasserstein-1 bounds to total variation distance via kernel smoothing. As a consequence, our results provide quantitative control of relative density mismatch, clarify the role of discretization bias in strong divergences, and yield principled guarantees relevant both for unadjusted sampling and for generating warm starts for Metropolis-adjusted Markov chains.

Nawaf Bou-Rabee, Andreas Eberle

We provide quantitative upper bounds on the total variation mixing time of the Markov chain corresponding to the unadjusted Hamiltonian Monte Carlo (uHMC) algorithm. For two general classes of models and fixed time discretization step size $h$, the mixing time is shown to depend only logarithmically on the dimension. Moreover, we provide quantitative upper bounds on the total variation distance between the invariant measure of the uHMC chain and the true target measure. As a consequence, we show that an $\varepsilon$-accurate approximation of the target distribution $μ$ in total variation distance can be achieved by uHMC for a broad class of models with $O\left(d^{3/4}\varepsilon^{-1/2}\log (d/\varepsilon )\right)$ gradient evaluations, and for mean field models with weak interactions with $O\left(d^{1/2}\varepsilon^{-1/2}\log (d/\varepsilon )\right)$ gradient evaluations. The proofs are based on the construction of successful couplings for uHMC that realize the upper bounds.

Nawaf Bou-Rabee, Andreas Eberle

Andersen dynamics is a standard method for molecular simulations, and a precursor of the Hamiltonian Monte Carlo algorithm used in MCMC inference. The stochastic process corresponding to Andersen dynamics is a PDMP (piecewise deterministic Markov process) that iterates between Hamiltonian flows and velocity randomizations of randomly selected particles. Both from the viewpoint of molecular dynamics and MCMC inference, a basic question is to understand the convergence to equilibrium of this PDMP particularly in high dimension. Here we present couplings to obtain sharp convergence bounds in the Wasserstein sense that do not require global convexity of the underlying potential energy.

Nawaf Bou-Rabee, Andreas Eberle, Raphael Zimmer

Based on a new coupling approach, we prove that the transition step of the Hamiltonian Monte Carlo algorithm is contractive w.r.t. a carefully designed Kantorovich (L1 Wasserstein) distance. The lower bound for the contraction rate is explicit. Global convexity of the potential is not required, and thus multimodal target distributions are included. Explicit quantitative bounds for the number of steps required to approximate the stationary distribution up to a given error are a direct consequence of contractivity. These bounds show that HMC can overcome diffusive behaviour if the duration of the Hamiltonian dynamics is adjusted appropriately.

Nawaf Bou-Rabee

We give accurate and ergodic numerical methods for semilinear, second-order Langevin stochastic partial differential equations (SPDE). As a byproduct, we also give good geometric numerical methods for their infinite-dimensional Hamiltonian counterpart. These methods are suitable for Hamiltonian Monte Carlo on Hilbert spaces without preconditioning the underlying Hamiltonian dynamics. A key tool in our approach is Krein's theory on strong stability of symplectic maps, which gives us sufficient conditions for stability of symplectic splitting schemes in highly oscillatory Hamiltonian problems.

Nawaf Bou-Rabee, Eric Vanden-Eijnden

Metropolized integrators for ergodic stochastic differential equations (SDE) are proposed which (i) are ergodic with respect to the (known) equilibrium distribution of the SDE and (ii) approximate pathwise the solutions of the SDE on finite time intervals. Both these properties are demonstrated in the paper and precise strong error estimates are obtained. It is also shown that the Metropolized integrator retains these properties even in situations where the drift in the SDE is nonglobally Lipschitz, and vanilla explicit integrators for SDEs typically become unstable and fail to be ergodic.