A Function Space HMC Algorithm With Second Order Langevin Diffusion Limit

Introduction

Markov chain Monte Carlo (MCMC) algorithms for sampling from high dimensional probability distributions constitute an important part of Bayesian statistical inference. This paper is focussed on the design and analysis of such algorithms to sample a probability distribution on an infinite dimensional Hilbert space ${\mathcal{H}}$ defined via a density with respect to a Gaussian; such problems arise in the Bayesian approach to inverse problems (or Bayesian nonparametrics) [Stua:10] and in the theory of conditioned diffusion processes [Hair:Stua:Voss:10]. Metropolis-Hastings algorithms [Hast:70] constitute a popular class of MCMC methods for sampling an arbitrary probability measure. They proceed by constructing an irreducible, reversible Markov chain by first proposing a candidate move and then accepting it with a certain probability. The acceptance probability is chosen so as to preserve the detailed balance condition ensuring reversibility. In this work, we build on the generalized Hybrid Monte Carlo (HMC) method of [H91] to construct a new non-reversible MCMC method appropriate for sampling measures defined via density with respect to a Gaussian measure on a Hilbert space. We also demonstrate that, for a particular set of parameter values in the algorithm, there is a natural diffusion limit to the second order Langevin (SOL) equation with invariant measure given by the target. We thus name the new method the SOL-HMC algorithm. Our construction is motivated by the following two key design principles:

non-reversible MCMC algorithms, which are hence not from the Metropolis-Hastings class, can have better sampling properties in comparison with their reversible counterparts.

The idea behind the first principle is explained in [David] which surveys a range of algorithms designed specifically to sample measures defined via a density with respect to a Gaussian; the unifying theme is that the proposal is reversible with respect to the underlying Gaussian so that the accept-reject mechanism depends only on the likelihood function and not the prior distribution. The second principle above is also well-documented: non-reversible Markov chains, often constructed by performing individual time-reversible For a definition of time-reversibility see Section 2.3. steps successively [hwang1993, hwang2005], or by building on Hamiltonian mechanics [H91, diac:etal:2000, neal2010mcmc]), may have better mixing properties.

Here we employ similar techniques as that of [PST13] to study our new non-reversible MCMC method, and show that, after appropriate rescaling, it converges to a second order non-reversible Langevin diffusion. Our new algorithm is inspired by similar algorithms in finite dimensions, starting with the work of [H91], who showed how the momentum updates could be correlated in the original HMC method of [Duane1987216], and the more recent work [bou2011patch] which made the explicit connection to second order Langevin diffusions; a helpful overview and discussion may be found in [neal2010mcmc]. Diffusion limit results similar to ours are proved in [bou2011patch, Bou] for finite dimensional problems. In those papers an accept-reject mechanism is appended to various standard integrators for the first and second order Langevin equations, and shown not to destroy the strong pathwise convergence of the underlying methods. The reason for this is that rejections are rare when small time-steps are used. The same reasoning underlies the results we present here, although we consider an infinite dimensional setting and use only weak convergence. Another existing work underpinning that presented here is the paper [BPSS11] which generalizes the hybrid Monte Carlo method for measures defined via density with respect to a Gaussian so that it applies on Hilbert space. Indeed the algorithm we introduce in this paper includes the one from [BPSS11] as a special case and uses the split-step (non-Verlet) integrator first used there. The key idea of the splitting employed is to split according to linear and nonlinear dynamics within the numerical Hamiltonian integration step of the algorithm, rather than according to position and momentum. This allows for an algorithm which exactly preserves the underlying Gaussian reference measure, without rejections, and is key to the fact that the methods are defined on Hilbert space even in the the non-Gaussian case.

for a real valued functional $\Psi$ (which denotes the negative log-likelihood in the case of Bayesian inference) and $M_{\Psi}$ a normalizing constant. Although the above formulation may appear quite abstract, we emphasize that this points to the wide-ranging applicability of our theory: the setting encompasses a large class of models arising in practice, including nonparametric regression using Gaussian random fields and statistical inference for diffusion processes and bridge sampling [Hair:Stua:Voss:10, Stua:10].

In Section 2 we introduce our new algorithm. We start in a finite dimensional context and then explain parametric choices made with reference to the high or infinite dimensional setting. We demonstrate that various other algorithms defined on Hilbert space, such as the function space MALA [Besk:etal:08] and function space HMC algorithms [BPSS11], are special cases. In Section 3 we describe the infinite dimensional setting in full and, in particular, detail the relationship between the change of measure, encapsulated in $\Psi$, and the properties of the Gaussian prior $\pi_{0}.$ Section LABEL:sec:infinite contains the theory of the SPDE which both motivates our class of algorithms, and acts as a limiting process for a specific instance of our algorithm applied on a sequence of spaces of increasing dimension $N$. We prove existence and uniqueness of solutions to the SPDE and characterize its invariant measure. Section LABEL:sec:difflim contains statement of the key diffusion limit Theorem LABEL:thm:difflim. Whilst the structure of the proof is outlined in some detail, various technical estimates are left for Appendices A and B. Section LABEL:sec:num contains some numerics illustrating the new algorithm in the context of a problem from the theory of conditioned diffusions. We make some brief concluding remarks in Section LABEL:sec:conc.

The new algorithm proposed and analyzed in this paper is of interest for two primary reasons. Firstly, it contains a number of existing function space algorithms as special cases and hence plays a useful conceptual role in unifying these methods. Secondly numerical evidence demonstrates that the method is comparable in efficiency to the function space HMC method introduced in [BPSS11] for a test problem arising in conditioned diffusions; until now, the function space HMC method was the clear best choice as demonstrated numerically in [BPSS11]. Furthermore, our numerical results indicate that for certain parameter choices in the SOL-HMC algorithm, and for certain target measures, we are able to improve upon the performance of the function space HMC algorithm, corroborating a similar observation made in [H91] for the finite dimensional samplers that form motivation for the new family of algorithms that we propose here. From a technical point of view the diffusion limit proved in this paper is similar to that proved for the function space MALA in [Pill:Stu:Thi:12], extending to the non-reversible case; however significant technical issues arise which are not present in the reversible case and, in particular, incorporating momentum flips into the analysis, which occur for every rejected step, requires new ideas.

The SOL-HMC Algorithm

where $\mathcal{L}=\mathcal{C}^{-1}$. The corresponding canonical Hamiltonian differential equation is given by

This equation preserves any smooth function of $\mathsf{H}(q,p)$ and, as a consequence, the Liouville equation corresponding to (2.1) preserves the probability density of $\Pi^{\prime}(q,p)$, which is proportional to \exp\big{(}-\mathsf{H}(q,p)\bigr{)}. This fact is the basis for HMC methods [Duane1987216] which randomly sample momentum from the Gaussian $N(0,{\mathcal{M}})$ and then run the Hamiltonian flow for $T$ time units; the resulting Markov chain on $q$ is $\pi(q)$ invariant. In practice the Hamiltonian flow must be integrated numerically, but if a suitable integrator is used (volume-preserving and time-reversible) then a simple accept-reject compensation corrects for numerical error.

and {equs}J = ( 0I-I0 ). Then the Hamiltonian system can be written as

where, abusing notation, $\mathsf{H}(z):=\mathsf{H}(q,p).$ The equation (2.3) preserves the measure $\Pi^{\prime}$.

whilst the equation for $p$ is simply the Ornstein-Uhlenbeck process:

Discretizing the Langevin equation (respectively the random walk found by ignoring the drift) and adding an accept-reject mechanism, leads to the Metropolis-Adjusted Langevin (MALA) (respectively the Random Walk Metropolis (RWM) algorithm).

A natural idea is to try and combine benefits of the HMC algorithm, which couples the position and momentum coordinates, with the MALA and RWM methods. This thought experiment suggests considering the second order Langevin equationPhysicists often refer to this as the Langevin equation for the choice ${\mathcal{K}}_{1}\equiv 0$ which leads to noise only appearing in the momentum equation.

which also preserves $\Pi^{\prime}$ as a straightforward calculation with the Fokker-Planck equation shows.

2. Velocity Rather Than Momentum

Our paper is concerned with using the equation (2.4) to motivate proposals for MCMC. In particular we will be interested in choices of the matrices ${\mathcal{M}}$, $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ which lead to well-behaved algorithms in the limit of large $N$. To this end we write the equation (2.4) in position and momentum coordinates as

In our subsequent analysis, which concerns the large $N$ limit, it turns out to be useful to work with velocity rather than momentum coordinates; this is because the optimal algorithms in this limit are based on ensuring that the velocity and position coordinates all vary on the same scale. For this reason we introduce $v={\mathcal{M}}^{-1}p$ and rewrite the equations as

In the infinite dimensional setting, i.e., when $\cal H$ is an infinite dimensional Hilbert space, this equation is still well posed (see (2.6) below and Theorem LABEL:t:ode). However in this case $W_{1}$ and $W_{2}$ are cylindrical Wiener processes on $\cal H$ (see Section 3.1) and $\mathcal{L}=\mathcal{C}^{-1}$ is necessarily an unbounded operator on ${\cal H}$ because the covariance operator $\cal{C}$ is trace class on $\cal{H}$. The unbounded operators introduce undesirable behaviour in the large $N$ limit when we approximate them; thus we choose $\mathcal{M}$ and the $\mathcal{K}_{i}$ to remove the appearance of unbounded operators. To this end we set ${\mathcal{M}}=\mathcal{L}={\mathcal{C}}^{-1}$, $\mathcal{K}_{1}=\Gamma_{1}{\mathcal{C}}$ and $\mathcal{K}_{2}=\Gamma_{2}{\mathcal{C}}^{-1}$ and assume that $\Gamma_{1}$ and $\Gamma_{2}$ commute with ${\mathcal{C}}$ to obtain the equations

In the above $B_{1}$ and $B_{2}$ are $\cal H$-valued Brownian motions with covariance operator ${\mathcal{C}}$. This equation is well-behaved in infinite dimensions provided that the $\Gamma_{i}$ are bounded operators, and under natural assumptions relating the reference measure, via its covariance ${\mathcal{C}}$, and the log density $\Psi$, which is a real valued functional defined on an appropriate subspace of $\cal H$. Detailed definitions and assumptions regarding (2.6) are contained in the next Section 3. Under such assumptions the function

has desirable properties (see Lemma LABEL:lem:lipschitz+taylor), making the existence theory for (2.6) straightforward. We develop such theory in Section LABEL:sec:infinite – see Theorem LABEL:t:ode. Furthermore, in Theorem LABEL:t:flowprev we will also prove that equation (2.6) preserves the measure $\Pi(dq,dv)$ defined by

where $\Pi_{0}$ is the independent product of $N(0,{\mathcal{C}})$ with itself. The measure $\Pi$ (resp. $\Pi_{0}$) is simply the measure $\Pi^{\prime}$ (resp. $\Pi_{0}^{\prime}$) in the case ${\mathcal{M}}={\mathcal{C}}^{-1}$ and rewritten in $(q,v)$ coordinates instead of $(q,p)$. In finite dimensions the invariance of $\Pi$ follows from the discussions concerning the invariance of $\Pi^{\prime}$.

3. Function Space Algorithm

We note that the choice $\Gamma_{1}\equiv 0$ gives the standard (physicists) Langevin equation

In this section we describe an MCMC method designed to sample the measure $\Pi$ given by (2.8) and hence, by marginalization, the measure $\pi$ given by (1.1). The method is based on discretization of the second order Langevin equation (2.9), written as the hypo-elliptic first order equation (2.10) below. In the finite dimensional setting a method closely related to the one that we introduce was proposed in [H91]; however we will introduce different Hamiltonian solvers which are tuned to the specific structure of our measure, in particular to the fact that it is defined via density with respect to a Gaussian. We will be particularly interested in choices of parameters in the algorithm which ensure that the output (suitability interpolated to continuous time) behaves like (2.9) whilst, as is natural for MCMC methods, exactly preserving the invariant measure. This perspective on discretization of the (physicists) Langevin equation in finite dimensions was introduced in [bou2011patch, Bou].

In position/velocity coordinates, and using (2.7), (2.6) becomes

The algorithm we use is based on splitting (2.10) into an Ornstein-Uhlenbeck (OU) process and a Hamiltonian ODE. The OU process is

To construct volume-preserving and time-reversible integrators the Hamiltonian integration will be performed by a further splitting of (2.12). The usual splitting for the widely used Verlet method is via the velocity and the position coordinates [H91]. Motivated by our infinite dimensional setting, we replace the Verlet integration by the splitting method proposed in [BPSS11]; this leads to an algorithm which is exact (no rejections) in the purely Gaussian case where $\Psi\equiv 0.$ The splitting method proposed in [BPSS11] is via the linear and nonlinear parts of the problem, leading us to consider the two equations

with solution denoted as $(q(t),v(t))=\mathsf{R}^{t}(q(0),v(0))$; and

with solution denoted as $(q(t),v(t))=\Theta^{t}_{1}(q(0),v(0))$. We note that the map {equs}χ^t &= Θ^t/2_1 ∘R^t ∘Θ^t/2_1 is a volume-preserving and time-reversible second order accurate approximation of the Hamiltonian ODE (2.12). We introduce the notation {equs}χ^t_τ &= (χ^t ∘⋯∘χ^t), ⌊τt ⌋ times to denote integration, using this method, up to time ${\tau}$. This integrator can be made to preserve the measure $\Pi$ if appended with a suitable accept-reject rule as detailed below. On the other hand the stochastic map $\Theta^{t}_{0}$ preserves $\Pi$ since it leaves $q$ invariant and since the OU process, which is solved exactly, preserves $\Pi_{0}.$ We now take this idea to define our MCMC method. The infinite dimensional Hilbert space $\mathcal{H}^{s}\times\mathcal{H}^{s}$ in which the chain is constructed will be properly defined in the next section. Here we focus on the algorithm, which will be explained in more details and analyzed in Section LABEL:sec:difflim.

Define the operation ′ so that $v^{\prime}$ is the velocity component of $\Theta^{\delta}_{0}(q,v).$ The preceding considerations suggest that from point $(q^{0},v^{0})\in{\mathcal{H}}^{s}\times{\mathcal{H}}^{s}$ we make the proposal {equs}(q_*^1,v_*^1) = χ^h_τ ∘Θ^δ_0 (q^0,v^0) and that the acceptance probability is given by {equs}α(x^0,ξ^δ): = 1 ∧exp(H(q^0,(v^0)’) - H(q_*^1,v_*^1)), where {equs} H(q,v)=12⟨q, C^-1q⟩+ 12⟨v, C^-1v⟩+ Ψ(q) , $\langle\cdot,\cdot\rangle$ denoting scalar product in $\mathcal{H}$. One step of the resulting MCMC method is then defined by setting

We will make further comments on this algorithm and on the expression (2.3) for the Hamiltonian in Section LABEL:sec:difflim, see Remark LABEL:rem:wellposaccprob. Here it suffices simply to note that whilst $H$ will be almost surely infinite, the energy difference is well-defined for the algorithms we employ. We stress that when the proposal is rejected the chain does not remain in $(q^{0},v^{0})$ but it moves to $(q^{0},-(v^{0})^{\prime})$; that is, the position coordinate stays the same while the velocity coordinate is first evolved according to (2.11) and then the sign is flipped. This flipping of the sign, needed to preserve reversibility, leads to some of the main technical differences with respect to [PST13]; see Remark LABEL:rem:flip. For the finite dimensional case with Verlet integration the form of the accept-reject mechanism and, in particular, the sign-reversal in the velocity, was first derived in [H91] and is discussed in further detail in section 5.3 of [neal2010mcmc]. The algorithm (2.15) preserves $\Pi$ and we refer to it as the SOL-HMC algorithm. Recalling that $v^{\prime}$ denotes the velocity component of $\Theta_{0}^{\delta}(q,v)$, we can equivalently use the notations $\alpha(x,\xi^{\delta})$ and $\alpha(q,{v}^{\prime})$, for $x=(q,v)$ (indeed, by the definition of $\Theta_{0}^{\delta}$, $v^{\prime}$ depends on $\xi^{\delta}$). With this in mind, the pseudo-code for the SOL-HMC is as follows.

Pick $(q^{0},v^{0})\in\mathcal{H}^{s}\times\mathcal{H}^{s}$ and set $k=0$;

given $(q^{k},v^{k})$, define $(v^{k})^{\prime}$ to be the $v$-component of $\Theta^{\delta}_{0}(q^{k},v^{k})$ and calculate the proposal

define the acceptance probability $\alpha(q^{k},(v^{k})^{\prime})$;

set $(q^{k+1},v^{k+1})=(q^{k+1}_{*},v^{k+1}_{*})$ with probability $\alpha(q^{k},(v^{k})^{\prime})$; otherwise set $(q^{k+1},v^{k+1})=(q^{k},-(v^{k})^{\prime})$;

Let Assumption LABEL:ass:1 hold. For any $\delta,h,{\tau}>0$, the Markov chain defined by (2.15) is invariant with respect to $\Pi$ given by (2.8).

We first note that if $\delta\rightarrow\infty$ then the algorithm (2.15) is that introduced in the paper [BPSS11]. From this it follows that, if $\delta=\infty$ and ${\tau}=h$, then the algorithm is simply the funtion-space Langevin introduced in [Besk:etal:08].

Secondly we mention that, in the numerical experiments reported later, we will choose $\Gamma_{2}=I$. The solution of the OU process (2.11) for $v$ is thus given as

where $w\sim N(0,C)$ and $e^{-2\delta}=(1-\iota^{2}).$ The numerical experiments will be described in terms of the parameter $\iota$ rather than $\delta.$

Preliminaries

In this section we detail the notation and the assumptions (Section 3.1 and Section LABEL:sec:assumptions, respectively) that we will use in the rest of the paper.

if $\sum_{j}\lambda_{j}^{2}j^{2r}<\infty$ then $y$ can be equivalently written as

Because ${\mathcal{C}}$ is defined on ${\mathcal{H}}$, the covariance operator