Skip to main content

Research

My current research interests are situated in the overlap between probability and applied mathematics. In addition to the works highlighted on this page, these interests are reflected on my MathOverflow profile page. As highlighted below, my interests span the following topics: Markov Chain Monte Carlo, Markov Chain Mixing Time, Couplings, Stochastic Differential Equations, and Piecewise-deterministic Markov Processes, among other topics. The stochastic processes considered in my research works often display characteristics that frequently present themselves in applications including high-dimensionality, multi-modality, and multi-scale behavior.

For a complete list of my publications, see my Google Scholar profile page.

Tuning Hamiltonian Monte Carlo Samplers

The No-U-Turn Sampler employs a remarkably clever way to locally adapt the path length tuning parameter in HMC to eponymously avoid U-turns in the underlying leapfrog integration legs, and at the same time, preserve detailed balance. Thanks to this self-tuning feature, the No-U-Turn Sampler has become the de facto sampler in many probabilistic programming languages like STAN, PyMC3, NIMBLE, Turing, and NumPyro. The invention of the No-U-Turn Sampler inevitably raises the question: is it possible to generalize this local adaptation strategy to the algorithm’s other tuning parameters? We address this question by developing a simple and flexible framework for self-tuning HMC termed GIST, signifying Gibbs self-tuning for HMC. Remarkably, this framework generalizes/unifies almost all locally adaptive HMC samplers in current use including the No-U-Turn Sampler, the Apogee-to-Apogee Path Sampler, and randomized Hamiltonian Monte Carlo. Besides this, the GIST sampling framework opens the door to: (i) simpler alternatives to the No-U-Turn Sampler for path length tuning; but also (ii) local adaptation of the algorithm’s other parameters, which will be the subject of forthcoming work.

Hamiltonian Monte Carlo Mixing Time

The Hamiltonian Monte Carlo (HMC) sampler is based on a combination of Hamiltonian dynamics and momentum randomizations. The Hamiltonian dynamics is discretized, and the discretization bias can either be borne or corrected by a Metropolis-Hastings accept-reject step. Despite its empirical success, until a few years ago there have been almost no convergence bounds for the algorithm. This situation has changed in the last years where approaches to quantify convergence to equilibrium based on coupling, conductance and hypocoercivity have been developed. In a series of works (see below), we developed a variety of novel coupling approaches to precisely prove the mixing and complexity properties of uadjusted/adjusted HMC in representative models that display high-dimensionality and non-convexity. These couplings offer significant flexibility for extensions including as a potential diagnostic to evaluate HMC convergence.

Randomized Time Integrators for Hamiltonian MCMC

Measure-preserving ODEs, SDEs, and PDMPs play an important role in constructing non-reversible MCMC kernels, and among time integrators for the corresponding flows, randomized time integrators have stood out in terms of their complexity properties. The essential idea behind these randomized time integrators is to replace deterministic quadrature rules used per time integration step with a suitable Monte Carlo integration rule. The resulting “stratified Monte Carlo” time integrator not only requires less regularity on the underlying ODE, SDE or PDMP, but remarkably, also leads to significantly better asymptotic bias properties for the corresponding Markov kernel, and hence, complexity. Variants of these schemes that are Metropolis-adjustable have been developed in Bou-Rabee & Marsden (2022). Moreover, under minimal regularity of the underlying drift coefficients, 2.5/3.5-order randomized Runge-Kutta-Nyström (rRKN) schemes were developed in Bou-Rabee & Kleppe (2023).

Piecewise Deterministic Markov Processes

Andersen dynamics describes a piecewise deterministic Markov process that follows Hamiltonian dynamics and randomizes the velocity of selected components after exponential waiting times. It is closely related to the kinetic Langevin diffusion except that it has neither diffusion nor explicit dissipation and involves a non-local jump operator. It includes as a special case Randomized Hamiltonian Monte Carlo. Building on the geometric ergodicity verification in Bou-Rabee & Sanz-Serna (2017), we used a maximal coupling of the underlying velocity randomizations to obtain an optimal dimension dependence for Andersen dynamics with mean field interactions in Bou-Rabee & Eberle (2022).

Path Integral Molecular Dynamics

Path Integral Molecular Dynamics (PIMD) is based on Feynmann’s path-integral formulation of quantum statistical mechanics, and provides a practical tool to compute static and dynamic quantum statistics, especially when quantum nuclear effects are relevant. Applications include calculations of chemical reaction rates, diffusion coefficients, and absorption spectra. The appeal of PIMD is that it uses a classical ring-polymer-bead system to solve for properties of the quantum system. At the heart of every numerical integration for PIMD is a free ring-polymer update. However, due to fast harmonic motions present in the free ring-polymer system, a strongly stable approximation to the free ring-polymer update is essential. Another requirement is that the integration scheme remains accurate in the infinite bead limit. However, an unfortunate property of existing numerical integration schemes for PIMD — including essentially all existing ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (T-RPMD) methods — is that for a given MD timestep, the overlap between the exact ring-polymer Boltzmann distribution and that sampled using the approximation becomes zero in the continuum or infinite bead limit. A key outcome of this collaboration is that these and other problems can be avoided through the introduction of strongly stable and “dimension-free” numerical integration schemes for which the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite bead limit in relevant models; these conclusions are borne out numerically for quantum liquid water. These properties of strong stability and dimensionality freedom allow for remarkably accurate and stable ergodic time-averaging, and involve a very simple “Cayley” modification of the standard “BAOAB” splitting scheme from molecular dynamics, and therefore, “BCOCB” is expected to be widely used for future PIMD simulations.

Numerical Solution of Stochastic Differential Equations

Sticky diffusion processes are ubiquitous in applications arising in a variety of models of physical and natural processes which are naturally described by a set of variables that can change dimension. However, until a few years ago, there were no methods to simulate them. This situation has changed. With Eric Vanden-Eijnden (NYU), we developed a discretization of SDEs in \emph{space}, rather than time, leading to a Markov jump process (MJP) approximation. Since the increments of the MJP can be uniformly bounded, these algorithms are numerically stable by construction. Another nice feature of the MJP is that it naturally handles boundary conditions on the solution (e.g. positivity constraints in population or interest rate models), which is not at all obvious how to do using standard methods. With Miranda Holmes-Cerfon (UBC), we showed how one can naturally incorporate “sticky” boundary conditions into the MJP and demonstrated that MJP can accelerate simulations considerably, in some cases by orders of magnitude. Beyond sticky diffusions, this MJP allows one to approximate SDEs and SPDEs with trickier characteristics (such as interface conditions) or even in situations where the operators of the underlying processes are not yet well understood.