Using Perturbed Underdamped Langevin Dynamics to Efficiently Sample from Probability Distributions

TL;DR

This paper introduces perturbed underdamped Langevin dynamics for sampling, which maintains the original invariant measure while reducing asymptotic variance, leading to more efficient sampling methods.

Contribution

The paper proposes a novel class of Langevin samplers with perturbations that preserve the invariant measure and improve sampling efficiency, supported by theoretical and numerical analysis.

Findings

Perturbed Langevin samplers maintain the target distribution.

Perturbations can reduce asymptotic variance.

Numerical experiments confirm improved performance.

Abstract

In this paper we introduce and analyse Langevin samplers that consist of perturbations of the standard underdamped Langevin dynamics. The perturbed dynamics is such that its invariant measure is the same as that of the unperturbed dynamics. We show that appropriate choices of the perturbations can lead to samplers that have improved properties, at least in terms of reducing the asymptotic variance. We present a detailed analysis of the new Langevin sampler for Gaussian target distributions. Our theoretical results are supported by numerical experiments with non-Gaussian target measures.

Figures

Equations

$$\pi (f):={\int }_{{\mathbb{R}}^d}f(x)\pi (\mathrm{d}x).$$

$${\pi }_T(f):=\frac{1}{T}{\int }_0^{T}f(X_s)\mathrm{d}s$$

$$\mathrm{d}{X}_t=-\nabla V(X_t)\mathrm{d}t+\sqrt{2}\mathrm{d}{W}_t,$$

$$\mathrm{d}{q}_t$$

$$\mathrm{d}{p}_t$$

$$\hat{\pi }(q,p)=\frac{1}{\hat{Z}}\exp \left(-V(q)-\frac{1}{2}p\cdot M^{-1}p\right),$$

$$\frac{1}{\sqrt{T}}\left({\pi }_T(f)-\pi (f)\right)\underset{T\to \infty }{\overset{d}{\to }}\mathcal{N}(0,2{\sigma }_f^2),$$

$$-\mathcal{L}\varphi =f-\pi (f),{\sigma }_f^2={\int }_{{\mathbb{R}}^d}\varphi (-\mathcal{L}\varphi )\pi (\mathrm{d}x).$$

$$d{X}_t=-P\nabla V(X_t)dt+\sqrt{2P}d{W}_t.$$

$$d{X}_t=-\left(\nabla V(X_t)-\gamma (X_t)\right)dt+\sqrt{2}d{W}_t,$$

$$\mathrm{d}{q}_t$$

$$\mathrm{d}{p}_t$$

$$\mathrm{d}{X}_t=a(X_t)\mathrm{d}t+\sqrt{2}b(X_t)\mathrm{d}{W}_t,$$

$$Lf=a\cdot \nabla f+\Sigma :D^{2}f,f\in C_c^2(R^d)$$

$$P_{t}f(x)={\int }_{R^d}p(t,x,y)f(y)\mathrm{d}y,t\ge 0.$$

$${\int }_{{\mathbb{R}}^d}{P}_{t}f(x)\pi (x)\mathrm{d}x={\int }_{{\mathbb{R}}^d}f(x)\pi (x)\mathrm{d}x,t\ge 0,$$

$$H^1(\pi )=\{f\in L^2(\pi ):\Vert \nabla f{\Vert }_{L^2(\pi )}^2<\infty \}$$

$$\Vert f{\Vert }_{H^1(\pi )}^2=\Vert f{\Vert }_{L^2(\pi )}^2+\Vert \nabla f{\Vert }_{L^2(\pi ).}^2$$

$$a=\Sigma \nabla \log \pi +\nabla \cdot \Sigma +\gamma ,$$

$$\nabla \cdot \left(\pi \gamma \right)=0.$$

$$\gamma =\frac{1}{\pi }\nabla \cdot \left(\pi C\right).$$

$$\mathcal{L}f=\frac{1}{\pi }\nabla \cdot \left((\Sigma +C)\pi \nabla f\right),f\in C_c^2({\mathbb{R}}^d).$$

$$\mathrm{d}{Z}_t=b(Z_t)\mathrm{d}t+\sqrt{2}\sigma (Z_t)\mathrm{d}{W}_t,$$

$$\mathcal{L}=b(z)\cdot {\nabla }_z+\Sigma (z):D_z^2,$$

$${\int }_{R^m}{\pi }_z(x,y)\mathrm{d}y=\pi (x).$$

$$\mathrm{d}{X}_t=-(I+J)\nabla V(X_t)\mathrm{d}t+\sqrt{2}\mathrm{d}{W}_t,$$

$$\gamma((q,p))=\left(\begin{array}[]{c}M^{-1}p\\ -\nabla V(q)\end{array}\right)$$

$$b=\left(\begin{array}[]{c}\boldsymbol{0}\\ \sqrt{\Gamma}\end{array}\right)\in\mathbb{R}^{2d\times d},$$

$$\gamma((q,p))=\left(\begin{array}[]{c}M^{-1}p-\mu J_{1}\nabla V(q)\\ -\nabla V(q)-\nu J_{2}M^{-1}p\end{array}\right)$$

$$b=\left(\begin{array}[]{c}\boldsymbol{0}\\ \sqrt{\Gamma}\end{array}\right)\in\mathbb{R}^{2d\times d},$$

$$\hat{\hat{\pi }}(q,p,u_1,\dots ,u_m)\propto e^{-\frac{|p{|}^2}{2}-\frac{u_1^2+\dots +u_m^2}{2}-V(q)},$$

Full text

1 Introduction and Motivation

Sampling from probability measures in high-dimensional spaces is a problem that appears frequently in applications, e.g. in computational statistical mechanics and in Bayesian statistics. In particular, we are faced with the problem of computing expectations with respect to a probability measure $\pi$ on $\mathbb{R}^{d}$, i.e. we wish to evaluate integrals of the form:

[TABLE]

As is typical in many applications, particularly in molecular dynamics and Bayesian inference, the density (for convenience denoted by the same symbol $\pi$) is known only up to a normalization constant; furthermore, the dimension of the underlying space is quite often large enough to render deterministic quadrature schemes computationally infeasible.

A standard approach to approximating such integrals is Markov Chain Monte Carlo (MCMC) techniques GCS*+* (14); Liu (08); RC (13), where a Markov process $(X_{t})_{t\geq 0}$ is constructed which is ergodic with respect to the probability measure $\pi$. Then, defining the long-time average

[TABLE]

for $f\in L^{1}(\pi)$, the ergodic theorem guarantees almost sure convergence of the long-time average $\pi_{T}(f)$ to $\pi(f)$.

There are infinitely many Markov, and, for the purposes of this paper diffusion, processes that can be constructed in such a way that they are ergodic with respect to the target distribution. A natural question is then how to choose the ergodic diffusion process $(X_{t})_{t\geq 0}$. Naturally the choice should be dictated by the requirement that the computational cost of (approximately) calculating (1) is minimized. A standard example is given by the overdamped Langevin dynamics defined to be the unique (strong) solution $(X_{t})_{t\geq 0}$ of the following stochastic differential equation (SDE):

[TABLE]

where $V=-\log\pi$ is the potential associated with the smooth positive density $\pi$. Under appropriate assumptions on $V$, i.e. on the measure $\pi(\mathrm{d}x)$, the process $(X_{t})_{t\geq 0}$ is ergodic and in fact reversible with respect to the target distribution.

Another well-known example is the underdamped Langevin dynamics given by $(X_{t})_{t\geq 0}=(q_{t},p_{t})_{t\geq 0}$ defined on the extended space (phase space) $\mathbb{R}^{d}\times\mathbb{\mathbb{R}}^{d}$ by the following pair of coupled SDEs:

[TABLE]

where mass and friction tensors $M$ and $\Gamma$, respectively, are assumed to be symmetric positive definite matrices. It is well-known Pav (14); LS (16) that $(q_{t},p_{t})_{t\geq 0}$ is ergodic with respect to the measure $\widehat{\pi}:=\pi\otimes\mathcal{N}(0,M)$, having density with respect to the Lebesgue measure on $\mathbb{R}^{2d}$ given by

[TABLE]

where $\widehat{Z}$ is a normalization constant. Note that $\widehat{\pi}$ has marginal $\pi$ with respect to $p$ and thus for functions $f\in L^{1}(\pi)$, we have that $\frac{1}{t}\int_{0}^{t}f(q_{t})\,\mathrm{d}t\rightarrow\pi(f)$ almost surely. Notice also that the dynamics restricted to the $q$-variables is no longer Markovian. The $p$-variables can thus be interpreted as giving some instantaneous memory to the system, facilitating efficient exploration of the state space. Higher order Markovian models, based on a finite dimensional (Markovian) approximation of the generalized Langevin equation can also be used CBP (09).

As there is a lot of freedom in choosing the dynamics in (2), see the discussion in Section 2, it is desirable to choose the diffusion process $(X_{t})_{t\geq 0}$ in such a way that $\pi_{T}(f)$ can provide a good estimation of $\pi(f)$. The performance of the estimator (2) can be quantified in various manners. The ultimate goal, of course, is to choose the dynamics as well as the numerical discretization in such a way that the computational cost of the longtime-average estimator is minimized, for a given tolerance. The minimization of the computational cost consists of three steps: bias correction, variance reduction and choice of an appropriate discretization scheme. For the latter step see Section 5 and (DLP, 16, Sec. 6).

Under appropriate conditions on the potential $V$ it can be shown that both (3) and (4) converge to equilibrium exponentially fast, e.g. in relative entropy. One performance objective would then be to choose the process $(X_{t})_{t\geq 0}$ so that this rate of convergence is maximised. Conditions on the potential $V$ which guarantee exponential convergence to equilibrium, both in $L^{2}(\pi)$ and in relative entropy can be found in MV (00); BGL (13). A powerful technique for proving exponentially fast convergence to equilibrium that will be used in this paper is C. Villani’s theory of hypocoercivity Vil (09). In the case when the target measure $\pi$ is Gaussian, both the overdamped  (3) and the underdamped (4) dynamics become generalized Ornstein-Uhlenbeck processes. For such processes the entire spectrum of the generator – or, equivalently, the Fokker-Planck operator – can be computed analytically and, in particular, an explicit formula for the $L^{2}$-spectral gap can be obtained MPP (02); OPPS (12, 15). A detailed analysis of the convergence to equilibrium in relative entropy for stochastic differential equations with linear drift, i.e. generalized Ornstein-Uhlenbeck processes, has been carried out in AE (14).

In addition to speeding up convergence to equilibrium, i.e. reducing the bias of the estimator (2), one is naturally also interested in reducing the asymptotic variance. Under appropriate conditions on the target measure $\pi$ and the observable $f$, the estimator $\pi_{T}(f)$ satisfies a central limit theorem (CLT) KLO (12), that is,

[TABLE]

where $\sigma^{2}_{f}<\infty$ is the asymptotic variance of the estimator $\pi_{T}(f)$. The asymptotic variance characterises how quickly fluctuations of $\pi_{T}(f)$ around $\pi(f)$ contract to [math]. Consequently, another natural objective is to choose the process $(X_{t})_{t\geq 0}$ such that $\sigma^{2}_{f}$ is as small as possible. It is well known that the asymptotic variance can be expressed in terms of the solution to an appropriate Poisson equation for the generator of the dynamics KLO (12)

[TABLE]

Techniques from the theory of partial differential equations can then be used in order to study the problem of minimizing the asymptotic variance. This is the approach that was taken in DLP (16), see also HNW (15), and it will also be used in this paper.

Other measures of performance have also been considered. For example, in RBS15b ; RBS15a , performance of the estimator is quantified in terms of the rate functional of the ensemble measure $\frac{1}{t}\int_{0}^{t}\delta_{X(t)}(dx)$. See also JO (10) for a study of the nonasymptotic behaviour of MCMC techniques, including the case of overdamped Langevin dynamics.

Similar analyses have been carried out for various modifications of (3). Of particular interest to us are the Riemannian manifold MCMC GC (11) (see the discussion in Section 2) and the nonreversible Langevin samplers HHMS (93, 05). As a particular example of the general framework that was introduced in GC (11), we mention the preconditioned overdamped Langevin dynamics that was presented in AMO (16)

[TABLE]

In this paper, the long-time behaviour of as well as the asymptotic variance of the corresponding estimator $\pi_{T}(f)$ are studied and applied to equilibrium sampling in molecular dynamics. A variant of the standard underdamped Langevin dynamics that can be thought of as a form of preconditioning and that has been used by practitioners is the mass-tensor molecular dynamics Ben (75).

The nonreversible overdamped Langevin dynamics

[TABLE]

where the vector field $\gamma$ satisfies $\nabla\cdot(\pi\gamma)=0$ is ergodic (but not reversible) with respect to the target measure $\pi$ for all choices of the divergence-free vector field $\gamma$. The asymptotic behaviour of this process was considered for Gaussian diffusions in HHMS (93), where the rate of convergence of the covariance to equilibrium was quantified in terms of the choice of $\gamma$. This work was extended to the case of non-Gaussian target densities, and consequently for nonlinear SDEs of the form (8) in HHMS (05). The problem of constructing the optimal nonreversible perturbation, in terms of the $L^{2}(\pi)$ spectral gap for Gaussian target densities was studied in  LNP (13) see also WHC (14). Optimal nonreversible perturbations with respect to miniziming the asymptotic variance were studied in DLP (16); HNW (15). In all these works it was shown that, in theory (i.e. without taking into account the computational cost of the discretization of the dynamics (8)), the nonreversible Langevin sampler (8) always outperforms the reversible one (3), both in terms of converging faster to the target distribution as well as in terms of having a lower asymptotic variance. It should be emphasized that the two optimality criteria, maximizing the spectral gap and minimizing the asymptotic variance, lead to different choices for the nonreversible drift $\gamma(x)$.

The goal of this paper is to extend the analysis presented in DLP (16); LNP (13) by introducing the following modification of the standard underdamped Langevin dynamics:

[TABLE]

where $M,\Gamma\in\mathbb{R}^{d\times d}$ are constant strictly positive definite matrices, $\mu$ and $\nu$ are scalar constants and $J_{1},J_{2}\in\mathbb{R}^{d\times d}$ are constant skew-symmetric matrices. As demonstrated in Section 3, the process defined by (9) will be ergodic with respect to the Gibbs measure $\widehat{\pi}$ defined in (5).

Our objective is to investigate the use of these dynamics for computing ergodic averages of the form (2). To this end, we study the long time behaviour of (9) and, using hypocoercivity techniques, prove that the process converges exponentially fast to equilibrium. This perturbed underdamped Langevin process introduces a number of parameters in addition to the mass and friction tensors which must be tuned to ensure that the process is an efficient sampler. For Gaussian target densities, we derive estimates for the spectral gap and the asymptotic variance, valid in certain parameter regimes. Moreover, for certain classes of observables, we are able to identify the choices of parameters which lead to the optimal performance in terms of asymptotic variance. While these results are valid for Gaussian target densities, we advocate these particular parameter choices also for more complex target densities. To demonstrate their efficacy, we perform a number of numerical experiments on more complex, multimodal distributions. In particular, we use the Langevin sampler (9) in order to study the problem of diffusion bridge sampling.

The rest of the paper is organized as follows. In Section 2 we present some background material on Langevin dynamics, we construct general classes of Langevin samplers and we introduce criteria for assessing the performance of the samplers. In Section 3 we study qualitative properties of the perturbed underdamped Langevin dynamics (9) including exponentially fast convergence to equilibrium and the overdamped limit. In Section 4 we study in detail the performance of the Langevin sampler (9) for the case of Gaussian target distributions. In Section 5 we introduce a numerical scheme for simulating the perturbed dynamics (9) and we present numerical experiments on the implementation of the proposed samplers for the problem of diffusion bridge sampling. Section 6 is reserved for conclusions and suggestions for further work. Finally, the appendices contain the proofs of the main results presented in this paper and of several technical results.

2 Construction of General Langevin Samplers

2.1 Background and Preliminaries

In this section we consider estimators of the form (2) where $(X_{t})_{t\geq 0}$ is a diffusion process given by the solution of the following Itô SDE:

[TABLE]

with drift coefficient $a:\operatorname{\mathbb{R}}^{d}\rightarrow\operatorname{\mathbb{R}}^{d}$ and diffusion coefficient $b:\operatorname{\mathbb{R}}^{d}\rightarrow\operatorname{\mathbb{R}}^{d\times m}$ both having smooth components, and where $(W_{t})_{t\geq 0}$ is a standard $\operatorname{\mathbb{R}}^{m}$–valued Brownian motion. Associated with (10) is the infinitesimal generator $\operatorname{\mathcal{L}}$ formally given by

[TABLE]

where $\Sigma=bb^{\top}$, $D^{2}f$ denotes the Hessian of the function $f$ and $\,:\,$ denotes the Frobenius inner product. In general, $\Sigma$ is nonnegative definite, and could possibly be degenerate. In particular, the infinitesimal generator (11) need not be uniformly elliptic. To ensure that the corresponding semigroup exhibits sufficient smoothing behaviour, we shall require that the process (10) is hypoelliptic in the sense of Hörmander. If this condition holds, then irreducibility of the process $(X_{t})_{t\geq 0}$ will be an immediate consequence of the existence of a strictly positive invariant distribution $\pi(x)\mathrm{d}x$, see Kli (87).

Suppose that $(X_{t})_{t\geq 0}$ is nonexplosive. It follows from the hypoellipticity assumption that the process $(X_{t})_{t\geq 0}$ possesses a smooth transition density $p(t,x,y)$ which is defined for all $t\geq 0$ and $x,y\in\operatorname{\mathbb{R}}^{d}$, (Bas, 98, Theorem VII.5.6). The associated strongly continuous Markov semigroup $(P_{t})_{t\geq 0}$ is defined by

[TABLE]

Suppose that $(P_{t})_{t\geq 0}$ is invariant with respect to the target distribution $\pi(x)\,\mathrm{d}x$, i.e.

[TABLE]

for all bounded continuous functions $f$. Then $(P_{t})_{t\geq 0}$ can be extended to a positivity preserving contraction semigroup on $L^{2}(\pi)$ which is strongly continuous. Moreover, the infinitesimal generator corresponding to $(P_{t})_{t\geq 0}$ is given by an extension of $(\operatorname{\mathcal{L}},C^{2}_{c}(\operatorname{\mathbb{R}}^{d}))$, also denoted by $\operatorname{\mathcal{L}}$.

Due to hypoellipticity, the probability measure $\pi$ on $\mathbb{R}^{d}$ has a smooth and positive density with respect to the Lebesgue measure, and (slightly abusing the notation) we will denote this density also by $\pi$. Let $L^{2}(\pi)$ be the Hilbert space of $\pi$-square integrable functions equipped with inner product $\left\langle\cdot,\cdot\right\rangle_{L^{2}(\pi)}$ and norm $\left\lVert{\cdot}\right\rVert_{L^{2}(\pi)}$. We will also make use of the Sobolev space

[TABLE]

of $L^{2}(\pi)$-functions with weak derivatives in $L^{2}(\pi)$, equipped with norm

[TABLE]

2.2 A General Characterisation of Ergodic Diffusions

A natural question is what conditions on the coefficients $a$ and $b$ of (10) are required to ensure that $(X_{t})_{t\geq 0}$ is invariant with respect to the distribution $\pi(x)\,\mathrm{d}x$. The following result provides a necessary and sufficient condition for a diffusion process to be invariant with respect to a given target distribution.

Theorem 2.1.

Consider a diffusion process $(X_{t})_{t\geq 0}$ on $\mathbb{R}^{d}$ defined by the unique, non-explosive solution to the Itô SDE (10) with drift $a\in C^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$ and diffusion coefficient $b\in C^{1}(\mathbb{R}^{d};\mathbb{R}^{d\times m})$. Then $(X_{t})_{t\geq 0}$ is invariant with respect to ${\pi}$ if and only if

[TABLE]

where $\Sigma=bb^{\top}$ and $\gamma:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ is a continuously differentiable vector field satisfying

[TABLE]

If additionally $\gamma\in L^{1}({\pi})$, then there exists a skew-symmetric matrix function $C:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ such that

[TABLE]

In this case the infinitesimal generator can be written as an $L^{2}(\pi)$-extension of

[TABLE]

The proof of this result can be found in (Pav, 14, Ch. 4); similar versions of this characterisation can be found in Vil (09) and HHMS (05). See also MCF (15).

Remark 1.

If (14) holds and $\mathcal{L}$ is hypoelliptic it follows immediately that $(X_{t})_{t\geq 0}$ is ergodic with unique invariant distribution $\pi(x)\,\mathrm{d}x$.

More generally, we can consider Itô diffusions in an extended phase space:

[TABLE]

where $(W_{t})_{t\geq 0}$ is a standard Brownian motion in $\operatorname{\mathbb{R}}^{N}$, $N\geq d$. This is a Markov process with generator

[TABLE]

where \Sigma(z)=\big{(}\sigma\sigma^{T}\big{)}(z). We will consider dynamics $(Z_{t})_{t\geq 0}$ that is ergodic with respect to $\pi_{z}(z)\,\mathrm{d}z$ such that

[TABLE]

where $z=(x,\,y),\;x\in\operatorname{\mathbb{R}}^{d},\,y\in\operatorname{\mathbb{R}}^{m},\;d+m=N$.

There are various well-known choices of dynamics which are invariant (and indeed ergodic) with respect to the target distribution $\pi(x)\mathrm{d}x$.

Choosing $b=I$ and $\gamma=0$ we immediately recover the overdamped Langevin dynamics (3). 2. 2.

Choosing $b=I$, and $\gamma\neq 0$ such that (15) holds gives rise to the nonreversible overdamped equation defined by (8). As it satisfies the conditions of Theorem 2.1, it is ergodic with respect to $\pi$. In particular choosing $\gamma(x)=J\nabla V(x)$ for a constant skew-symmetric matrix $J$ we obtain

[TABLE]

which has been studied in previous works. 3. 3.

Given a target density $\pi>0$ on $\operatorname{\mathbb{R}}^{d}$, if we consider the augmented target density $\widehat{\pi}$ on $\operatorname{\mathbb{R}}^{2d}$ given in (5), then choosing

[TABLE]

and

[TABLE]

where $M$ and $\Gamma$ are positive definite symmetric matrices, the conditions of Theorem 2.1 are satisfied for the target density $\widehat{\pi}$. The resulting dynamics $(q_{t},p_{t})_{t\geq 0}$ is determined by the underdamped Langevin equation (4). It is straightforward to verify that the generator is hypoelliptic, (LRS, 10, Sec 2.2.3.1), and thus $(q_{t},p_{t})_{t\geq 0}$ is ergodic. 4. 4.

More generally, consider the augmented target density $\widehat{\pi}$ on $\mathbb{R}^{2d}$ as above, and choose

[TABLE]

and

[TABLE]

where $\mu$ and $\nu$ are scalar constants and $J_{1},J_{2}\in\mathbb{R}^{d\times d}$ are constant skew-symmetric matrices. With this choice we recover the perturbed Langevin dynamics (9). It is straightforward to check that (22) satisfies the invariance condition (15), and thus Theorem 2.1 guarantees that (9) is invariant with respect to $\widehat{\pi}$. 5. 5.

In a similar fashion, one can introduce an augmented target density on $\mathbb{R}^{(m+2)d}$, with

[TABLE]

where $p,q,u_{i}\in\mathbb{R}^{d}$, for $i=1,\ldots,m$. Clearly $\int_{\mathbb{R}^{d}\times\mathbb{R}^{md}}\widehat{\widehat{\pi}}(q,p,u_{1},\ldots,u_{m})\,\mathrm{d}p\,\mathrm{d}u_{1}\,\ldots\mathrm{d}u_{m}=\pi(q)$. We now define $\gamma:\mathbb{R}^{(m+2)d}\rightarrow\mathbb{R}^{(m+2)d}$ by

[TABLE]

and $b:\mathbb{R}^{(m+2)d}\rightarrow\mathbb{R}^{(m+2)d\times(m+2)d}$ by

[TABLE]

where $\lambda_{i}\in\mathbb{R}$ and $\alpha_{i}>0$, for $i=1,\ldots,m$. The resulting process (10) is given by

[TABLE]

where $(W^{1}_{t})_{t\geq 0},\ldots(W^{m}_{t})_{t\geq 0}$ are independent $\mathbb{R}^{d}$–valued Brownian motions. This process is ergodic with unique invariant distribution $\widehat{\widehat{\pi}}$, and under appropriate conditions on $V$, converges exponentially fast to equilibrium in relative entropy OP (11). Equation (24) is a Markovian representation of a generalised Langevin equation of the form

[TABLE]

where $N(t)$ is a mean-zero stationary Gaussian process with autocorrelation function $F(t)$, i.e.

[TABLE]

and

[TABLE] 6. 6.

Let $\widetilde{\pi}(z)\propto\exp(-\Phi(z))$ be a positive density on $\mathbb{R}^{N}$ where $N>d$ such that

[TABLE]

where $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{N-d}$. Then choosing $b=I_{D\times D}$ and $\gamma=0$ we obtain the dynamics

[TABLE]

then $(X_{t},Y_{t})_{t\geq 0}$ is immediately ergodic with respect to $\widetilde{\pi}$.

2.3 Comparison Criteria

For a fixed observable $f$, a natural measure of accuracy of the estimator $\pi_{T}(f)=t^{-1}\int_{0}^{t}f(X_{s})\,\mathrm{d}s$ is the mean square error (MSE) defined by

[TABLE]

where $\mathbb{E}_{x}$ denotes the expectation conditioned on the process $(X_{t})_{t\geq 0}$ starting at $x$. It is instructive to introduce the decomposition $MSE(f,T)=\mu^{2}(f,T)+\sigma^{2}(f,T)$, where

[TABLE]

Here $\mu(f,T)$ measures the bias of the estimator $\pi_{T}(f)$ and $\sigma^{2}(f,T)$ measures the variance of fluctuations of $\pi_{T}(f)$ around the mean.

The speed of convergence to equilibrium of the process $(X_{t})_{t\geq 0}$ will control both the bias term $\mu(f,T)$ and the variance $\sigma^{2}(f,T)$. To make this claim more precise, suppose that the semigroup $(P_{t})_{t\geq 0}$ associated with $(X_{t})_{t\geq 0}$ decays exponentially fast in $L^{2}(\pi)$, i.e. there exist constants $\lambda>0$ and $C\geq 1$ such that

[TABLE]

Remark 2.

If (27) holds with $C=1$, this estimate is equivalent to $-\operatorname{\mathcal{L}}$ having a spectral gap in $L^{2}(\pi)$. Allowing for a constant $C>1$ is essential for our purposes though in order to treat nonreversible and degenerate diffusion processes by the theory of hypocoercivity as outlined in Vil (09).

The following lemma characterises the decay of the bias $\mu(f,T)$ as $T\rightarrow\infty$ in terms of $\lambda$ and $C$. The proof can be found in Appendix A.

Lemma 1.

Let $(X_{t})_{t\geq 0}$ be the unique, non-explosive solution of (10), such that $X_{0}\sim\pi_{0}\ll\pi$ and $\frac{d\pi_{0}}{d\pi}\in L^{2}(\pi)$, where $\frac{d\pi_{0}}{d\pi}$ denotes the Radon-Nikodym derivative of $\pi_{0}$ with respect to $\pi$. Suppose that the process is ergodic with respect to $\pi$ such that the Markov semigroup $(P_{t})_{t\geq 0}$ satisfies (27). Then for $f\in L^{\infty}(\pi)$,

[TABLE]

The study of the behaviour of the variance $\sigma^{2}(f,T)$ involves deriving a central limit theorem for the additive functional $\int_{0}^{t}f(X_{t})-\pi(f)\,\mathrm{d}t$. As discussed in CCG (12), we reduce this problem to proving well-posedness of the Poisson equation

[TABLE]

The only complications in this approach arise from the fact that the generator $\operatorname{\mathcal{L}}$ need not be symmetric in $L^{2}(\pi)$ nor uniformly elliptic. The following result summarises conditions for the well-posedness of the Poisson equation and it also provides with us with a formula for the asymptotic variance. The proof can be found in Appendix A.

Lemma 2.

Let $(X_{t})_{t\geq 0}$ be the unique, non-explosive solution of (10) with smooth drift and diffusion coefficients, such that the corresponding infinitesimal generator is hypoelliptic. Syppose that $(X_{t})_{t\geq 0}$ is ergodic with respect to $\pi$ and moreover, $(P_{t})_{t\geq 0}$ decays exponentially fast in $L^{2}(\pi)$ as in (27). Then for all $f\in L^{2}(\pi)$, there exists a unique mean zero solution $\chi$ to the Poisson equation (28). If $X_{0}\sim\pi$, then for all $f\in C^{\infty}(\mathbb{R}^{d})\cap L^{2}(\pi)$

[TABLE]

where $\sigma^{2}_{f}$ is the asymptotic variance defined by

[TABLE]

Moreover, if $X_{0}\sim\pi_{0}$ where $\pi_{0}\ll\pi$ and $\frac{d\pi_{0}}{d\pi}\in L^{2}(\pi)$ then (29) holds for all $f\in C^{\infty}(\mathbb{R}^{d})\cap L^{\infty}(\pi)$.

Clearly, observables that only differ by a constant have the same asymptotic variance. In the sequel, we will hence restrict our attention to observables $f\in L^{2}(\pi)$ satisfying $\pi(f)=0$, simplifying expressions (28) and (29). The corresponding subspace of $L^{2}(\pi)$ will be denoted by

[TABLE]

If the exponential decay estimate (27) is satisfied, then Lemma 2 shows that $-\operatorname{\mathcal{L}}$ is invertible on $L^{2}_{0}(\pi)$, so we can express the asymptoptic variance as

[TABLE]

Let us also remark that from the proof of Lemma 2 it follows that the inverse of $\mathcal{L}$ is given by

[TABLE]

We note that the constants $C$ and $\lambda$ appearing in the exponential decay estimate (27) also control the speed of convergence of $\sigma^{2}(f,T)$ to zero. Indeed, it is straightforward to show that if (27) is satisfied, then the solution $\chi$ of (28) satisfies

[TABLE]

Lemmas 1 and 2 would suggest that choosing the coefficients $\Sigma$ and $\gamma$ to optimize the constants $C$ and $\lambda$ in (34) would be an effective means of improving the performance of the estimator $\pi_{T}(f)$, especially since the improvement in performance would be uniform over an entire class of observables. When this is possible, this is indeed the case. However, as has been observed in LNP (13); HHMS (93, 05), maximising the speed of convergence to equilibrium is a delicate task. As the leading order term in $MSE(f,T)$, it is typically sufficient to focus specifically on the asymptotic variance $\sigma^{2}_{f}$ and study how the parameters of the SDE (10) can be chosen to minimise $\sigma^{2}_{f}$. This study was undertaken in DLP (16) for processes of the form (8).

3 Perturbation of Underdamped Langevin Dynamics

The primary objective of this work is to compare the performances of the perturbed underdamped Langevin dynamics (9) and the unperturbed dynamics (4) according to the criteria outlined in Section 2.3 and to find suitable choices for the matrices $J_{1}$, $J_{2}$, $M$ and $\Gamma$ that improve the performance of the sampler. We begin our investigations of (9) by establishing ergodicity and exponentially fast return to equilibrium, and by studying the overdamped limit of (9). As the latter turns out to be nonreversible and therefore in principle superior to the usual overdamped limit (3),e.g. HHMS (05), this calculation provides us with further motivation to study the proposed dynamics.

For the bulk of this work, we focus on the particular case when the target measure is Gaussian, i.e. when the potential is given by $V(q)=\frac{1}{2}q^{T}Sq$ with a symmetric and positive definite precision matrix $S$ (i.e. the covariance matrix is given by $S^{-1}$). In this case, we advocate the following conditions for the choice of parameters:

[TABLE]

Under the above choices (35), we show that the large perturbation limit $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}$ exists and is finite and we provide an explicit expression for it (see Theorem 5). From this expression, we derive an algorithm for finding optimal choices for $J_{1}$ in the case of quadratic observables (see Algorithm 2).

If the friction coefficient is not too small ($\gamma>\sqrt{2}$), and under certain mild nondegeneracy conditions, we prove that adding a small perturbation will always decrease the asymptotic variance for observables of the form $f(q)=q\cdot Kq+l\cdot q+C$:

[TABLE]

see Theorem 4.1. In fact, we conjecture that this statement is true for arbitrary observables $f\in L^{2}(\pi)$, but we have not been able to prove this. The dynamics (9) (used in conjunction with the conditions (35a)-(35c)) proves to be especially effective when the observable is antisymmetric (i.e. when it is invariant under the substitution $q\mapsto-q$) or when it has a significant antisymmetric part. In particular, in Proposition 5 we show that under certain conditions on the spectrum of $J_{1}$, for any antisymmetric observable $f\in L^{2}(\pi)$ it holds that $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}=0$.

Numerical experiments and analysis show that departing significantly from 35c in fact possibly decreases the performance of the sampler. This is in stark contrast to (8), where it is not possible to increase the asymptotic variance by any perturbation. For that reason, until now it seems practical to use (9) as a sampler only when a reasonable estimate of the global covariance of the target distribution is available. In the case of Bayesian inverse problems and diffusion bridge sampling, the target measure $\pi$ is given with respect to a Gaussian prior. We demonstrate the effectiveness of our approach in these applications, taking the prior Gaussian covariance as $S$ in (35a)-(35c).

Remark 3.

In (LNP, 13, Rem. 3) another modification of (4) was suggested (albeit with the simplifications $\Gamma=\gamma\cdot I$ and $M=I$):

[TABLE]

$J$ again denoting an antisymmetric matrix. However, under the change of variables $p\mapsto(1+J)\tilde{p}$ the above equations transform into

[TABLE]

where $\tilde{M}=(1+J)^{-1}M(1-J)^{-1}$ and $\tilde{\Gamma}=(1+J)^{-1}\Gamma(1-J)^{-1}$. Since any observable $f$ depends only on $q$ (the $p$-variables are merely auxiliary), the estimator $\pi_{T}(f)$ as well as its associated convergence characteristics (i.e. asymptotic variance and speed of convergence to equilibrium) are invariant under this transformation. Therefore, (36) reduces to the underdamped Langevin dynamics (4) and does not represent an independent approach to sampling. Suitable choices of $M$ and $\Gamma$ will be discussed in Section 4.5.

3.1 Properties of Perturbed Underdamped Langevin Dynamics

In this section we study some of the properties of the perturbed underdamped dynamics (9). First, note that its generator is given by

[TABLE]

decomposed into the perturbation $\mathcal{L}_{pert}$ and the unperturbed operator $\mathcal{L}_{0}$, which can be further split into the Hamiltonian part $\mathcal{L}_{ham}$ and the thermostat (Ornstein-Uhlenbeck) part $\mathcal{L}_{therm}$, see Pav (14); LRS (10); LS (16).

Lemma 3.

The infinitesimal generator $\operatorname{\mathcal{L}}$ (37) is hypoelliptic.

Proof.

See Appendix B.

An immediate corollary of this result and of Theorem 2.1 is that the perturbed underdamped Langevin process (9) is ergodic with unique invariant distribution $\widehat{\pi}$ given by (5).

As explained in Section 2.3, the exponential decay estimate (27) is crucial for our approach, as in particular it guarantees the well-posedness of the Poisson equation (28). From now on, we will therefore make the following assumption on the potential $V,$ required to prove exponential decay in $L^{2}(\pi)$:

Assumption 1**.**

Assume that the Hessian of $V$ is bounded and that the target measure $\pi(\mathrm{d}q)=\frac{1}{Z}e^{-V}\mathrm{d}q$ satisfies a Poincare inequality, i.e. there exists a constant $\rho>0$ such that

[TABLE]

*holds for all $\phi\in L_{0}^{2}(\pi)\cap H^{1}(\pi)$. *

Sufficient conditions on the potential so that Poincaré’s inequality holds, e.g. the Bakry-Emery criterion, are presented in BGL (13).

Theorem 3.1.

Under Assumption 1 there exist constants $C\geq 1$ and $\lambda>0$ such that the semigroup $(P_{t})_{t\geq 0}$ generated by $\operatorname{\mathcal{L}}$ satisfies exponential decay in $L^{2}(\pi)$ as in (27).

Proof.

See Appendix B.

Remark 4.

The proof uses the machinery of hypocoercivity developed in Vil (09). However, it seems likely that using the framework of DMS (15), the assumption on the boundedness of the Hessian of $V$ can be substantially weakened.

3.2 The Overdamped Limit

In this section we develop a connection between the perturbed underdamped Langevin dynamics (9) and the nonreversible overdamped Langevin dynamics (8). The analysis is very similar to the one presented in (LRS, 10, Section 2.2.2) and we will be brief. For convenience in this section we will perform the analysis on the $d$-dimensional torus $\mathbb{T}^{d}\cong(\mathbb{R}/\mathbb{Z})^{d}$, i.e. we will assume $q\in\mathbb{T}^{d}$. Consider the following scaling of (9):

[TABLE]

valid for the small mass/small momentum regime

[TABLE]

Equivalently, those modifications can be obtained from subsituting $\Gamma\rightarrow\epsilon^{-1}\Gamma$ and $t\mapsto\epsilon^{-1}t$, and so in the limit as $\epsilon\rightarrow 0$ the dynamics (39) describes the limit of large friction with rescaled time. It turns out that as $\epsilon\rightarrow 0$, the dynamics (39) converges to the limiting SDE

[TABLE]

The following proposition makes this statement precise.

Proposition 1.

Denote by $(q_{t}^{\epsilon},p_{t}^{\epsilon})$ the solution to (39) with (deterministic) initial conditions $(q_{0}^{\epsilon},p_{0}^{\epsilon})=(q_{init},p_{init})$ and by $q_{t}^{0}$ the solution to (40) with initial condition $q_{0}^{0}=q_{init}.$ For any $T>0$, $(q_{t}^{\epsilon})_{0\leq t\leq T}$ converges to $(q_{t}^{0})_{0\leq t\leq T}$ in $L^{2}(\Omega,C([0,T]),\mathbb{T}^{d})$ as $\epsilon\rightarrow 0$, i.e.

[TABLE]

Remark 5.

By a refined analysis, it is possible to get information on the rate of convergence; see, e.g. PS (03, 05).

The limiting SDE (40) is nonreversible due to the term $-\mu J_{1}\nabla_{q}V(q_{t})\mathrm{d}t$ and also because the matrix $(\nu J_{2}+\Gamma)^{-1}$ is in general neither symmetric nor antisymmetric. This result, together with the fact that nonreversible perturbations of overdamped Langevin dynamics of the form (8) are by now well-known to have improved performance properties, motivates further investigation of the dynamics (9).

Remark 6.

The limit we described in this section respects the invariant distribution, in the sense that the limiting dynamics (40) is ergodic with respect to the measure $\pi(dq)=\frac{1}{Z}e^{-V}\mathrm{d}q.$ To see this, we have to check that (we are using the notation $\nabla$ instead of $\nabla_{q}$)

[TABLE]

where $\mathcal{L}^{\dagger}$ refers to the $L^{2}(\mathbb{R}^{d})$-adjoint of the generator of (40), i.e. to the associated Fokker-Planck operator. Indeed, the term $\nabla\cdot(\mu e^{-V}J_{1}\nabla V)$ vanishes because of the antisymmetry of $J_{1}.$ Therefore, it remains to show that

[TABLE]

i.e. that the matrix $(\nu J_{2}+\Gamma)^{-1}\Gamma(-\nu J_{2}+\Gamma)^{-1}-(\nu J_{2}+\Gamma)^{-1}$ is antisymmetric. Clearly, the first term is symmetric and furthermore it turns out to be equal to the symmetric part of the second term:

[TABLE]

so $\pi$ is indeed invariant under the limiting dynamics (40).

4 Sampling from a Gaussian Distribution

In this section we study in detail the performance of the Langevin sampler (9) for Gaussian target densities, first considering the case of unit covariance. In particular, we study the optimal choice for the parameters in the sampler, the exponential decay rate and the asymptotic variance. We then extend our results to Gaussian target densities with arbitrary covariance matrices.

4.1 Unit covariance - small perturbations

In our study of the dynamics given by (9) we first consider the simple case when $V(q)=\frac{1}{2}|q|^{2}$, i.e. the task of sampling from a Gaussian measure with unit covariance. We will assume $M=I$, $\Gamma=\gamma I$ and $J_{1}=J_{2}=:J$ (so that the $q-$ and $p-$dynamics are perturbed in the same way, albeit posssibly with different strengths $\mu$ and $\nu$). Using these simplifications, (9) reduces to the linear system

[TABLE]

The above dynamics are of Ornstein-Uhlenbeck type, i.e. we can write

[TABLE]

with $X=(q,p)^{T}$,

[TABLE]

[TABLE]

and $(\bar{W}_{t})_{t\geq 0}$ denoting a standard Wiener process on $\mathbb{R}^{2d}$. The generator of (42) is then given by

[TABLE]

We will consider quadratic observables of the form

[TABLE]

with $K\in\mathbb{R}_{sym}^{d\times d}$, $l\in\mathbb{R}^{d}$ and $C\in\mathbb{R}$, however it is worth recalling that the asymptotic variance $\sigma^{2}_{f}$ does not depend on $C$. We also stress that $f$ is assumed to be independent of $p$ as those extra degrees of freedom are merely auxiliary. Our aim will be to study the associated asymptotic variance $\sigma_{f}^{2}$, see equation (30), in particular its dependence on the parameters $\mu$ and $\nu$. This dependence is encoded in the function

[TABLE]

assuming a fixed observable $f$ and perturbation matrix $J$. In this section we will focus on small perturbations, i.e. on the behaviour of the function $\Theta$ in the neighbourhood of the origin. Our main theoretical tool will be the Poisson equation (28), see the proofs in Appendix C. Anticipating the forthcoming analysis, let us already state our main result, showing that in the neighbourhood of the origin, the function $\Theta$ has favourable properties along the diagonal $\mu=\nu$ (note that the perturbation strengths in the first and second line of (46) coincide):

Theorem 4.1.

Consider the dynamics

[TABLE]

with $\gamma>\sqrt{2}$ and an observable of the form $f(q)=q\cdot Kq+l\cdot q+C$. If at least one of the conditions $[J,K]\neq 0$ and $l\notin\ker J$ is satisfied, then the asymptotic variance of the unperturbed sampler is at a local maximum independently of $K$ and $J$ (and $\gamma$, as long as $\gamma>\sqrt{2}$), i.e.

[TABLE]

and

[TABLE]

4.1.1 Purely quadratic observables

Let us start with the case $l=0$, i.e. $f(q)=q\cdot Kq+C$. The following holds:

Proposition 2.

The function $\Theta$ satisfies

[TABLE]

and

[TABLE]

Proof.

See Appendix C.

The above proposition shows that the unperturbed dynamics represents a critical point of $\Theta$, independently of the choice of $K$, $J$ and $\gamma$. In general though, $\operatorname{Hess}\Theta|_{(\mu,\nu)=(0,0)}$ can have both positive and negative eigenvalues. In particular, this implies that an unfortunate choice of the perturbations will actually increase the asymptotic variance of the dynamics (in contrast to the situation of perturbed *overdamped *Langevin dynamics, where any nonreversible perturbation leads to an improvement in asymptotic variance as detailed in HNW (15) and DLP (16)). Furthermore, the nondiagonality of $\operatorname{Hess}\Theta|_{(\mu,\nu)=(0,0)}$ hints at the fact that the interplay of the perturbations $J_{1}$ and $J_{2}$ (or rather their relative strengths $\mu$ and $\nu$) is crucial for the performance of the sampler and, consequently, the effect of these perturbations cannot be satisfactorily studied independently.

Example 1.

Assuming $J^{2}=-I$ and $[J,K]=0$ it follows that

[TABLE]

for all nonzero $K$. Therefore in this case, a small perturbation of $J_{1}$ only or $J_{2}$ only will increase the asymptotic variance, uniformly over all choices of $K$ and $\gamma$.

However, it turns out that it is possible to construct an improved sampler by combining both perturbations in a suitable way. Indeed, the function $\Theta$ can be seen to have good properties along $\mu=\nu$. We set $\mu(s)=s$, $\nu(s):=s$ and compute

[TABLE]

The last inequality follows from

[TABLE]

and

[TABLE]

(both inequalities are proven in the Appendix, Lemma 10), where the last inequality is strict if $[J,K]\neq 0$. Consequently, choosing both perturbations to be of the same magnitude ($\mu=\nu$) and assuring that $J$ and $K$ do not commute always leads to a smaller asymptotic variance, independently of the choice of $K$, $J$ and $\gamma$. We state this result in the following corrolary:

Corollary 1.

Consider the dynamics

[TABLE]

and a quadratic observable $f(q)=q\cdot Kq+C$. If $[J,K]\neq 0,$ then the asymptotic variance of the unperturbed sampler is at a local maximum independently of K, $J$ and $\gamma$, i.e.

[TABLE]

and

[TABLE]

Remark 7.

As we will see in Section 4.3, more precisely Example 3, if $[J,K]=0,$ the asymptotic variance is constant as a function of $\mu$, i.e. the perturbation has no effect.

Example 2.

Let us set $\mu(s):=s$ and $\nu(s):=-s$ (this corresponds to a small perturbation with $J\nabla V(q_{t})\mathrm{d}t$ in $q$ and $-Jp_{t}\mathrm{d}t$ in $p$). In this case we get

[TABLE]

which changes its sign depending on $J$ and $K$ as the first term is negative and the second is positive. Whether the perturbation improves the performance of the sampler in terms of asymptotic variance therefore depends on the specifics of the observable and the perturbation in this case.

4.1.2 Linear observables

Here we consider the case $K=0$, i.e. $f(q)=l\cdot q+C$, where again $l\in\mathbb{R}^{d}$ and $C\in\mathbb{R}$. We have the following result:

Proposition 3.

The function $\Theta$ satisfies

[TABLE]

and

[TABLE]

Proof.

See Appendix C.

Let us assume that $l\notin\ker J$. Then $\partial_{\mu}^{2}\Theta|_{\mu,\nu=0}<0$, and hence Theorem 3 shows that a small perturbation by $\mu J\nabla V(q_{t})\mathrm{d}t$ alone always results in an improvement of the asymptotic variance. However, if we combine both perturbations $\mu J\nabla V(q_{t})\mathrm{d}t$ and $\nu Jp_{t}\mathrm{d}t$, then the effect depends on the sign of

[TABLE]

This will be negative if $\mu$ and $\nu$ have different signs, and also if they have the same sign and $\gamma$ is big enough. Following Section 4.1.1, we require $\mu=\nu$. We then end up with the requirement

[TABLE]

which is satisfied if $\gamma>\sqrt{2}$

Summarizing the results of this section, for observables of the form $f(q)=q\cdot Kq+l\cdot q+C$, choosing equal perturbations ($\mu=\nu$) with a sufficiently strong damping $(\gamma>\sqrt{2}$) always leads to an improvement in asymptotic variance under the conditions $[J,K]\neq 0$ and $l\notin\ker J$. This is finally the content of Theorem 4.1.

Let us illustrate the results of this section by plotting the asymptotic variance as a function of the perturbation strength $\mu$ (see Figure 1), making the choices $d=2$, $l=(1,1)^{T}$,

[TABLE]

The asymptotic variance has been computed according to (114), using (113a) and (113b) from Appendix C. The graphs confirm the results summarized in Corollary 4.1 concerning the asymptotic variance in the neighbourhood of the unperturbed dynamics ($\mu=0$). Additionally, they give an impression of the global behaviour, i.e. for larger values of $\mu$.

Figures 1(a), 1(b) and 1(c) show the asymptotic variance associated with the quadratic observable $f(q)=q\cdot Kq$. In accordance with Corollary 1, the asymptotic variance is at a local maximum at zero perturbation in the case $\mu=\nu$ (see Figure 1(a)). For increasing perturbation strength, the graph shows that it decays monotonically and reaches a limit for $\mu\rightarrow\infty$ (this limiting behaviour will be explored analytically in Section 4.3). If the condition $\mu=\nu$ is only approximately satisfied (Figure 1(b)), our numerical examples still exhibits decaying asymptotic variance in the neighbourhood of the critical point. In this case, however, the asymptotic variance diverges for growing values of the perturbation $\mu$. If the perturbations are opposed ($\mu=-\nu$) as in Example 2, it is possible for certain observables that the unperturbed dynamics represents a global minimum. Such a case is observed in Figure 1(c). In Figures 1(d) and 1(e) the observable $f(q)=l\cdot q$ is considered. If the damping is sufficiently strong ($\gamma>\sqrt{2}$), the unperturbed dynamics is at a local maximum of the asymptotic variance (Figure 1(d)). Furthermore, the asymptotic variance approaches zero as $\mu\rightarrow\infty$ (for a theoretical explanation see again Section 4.3). The graph in Figure 1(e) shows that the assumption of $\gamma$ not being too small cannot be dropped from Corollary 4.1. Even in this case though the example shows decay of the asymptotic variance for large values of $\mu$.

4.2 Exponential decay rate

Let us denote by $\lambda^{*}$ the optimal exponential decay rate in (27), i.e.

[TABLE]

Note that $\lambda^{*}$ is well-defined and positive by Theorem 3.1. We also define the spectral bound of the generator $\operatorname{\mathcal{L}}$ by

[TABLE]

In MPP (02) it is proven that the Ornstein-Uhlenbeck semigroup $(P_{t})_{t\geq 0}$ considered in this section is differentiable (see Proposition 2.1). In this case (see Corollary 3.12 of EN (00)), it is known that the exponential decay rate and the spectral bound coincide, i.e. $\lambda^{*}=s(\operatorname{\mathcal{L}})$, whereas in general only $\lambda^{*}\leq s(\operatorname{\mathcal{L}})$ holds. In this section we will therefore analyse the spectral properties of the generator (45). In particular, this leads to some intuition of why choosing equal perturbations ($\mu=\nu$) is crucial for the performance of the sampler.

In MPP (02) (see also OPPS (12)), it was proven that the spectrum of $\mathcal{L}$ as in (45) in $L^{2}(\widehat{\pi})$ is given by

[TABLE]

Note that $\sigma(\mathcal{L})$ only depends on the drift matrix $B$. In the case where $\mu=\nu$, the spectrum of $B$ can be computed explicitly.

Lemma 4.

Assume $\mu=\nu$. Then the spectrum of $B$ is given by

[TABLE]

Proof.

We will compute \sigma\big{(}B-\frac{\gamma}{2}I\big{)} and then use the identity

[TABLE]

We have

[TABLE]

where $I$ is understood to denote the identity matrix of appropriate dimension. The above quantity is zero if and only if

[TABLE]

or

[TABLE]

Together with (55), the claim follows.

Using formula (53), in Figure 2(a) we show a sketch of the spectrum $\sigma(-\mathcal{L}$) for the case of equal perturbations ($\mu=\nu)$ with the convenient choices $n=1$ and $\gamma=2.$ Of course, the eigenvalue at [math] is associated to the invariant measure since $\sigma(-\mathcal{L})=\sigma(-\mathcal{L}^{\dagger})$ and $\mathcal{L}^{\dagger}\widehat{\pi}=0$, where $\mathcal{L}^{\dagger}$ denotes the Fokker-Planck operator, i.e. the $L^{2}(\mathbb{R}^{2d})$-adjoint of $\mathcal{L}$. The arrows indicate the movement of the eigenvalues as the perturbation $\mu$ increases in accordance with Lemma 4. Clearly, the spectral bound of $\operatorname{\mathcal{L}}$ is not affected by the perturbation. Note that the eigenvalues on the real axis stay invariant under the perturbation. The subspace of $L_{0}^{2}(\widehat{\pi})$ associated to those will turn out to be crucial for the characterisation of the limiting asymptotic variance as $\mu\rightarrow\infty$.

To illustrate the suboptimal properties of the perturbed dynamics when the perturbations are not equal, we plot the spectrum of the drift matrix $\sigma(B)$ in the case when the dynamics is only perturbed by the term $\nu J_{2}p\mathrm{d}t$ (i.e. $\mu=0$) for $n=2$, $\gamma=2$ and

[TABLE]

(see Figure 2(b)). Note that the full spectrum $\sigma(-\mathcal{L})$ can be inferred from (53). For $\nu=0$ we have that the spectrum $\sigma(B)$ only consists of the (degenerate) eigenvalue $1$. For increasing $\nu$, the figure shows that the degenerate eigenvalue splits up into four eigenvalues, two of which get closer to the imaginary axis as $\nu$ increases, leading to a smaller spectral bound and therefore to a decrease in the speed of convergence to equilibrium. Figures (2(a)) and (2(b)) give an intuitive explanation of why the fine-tuning of the perturbation strengths is crucial.

4.3 Unit covariance - large perturbations

In the previous subsection we observed that for the particular perturbation $J_{1}=J_{2}$ and $\mu=\nu$, i.e.

[TABLE]

the perturbed Langevin dynamics demonstrated an improvement in performance for $\mu$ in a neighbourhood of [math], when the observable is linear or quadratic. Recall that this dynamics is ergodic with respect to a standard Gaussian measure $\widehat{\pi}$ on $\mathbb{R}^{2d}$ with marginal $\pi$ with respect to the $q$–variable. In the following we shall consider only observables that do not depend on $p$. Moreover, we assume without loss of generality that $\pi(f)=0$. For such an observable we will write $f\in L^{2}_{0}(\pi)$ and assume the canonical embedding $L^{2}_{0}(\pi)\subset L^{2}(\widehat{\pi})$. The infinitesimal generator of (57) is given by

[TABLE]

where we have introduced the notation $\mathcal{L}_{pert}=\mu\mathcal{A}$. In the sequel, the adjoint of an operator $B$ in $L^{2}(\widehat{\pi})$ will be denoted by $B^{*}$. In the rest of this section we will make repeated use of the Hermite polynomials

[TABLE]

invoking the notation $x=(q,p)\in\mathbb{R}^{2d}$. For $m\in\mathbb{N}_{0}$ define the spaces

[TABLE]

with induced scalar product

[TABLE]

The space $(H_{m},\langle\cdot,\cdot\rangle_{m})$ is then a real Hilbert space with (finite) dimension

[TABLE]

The following result (Theorem 4.2) holds for operators of the form

[TABLE]

where the quadratic drift and diffusion matrices $B$ and $Q$ are such that $\mathcal{L}$ is the generator of an ergodic stochastic process (see (AE, 14, Definition 2.1) for precise conditions on $B$ and $Q$ that ensure ergodicity). The generator of the SDE (57) is given by (60) with $B$ and $Q$ as in equations (43) and (44), respectively. The following result provides an orthogonal decomposition of $L^{2}(\widehat{\pi})$ into invariant subspaces of the operator $\mathcal{L}$.

Theorem 4.2.

(AE, 14, Section 5)**. The following holds:

1. (a)

The space $L^{2}(\widehat{\pi})$ has a decomposition into mutually orthogonal subspaces:

[TABLE] 2. (b)

For all $m\in\mathbb{N}_{0}$, $H_{m}$ is invariant under $\mathcal{L}$ as well as under the semigroup $(e^{-t\mathcal{L}})_{t\geq 0}$. 3. (c)

The spectrum of $\mathcal{L}$ has the following decomposition:

[TABLE]

where

[TABLE]

Remark 8.

Note that by the ergodicity of the dynamics, $\ker\mathcal{L}$ consists of constant functions and so $\ker\mathcal{L}=H_{0}$. Therefore, $L^{2}_{0}(\widehat{\pi})$ has the decomposition

[TABLE]

Our first main result of this section is an expression for the asymptotic variance in terms of the unperturbed operator $\mathcal{L}_{0}$ and the perturbation $\mathcal{A}$:

Proposition 4.

Let $f\in L_{0}^{2}(\pi)$ (so in particular $f=f(q)$). Then the associated asymptotic variance is given by

[TABLE]

Remark 9.

The proof of the preceding Proposition will show that $\mathcal{L}_{0}^{2}+\mu^{2}\mathcal{A}^{*}\mathcal{A}$ is invertible on $L^{2}_{0}(\widehat{\pi})$ and that $(\mathcal{L}_{0}^{2}+\mu^{2}\mathcal{A}^{*}\mathcal{A})^{-1}f\in\mathcal{D}(\mathcal{L}_{0})$ for all $f\in L^{2}_{0}(\widehat{\pi})$.

To prove Proposition 4 we will make use of the generator with reversed perturbation

[TABLE]

and the momentum flip operator

[TABLE]

Clearly, $P^{2}=I$ and $P^{*}=P$. Further properties of $\mathcal{L}_{0}$, $\mathcal{A}$ and the auxiliary operators $\mathcal{L}_{-}$ and $P$ are gathered in the following lemma:

Lemma 5.

For all $\phi,\psi\in C^{\infty}(\mathbb{R}^{2d})\cap L^{2}(\widehat{\pi})$ the following holds:

1. (a)

The generator $\mathcal{L}_{0}$ is symmetric in $L^{2}(\widehat{\pi})$ with respect to $P$:

[TABLE] 2. (b)

The perturbation $\mathcal{A}$ is skewadjoint in $L^{2}(\widehat{\pi})$:

[TABLE] 3. (c)

The operators $\mathcal{L}_{0}$ and $\mathcal{A}$ commute:

[TABLE] 4. (d)

The perturbation $\mathcal{A}$ satisfies

[TABLE] 5. (e)

$\mathcal{L}$* and $\mathcal{L}_{-}$ commute,*

[TABLE]

and the following relation holds:

[TABLE] 6. (f)

The operators $\mathcal{L}$, $\mathcal{L}_{0}$, $\mathcal{L}_{-}$, $\mathcal{A}$ and $P$ leave the Hermite spaces $H_{m}$ invariant.

Remark 10.

The claim (c) in the above lemma is crucial for our approach, which itself rests heavily on the fact that the $q-$ and $p-$perturbations match ($J_{1}=J_{2}$).

of Lemma 5.

To prove (a), consider the following decomposition of $\mathcal{L}_{0}$ as in (37):

[TABLE]

By partial integration it is straightforward to see that

[TABLE]

and

[TABLE]

for all $\phi,\psi\in C^{\infty}(\mathbb{R}^{2d})\cap L^{2}(\widehat{\pi})$, i.e. $\mathcal{L}_{ham}$ and $\mathcal{L}_{therm}$ are antisymmetric and symmetric in $L^{2}(\widehat{\pi})$ respectively. Furthermore, we immediately see that $P\mathcal{L}_{ham}P\phi=-\mathcal{L}_{ham}\phi$ and $P\mathcal{L}_{therm}P\phi=\mathcal{L}_{therm}\phi$, so that

[TABLE]

We note that this result holds in the more general setting of Section 3 for the infinitesimal generator (37). The claim (b) follows by noting that the flow vector field $b(q,p)=(-Jq,-Jp)$ associated to $\mathcal{A}$ is divergence-free with respect to $\widehat{\pi}$, i.e. $\nabla\cdot(\widehat{\pi}b)=0$. Therefore, $\mathcal{A}$ is the generator of a strongly continuous unitary semigroup on $L^{2}(\widehat{\pi})$ and hence skewadjoint by Stone’s Theorem. To prove (c) we use the decomposition $\mathcal{L}_{0}=\mathcal{L}_{ham}+\mathcal{L}_{therm}$ to obtain

[TABLE]

The first term of (64) gives

[TABLE]

The second term of (64) gives

[TABLE]

since $Jq\cdot\nabla_{q}$ commutes with $p\cdot\nabla_{p}+\Delta_{p}$. Both terms in (65) are clearly zero due the antisymmetry of $J$ and the symmetry of the Hessian $D^{2}_{p}\phi$.

The claim (d) follows from a short calculation similar to the proof of (a). To prove (e), note that the fact that $\mathcal{L}$ and $\mathcal{L}_{-}$ commute follows from (c), as

[TABLE]

while the property $\langle\phi,P\mathcal{L}_{0}P\psi\rangle_{L^{2}(\widehat{\pi})}=\langle\mathcal{L}_{-}\phi,\psi\rangle_{L^{2}(\widehat{\pi})}$ follows from properties (a), (b) and (d). Indeed,

[TABLE]

as required. To prove (f) first notice that $\mathcal{L}$, $\mathcal{L}_{0}$ and $\mathcal{L}_{-}$ are of the form (60) and therefore leave the spaces $H_{m}$ invariant by Theorem 4.2. It follows immediately that also $\mathcal{A}$ leaves those spaces invariant. The fact that $P$ leaves the spaces $H_{m}$ invariant follows directly by inspection of (59).

Now we proceed with the proof of Proposition 4:

of Proposition 4.

Since the potential $V$ is quadratic, Assumption 1 clearly holds and thus Lemma 2 ensures that $\mathcal{L}$ and $\mathcal{L}_{-}$ are invertible on $L^{2}_{0}(\widehat{\pi})$ with

[TABLE]

and analogously for $\mathcal{L}_{-}^{-1}$. In particular, the asymptotic variance can be written as

[TABLE]

Due to the respresentation (67) and Theorem 4.2, the inverses of $\mathcal{L}$ and $\mathcal{L}_{-}$ leave the Hermite spaces $H_{m}$ invariant. We will prove the claim from Proposition 4 under the assumption that $Pf=f$ which includes the case $f=f(q)$. For the following calculations we will assume $f\in H_{m}$ for fixed $m\geq 1$. Combining statement (f) with (a) and (e) of Lemma 5 (and noting that $H_{m}\subset C^{\infty}(\mathbb{R}^{2d})\cap L^{2}(\widehat{\pi})$) we see that

[TABLE]

and

[TABLE]

when restricted to $H_{m}$. Therefore, the following calculations are justified:

[TABLE]

where in the third line we have used the assumption $Pf=f$ and in the fourth line the properties $P^{2}=I$, $P^{*}=P$ and equation (68). Since $\mathcal{L}$ and $\mathcal{L}_{-}$ commute on $H_{m}$ according to Lemma 5(e),(f) we can write

[TABLE]

for the restrictions on $H_{m}$, using $\mathcal{L}+\mathcal{L}_{-}=2\mathcal{L}_{0}$. We also have

[TABLE]

since $\mathcal{L}_{0}$ and $\mathcal{A}$ commute. We thus arrive at the formula

[TABLE]

Now since $(\mathcal{L}_{0}^{2}+\mu^{2}\mathcal{A}^{*}\mathcal{A})^{-1}f=(\mathcal{L}\mathcal{L}_{-})^{-1}f\in\mathcal{D}(\mathcal{L}_{0})$ for all $f\in L^{2}(\widehat{\pi})$, it follows that the operator $-\mathcal{L}_{0}(\mathcal{L}_{0}^{2}+\mu^{2}\mathcal{A}^{*}\mathcal{A})^{-1}$ is bounded. We can therefore extend formula (70) to the whole of $L^{2}(\widehat{\pi})$ by continuity, using the fact that $L^{2}_{0}(\widehat{\pi})=\bigoplus_{m\geq 1}H_{m}$.

Applying Proposition 4 we can analyse the behaviour of $\sigma_{f}^{2}$ in the limit of large perturbation strength $\mu\rightarrow\infty$. To this end, we introduce the orthogonal decomposition

[TABLE]

where $Jq\cdot\nabla_{q}$ is understood as an unbounded operator acting on $L_{0}^{2}(\pi)$, obtained as the smallest closed extension of $Jq\cdot\nabla_{q}$ acting on $C^{\infty}_{c}(\mathbb{R}^{d})$. In particular, $\ker(Jq\cdot\nabla_{q})$ is a closed linear subspace of $L^{2}_{0}(\pi)$. Let $\Pi$ denote the $L_{0}^{2}(\pi)$-orthogonal projection onto $\ker(Jq\cdot\nabla_{q})$. We will write $\sigma_{f}^{2}(\mu)$ to stress the dependence of the asymptotic variance on the perturbation strength. The following result shows that for large perturbations, the limiting asymptotic variance is always smaller than the asymptotic variance in the unperturbed case. Furthermore, the limit is given as the asymptotic variance of the projected observable $\Pi f$ for the unperturbed dynamics.

Theorem 4.3.

Let $f\in L_{0}^{2}(\pi)$, then

[TABLE]

Remark 11.

Note that the fact that the limit exists and is finite is nontrivial. In particular, as Figures 1(b) and 1(c) demonstrate, it is often the case that $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}(\mu)=\infty$ if the condition $\mu=\nu$ is not satisfied.

Remark 12.

The projection $\Pi$ onto $\ker(Jq\cdot\nabla_{q})$ can be understood in terms of Figure 2(a). Indeed, the eigenvalues on the real axis (highlighted by diamonds) are not affected by the perturbations. Let us denote by $\tilde{\Pi}$ the projection onto the span of the eigenspaces of those eigenvalues. As $\mu\rightarrow\infty$, the limiting asymptotic variance is given as the asymptotic variance associated to the unperturbed dynamics of the projection $\tilde{\Pi}f$. If we denote by $\Pi_{0}$ the projection of $L^{2}(\widehat{\pi})$ onto $L^{2}_{0}(\pi)$, then we have that $\Pi\Pi_{0}=\Pi_{0}\tilde{\Pi}$.

of Theorem 4.3.

Note that $\mathcal{L}_{0}$ and $\mathcal{A}^{*}\mathcal{A}$ leave the Hermite spaces $H_{m}$ invariant and their restrictions to those spaces commute (see Lemma 5, (b), (c) and (f)). Furthermore, as the Hermite spaces $H_{m}$ are finite-dimensional, those operators have discrete spectrum. As $\mathcal{A}^{*}\mathcal{A}$ is nonnegative self-adjoint, there exists an orthogonal decomposition $L_{0}^{2}(\pi)=\bigoplus_{i}W_{i}$ into eigenspaces of the operator $-\mathcal{L}_{0}(\mathcal{L}_{0}^{2}+\mu^{2}\mathcal{A}^{*}\mathcal{A})^{-1}$, the decomposition $\bigoplus W_{i}$ being finer then $\bigoplus H_{m}$ in the sense that every $W_{i}$ is a subspace of some $H_{m}$. Moreover,

[TABLE]

where $\lambda_{i}\geq 0$ is the eigenvalue of $\mathcal{A}^{*}\mathcal{A}$ associated to the subspace $W_{i}$. Consequently, formula (62) can be written as

[TABLE]

where $f=\sum_{i}f_{i}$ and $f_{i}\in W_{i}$. Let us assume now without loss of generality that $W_{0}=\ker\mathcal{A}^{*}\mathcal{A}$, so in particular $\lambda_{0}=0$. Then clearly

[TABLE]

Now note that $W_{0}=\ker\mathcal{A}^{*}\mathcal{A}=\ker\mathcal{A}$ due to $\ker\mathcal{A}^{*}=(\operatorname{im}\mathcal{A})^{\perp}$. It remains to show that $\sigma_{\Pi f}^{2}(0)\leq\sigma_{f}^{2}(0)$. To see this, we write

[TABLE]

where

[TABLE]

Note that since we only consider observables that do not depend on $p$, $\Pi f\in\ker(Jq\cdot\nabla_{q})$ and $(1-\Pi)f\in\bigoplus_{i\geq 1}W_{i}$. Since $\mathcal{L}_{0}$ commutes with $\mathcal{A}$, it follows that $(-\mathcal{L}_{0})^{-1}$ leaves both $W_{0}$ and $\bigoplus_{i\geq 1}W_{i}$ invariant. Therefore, as the latter spaces are orthogonal to each other, it follows that $R=0$, from which the result follows.

From Theorem 4.3 it follows that in the limit as $\mu\rightarrow\infty$, the asymptotic variance $\sigma_{f}^{2}(\mu)$ is not decreased by the perturbation if $f\in\ker(Jq\cdot\nabla_{q})$. In fact, this result also holds true non-asymptotically, i.e. observables in $\ker(Jq\cdot\nabla_{q})$ are not affected at all by the perturbation:

Lemma 6.

Let $f\in\ker(Jq\cdot\nabla_{q})$. Then

[TABLE]

for all $\mu\in\mathbb{R}$.

Proof.

From $f\in\ker(Jq\cdot\nabla_{q})$ it follows immediately that $f\in\ker\mathcal{A}^{*}\mathcal{A}$. Then the claim follows from the expression (72).

Example 3.

Recall the case of observables of the form $f(q)=q\cdot Kq+l\cdot q+C$ with $K\in\mathbb{R}_{sym}^{d\times d}$, $l\in\mathbb{R}^{d}$ and $C\in\mathbb{R}$ from Section 4.1. If $[J,K]=0$ and $l\in\ker J$, then $f\in\ker(Jq\cdot\nabla_{q})$ as

[TABLE]

From the preceding lemma it follows that $\sigma_{f}^{2}(\mu)=\sigma_{f}^{2}(0)$ for all $\mu\in\mathbb{R},$ showing that the assumption in Theorem 4.1 does not exclude nontrivial cases.

The following result shows that the dynamics (57) is particularly effective for antisymmetric observables (at least in the limit of large perturbations):

Proposition 5.

Let $f\in L_{0}^{2}(\pi)$ satisfy $f(-q)=-f(q)$ and assume that $\ker J=\{0\}$. Furthermore, assume that the eigenvalues of $J$ are rationally independent, i.e.

[TABLE]

with $\lambda_{i}\in\mathbb{R}_{>0}$ and $\sum_{i}k_{i}\lambda_{i}\neq 0$ for all $(k_{1},\ldots,k_{d})\in\mathbb{Z}^{d}\setminus(0,\ldots,0)$. Then $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}(\mu)=0$.

of Proposition 5.

The claim would immediately follow from $f\in\ker(Jq\cdot\nabla)^{\perp}$ according to Theorem 4.3, but that does not seem to be so easy to prove directly. Instead, we again make use of the Hermite polynomials.

Recall from the proof of Proposition 4 that $\operatorname{\mathcal{L}}$ is invertible on $L_{0}^{2}(\widehat{\pi})$ and its inverse leaves the Hermite spaces $H_{m}$ invariant. Consequently, the asymptotic variance of an observable $f\in L_{0}^{2}(\widehat{\pi})$ can be written as

[TABLE]

where $\Pi_{m}:L_{0}^{2}(\widehat{\pi})\rightarrow H_{m}$ denotes the orthogonal projection onto $H_{m}$. From (59) it is clear that $g_{a}$ is symmetric for $|\alpha|$ even and antisymmetric for $|\alpha|$ odd. Therefore, from $f$ being antisymmetric it follows that

[TABLE]

In view of (54), (61) and (73) the spectrum of $\mathcal{L}_{|H_{m}}$ can be written as

[TABLE]

with appropriate real constants $C_{\alpha,\gamma}\in\mathbb{R}$ that depend on $\alpha$ and $\gamma$, but not on $\mu$. For $|\alpha|=\sum_{j=1}^{2d}\alpha_{j}=m$ odd, we have that

[TABLE]

Indeed, assume to the contrary that the above expression is zero. Then it follows that $\alpha_{j}=\alpha_{j+d}$ for all $j=1,\ldots,d$ by rational independence of $\lambda_{1},\ldots,\lambda_{d}$. From (75a) and (76) it is clear that

[TABLE]

where $B(0,r)$ denotes the ball of radius $r$ centered at the origin in $\mathbb{C}$. Consequently, the spectral radius of $(-\mathcal{L}|_{H_{m}})^{-1}$ and hence $(-\mathcal{L}|_{H_{m}})^{-1}$ itself converge to zero as $\mu\rightarrow\infty$. The result then follows from (74b).

Remark 13.

The idea of the preceding proof can be explained using Figure 2(a) and Remark 12. Since the real eigenvalues correspond to Hermite polynomials of even order, antisymmetric observables are orthogonal to the associated subspaces. The rational independence condition on the eigenvalues of $J$ prevents cancellations that would lead to further eigenvalues on the real axis.

The following corollary gives a version of the converse of Proposition 5 and provides further intuition into the mechanics of the variance reduction achieved by the perturbation.

Corollary 2.

Let $f\in L_{0}^{2}(\pi)$ and assume that $lim_{\mu\rightarrow\infty}\sigma_{f}^{2}(\mu)=0$. Then

[TABLE]

for all $r\in(0,\infty)$, where $B(0,r)$ denotes the ball centered at [math] with radius $r$.

Proof.

According to Theorem 4.3, $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}(\mu)=0$ implies $\sigma_{\Pi f}^{2}(0)=0$. We can write

[TABLE]

and recall from the proof of Proposition 4 that $(-\mathcal{L}_{0})^{-1}$ and $(-\mathcal{L}^{*}_{0})^{-1}$ leave the Hermite spaces $H_{m}$ invariant. Therefore

[TABLE]

in $L^{2}_{0}(\widehat{\pi})$, and in particular $\sigma_{\Pi f}^{2}(0)=0$ implies $\Pi f=0$, which in turn shows that $f\in\ker(Jq\cdot\nabla)^{\perp}$. Using $\ker(Jq\cdot\nabla)^{\perp}=\overline{\operatorname{im}(Jq\cdot\nabla)}$, it follows that there exists a sequence $(\phi_{n})_{n}\in C_{c}^{\infty}(\mathbb{R}^{d})$ such that $Jq\cdot\nabla\phi_{n}\rightarrow f$ in $L^{2}(\pi)$. Taking a subsequence if necessary, we can assume that the convergence is pointwise $\pi$-almost everywhere and that the sequence is pointwise bounded by a function in $L^{1}(\pi)$. Since $J$ is antisymmetric, we have that $Jq\cdot\nabla\phi_{n}=\nabla\cdot(\phi_{n}Jq)$. Now Gauss’s theorem yields

[TABLE]

where $n$ denotes the outward normal to the sphere $\partial B(0,r)$. This quantity is zero due to the orthogonality of $Jq$ and $n$, and so the result follows from Lebesgue’s dominated convergence theorem.

4.4 Optimal Choices of $J$ for Quadratic Observables

Assume $f\in L_{0}^{2}(\pi)$ is given by $f(q)=q\cdot Kq+l\cdot q-\operatorname{Tr}K$, with $K\in\mathbb{R}_{sym}^{d\times d}$ and $l\in\mathbb{R}^{d}$ (note that the constant term is chosen such that $\pi(f)=0$). Our objective is to choose $J$ in such a way that $\lim_{\mu\rightarrow\infty}\sigma_{f}^{2}(\mu)$ becomes as small as possible. To stress the dependence on the choice of $J$, we introduce the notation $\sigma_{f}^{2}(\mu,J)$. Also, we denote the orthogonal projection onto $(\ker J)^{\perp}$ by $\Pi^{\perp}_{\ker J}$.

Lemma 7.

(Zero variance limit for linear observables). Assume $K=0$ and $\Pi^{\perp}_{\ker J}l=0$. Then

[TABLE]

Proof.

According to Proposition 4.3, we have to show that $\Pi f=0$, where $\Pi$ is the $L^{2}(\pi)$-orthogonal projection onto $\ker(Jq\cdot\nabla)$. Let us thus prove that

[TABLE]

where the second identity uses the fact that $(Jq\cdot\nabla)^{*}=-Jq\cdot\nabla$. Indeed, since $\Pi^{\perp}_{\ker J}=0$, by Fredholm’s alternative there exists $u\in\mathbb{R}^{d}$ such that $Ju=l$. Now define $\phi\in L_{0}^{2}(\pi)$ by $\phi(q)=-u\cdot q,$ leading to

[TABLE]

so the result follows.

Lemma 8.

(Zero variance limit for purely quadratic observables.) Let $l=0$ and consider the decomposition $K=K_{0}+K_{1}$ into the traceless part $K_{0}=K-\frac{\operatorname{Tr}K}{d}\cdot I$ and the trace-part $K_{1}=\frac{\operatorname{Tr}K}{d}\cdot I.$ For the corresponding decomposition of the observable

[TABLE]

the following holds:

1. (a)

There exists an antisymmetric matrix $J$ such that $\lim_{\mu\rightarrow\infty}\sigma_{f_{0}}^{2}(\mu,J)=0,$ and there is an algorithmic way (see Algorithm 1) to compute an appropriate $J$ in terms of $K$. 2. (b)

The trace-part is not effected by the perturbation, i.e. $\sigma_{f_{1}}^{2}(\mu,J)=\sigma_{f_{1}}^{2}(0)$ for all $\mu\in\mathbb{R}$.

Proof.

To prove the first claim, according to Theorem 4.3 it is sufficient to show that $f_{0}\in\ker(Jq\cdot\nabla)^{\perp}=\overline{\operatorname{im}(Jq\cdot\nabla)}$. Let us consider the function $\phi(q)=q\cdot Aq$, with $A\in\mathbb{R}_{sym}^{d\times d}$. It holds that

[TABLE]

The task of finding an antisymmetric matrix $J$ such that

[TABLE]

can therefore be accomplished by constructing an antisymmetric matrix $J$ such that there exists a symmetric matrix $A$ with the property $K_{0}=[A,J]$. Given any traceless matrix $K_{0}$ there exists an orthogonal matrix $U\in O(\mathbb{R}^{d})$ such that $UK_{0}U^{T}$ has zero entries on the diagonal, and that $U$ can be obtained in an algorithmic manner (see for example Kaz (88) or (HJ, 13, Chapter 2, Section 2, Problem 3); for the reader’s convenience we have summarised the algorithm in Appendix D.) Assume thus that such a matrix $U\in O(\mathbb{R}^{d})$ has been found and choose real numbers $a_{1},\ldots,a_{d}\in\mathbb{R}$ such that $a_{i}\neq a_{j}$ if $i\neq j$. We now set

[TABLE]

and

[TABLE]

Observe that since $UK_{0}U^{T}$ is symmetric, $\bar{J}$ is antisymmetric. A short calculation shows that $[\bar{A},\bar{J}]=UK_{0}U^{T}$. We can thus define $A=U^{T}\bar{A}U$ and $J=U^{T}\bar{J}U$ to obtain $[A,J]=K_{0}$. Therefore, the $J$ constructed in this way indeed satisfies (79). For the second claim, note that $f_{1}\in\ker(Jq\cdot\nabla)$, since

[TABLE]

because of the antisymmetry of $J$. The result then follows from Lemma 6.

We would like to stress that the perturbation $J$ constructed in the previous lemma is far from unique due to the freedom of choice of $U$ and $a_{1},\ldots,a_{d}\in\mathbb{R}$ in its proof. However, it is asymptotically optimal:

Corollary 3.

In the setting of Lemma 8 the following holds:

[TABLE]

Proof.

The claim follows immediately since $f_{1}\in\ker(Jq\cdot\nabla)$ for arbitrary antisymmetric $J$ as shown in (82), and therefore the contribution of the trace part $f_{1}$ to the asymptotic variance cannot be reduced by any choice of $J$ according to Lemma 6.

As the proof of Lemma 8 is constructive, we obtain the following algorithm for determining optimal perturbations for quadratic observables:

Algorithm 1**.**

Given $K\in\mathbb{R}_{sym}^{d\times d}$, determine an optimal antisymmetric perturbation $J$ as follows:

Set $K_{0}=K-\frac{\operatorname{Tr}K}{d}\cdot I.$ 2. 2.

Find $U\in O(\mathbb{R}^{d})$ such that $UK_{0}U^{T}$ has zero entries on the diagonal (see Appendix D). 3. 3.

Choose $a_{i}\in\mathbb{R},$ $i=1,\ldots d$ such that $a_{i}\neq a_{j}$ for $i\neq j$ and set

[TABLE]

for $i\neq j$ and $\bar{J}_{ii}=0$ otherwise. 4. 4.

Set $J=U^{T}\bar{J}U$.

Remark 14.

In DLP (16), the authors consider the task of finding optimal perturbations $J$ for the nonreversible overdamped Langevin dynamics given in (19). In the Gaussian case this optimization problem turns out be equivalent to the one considered in this section. Indeed, equation (39) of DLP (16) can be rephrased as

[TABLE]

Therefore, Algorithm 1 and its generalization Algorithm 2 (described in Section 4.5) can be used without modifications to find optimal perturbations of overdamped Langevin dynamics.

4.5 Gaussians with Arbitrary Covariance and Preconditioning

In this section we extend the results of the preceding sections to the case when the target measure $\pi$ is given by a Gaussian with arbitrary covariance, i.e. $V(q)=\frac{1}{2}q\cdot Sq$ with $S\in\mathbb{R}_{sym}^{d\times d}$ symmetric and positive definite. The dynamics (9) then takes the form

[TABLE]

The key observation is now that the choices $M=S$ and $\Gamma=\gamma S$ together with the transformation $\widetilde{q}=S^{1/2}q$ and $\widetilde{p}=S^{-1/2}p$ lead to the dynamics

[TABLE]

which is of the form (41) if $J_{1}$ and $J_{2}$ obey the condition $SJ_{1}S=J_{2}$ (note that both $S^{1/2}J_{1}S^{1/2}$ and $S^{-1/2}J_{2}S^{-1/2}$ are of course antisymmetric). Clearly the dynamics (84) is ergodic with respect to a Gaussian measure with unit covariance, in the following denoted by $\widetilde{\pi}$. The connection between the asymptotic variances associated to (83) and (84) is as follows:

For an observable $f\in L_{0}^{2}(\pi)$ we can write

[TABLE]

where $\widetilde{f}(q)=f(S^{-1/2}q)$. Therefore, the asymptotic variances satisfy

[TABLE]

where $\widetilde{\sigma}_{\widetilde{f}}^{2}$ denotes the asymptotic variance of the process $(\widetilde{q}_{t})_{t\geq 0}$. Because of this, the results from the previous sections generalise to (83), subject to the condition that the choices $M=S$, $\Gamma=\gamma S$ and $SJ_{1}S=J_{2}$ are made. We formulate our results in this general setting as corollaries:

Corollary 4.

Consider the dynamics

[TABLE]

with $V(q)=\frac{1}{2}q\cdot Sq$. Assume that $M=S$, $\Gamma=\gamma S$ with $\gamma>\sqrt{2}$ and $SJ_{1}S=J_{2}$. Let $f\in L^{2}(\pi)$ be an observable of the form

[TABLE]

with $K\in\mathbb{R}_{sym}^{d\times d}$, $l\in\mathbb{R}^{d}$ and $C\in\mathbb{R}$. If at least one of the conditions $KJ_{1}S\neq SJ_{1}K$ and $l\notin\ker J$ is satisfied, then the asymptotic variance is at a local maximum for the unperturbed sampler, i.e.

[TABLE]

Proof.

Note that

[TABLE]

is again of the form (87) (where in the last equality, $\widetilde{K}=S^{-1/2}KS^{-1/2}$ and $\widetilde{l}=S^{-1/2}l$ have been defined). From (84), (85) and Theorem 4.1 the claim follows if at least one of the conditions $[\widetilde{K},S^{1/2}J_{1}S^{1/2}]\neq 0$ and $\widetilde{l}\notin\ker(S^{1/2}J_{1}S^{1/2})$ is satisfied. The first of those can easily seen to be equivalent to

[TABLE]

which is equivalent to $KJ_{1}S\neq SJ_{1}K$ since $S$ is nondegenerate. The second condition is equivalent to

[TABLE]

which is equivalent to $J_{1}l\neq 0,$ again by nondegeneracy of $S$.

Corollary 5.

Assume the setting from the previous corollary and denote by $\Pi$ the orthogonal projection onto $\ker(J_{1}Sq\cdot\nabla)$. For $f\in L^{2}(\pi)$ it holds that

[TABLE]

Proof.

Theorem 4.3 implies

[TABLE]

for the transformed system (84). Here $\widetilde{f}(q)=f(S^{-1/2}q)$ is the transformed observable and $\widetilde{\Pi}$ denotes $L^{2}(\pi)$-orthogonal projection onto $\ker(S^{1/2}J_{1}S^{1/2}q\cdot\nabla)$. According to (85), it is sufficient to show that $(\Pi f)\circ S^{-1/2}=\widetilde{\Pi}\widetilde{f}$. This however follows directly from the fact that the linear transformation $\phi\mapsto\phi\circ S^{1/2}$ maps $\ker(S^{1/2}J_{1}S^{1/2}q\cdot\nabla)$ bijectively onto $\ker(J_{1}Sq\cdot\nabla)$.

Let us also reformulate Algorithm 1 for the case of a Gaussian with arbitrary covariance.

Algorithm 2**.**

Given $K,S\in\mathbb{R}_{sym}^{d\times d}$ with $f(q)=q\cdot Kq$ and $V(q)=\frac{1}{2}q\cdot Sq$ (assuming $S$ is nondegenerate), determine optimal perturbations $J_{1}$ and $J_{2}$ as follows:

Set $\widetilde{K}=S^{-1/2}KS^{-1/2}$ and $\widetilde{K}_{0}=\widetilde{K}-\frac{\operatorname{Tr}\widetilde{K}}{d}\cdot I$. 2. 2.

Find $U\in O(\mathbb{R}^{d})$ such that $U\widetilde{K}_{0}U^{T}$ has zero entries on the diagonal (see Appendix D). 3. 3.

Choose $a_{i}\in\mathbb{R}$, $i=1,\ldots,d$ such that $a_{i}\neq a_{j}$ for $i\neq j$ and set

[TABLE] 4. 4.

Set $\widetilde{J}=U^{T}\bar{J}U$. 5. 5.

Put $J_{1}=S^{-1/2}\widetilde{J}S^{-1/2}$ and $J_{2}=S^{1/2}JS^{1/2}$.

Finally, we obtain the following optimality result from Lemma 7 and Corollary 3.

Corollary 6.

Let $f(q)=q\cdot Kq+l\cdot q-\operatorname{Tr}K$ and assume that $\Pi^{\perp}_{\ker J}l=0$. Then

[TABLE]

where $f_{1}(q)=q\cdot K_{1}q$, $K_{1}=\frac{\operatorname{Tr}(S^{-1}K)}{d}S$. Optimal choices for $J_{1}$ and $J_{2}$ can be obtained using Algorithm 2.

Remark 15.

Since in Section 4.1 we analysed the case where $J_{1}$ and $J_{2}$ are proportional, we are not able to drop the restriction $J_{2}=SJ_{1}S$ from the above optimality result. Analysis of completely arbitrary perturbations will be the subject of future work.

Remark 16.

The choices $M=S$ and $\Gamma=\gamma S$ have been introduced to make the perturbations considered in this article lead to samplers that perform well in terms of reducing the asymptotic variance. However, adjusting the mass and friction matrices according to the target covariance in this way (i.e. $M=S$ and $\Gamma=\gamma S$) is a popular way of preconditioning the dynamics, see for instance GC (11) and, in particular mass-tensor molecular dynamics Ben (75). Here we will present an argument why such a preconditioning is indeed beneficial in terms of the convergence rate of the dynamics. Let us first assume that $S$ is diagonal, i.e. $S=\operatorname{diag}(s^{(1)},\ldots,s^{(d)})$ and that $M=\operatorname{diag}(m^{(d)},\ldots,m^{(d)})$ and $\Gamma=\operatorname{diag}(\gamma^{(d)},\ldots,\gamma^{(d)})$ are chosen diagonally as well. Then (83) decouples into one-dimensional SDEs of the following form:

[TABLE]

Let us write those Ornstein-Uhlenbeck processes as

[TABLE]

with

[TABLE]

As in Section 4.2, the rate of the exponential decay of (89) is equal to $\min\text{Re}\,\sigma(B^{(i)})$. A short calculation shows that the eigenvalues of $B^{(i)}$ are given by

[TABLE]

Therefore, the rate of exponential decay is maximal when

[TABLE]

in which case it is given by

[TABLE]

Naturally, it is reasonable to choose $m^{(i)}$ in such a way that the exponential rate $(\lambda^{(i)})^{*}$ is the same for all $i$, leading to the restriction $M=cS$ with $c>0$. Choosing $c$ small will result in fast convergence to equilibrium, but also make the dynamics (88) quite stiff, requiring a very small timestep $\Delta t$ in a discretisation scheme. The choice of $c$ will therefore need to strike a balance between those two competing effects. The constraint (90) then implies $\Gamma=2cS$. By a coordinate transformation, the preceding argument also applies if $S$, $M$ and $\Gamma$ are diagonal in the same basis, and of course $M$ and $\Gamma$ can always be chosen that way. Numerical experiments show that it is possible to increase the rate of convergence to equilibrium even further by choosing $M$ and $\Gamma$ nondiagonally with respect to $S$ (although only by a small margin). A clearer understanding of this is a topic of further investigation.

5 Numerical Experiments: Diffusion Bridge Sampling

5.1 Numerical Scheme

In this section we introduce a splitting scheme for simulating the perturbed underdamped Langevin dynamics given by equation (9). In the unpertubed case, i.e. when $J_{1}=J_{2}=0$, the right-hand side can be decomposed into parts $A$, $B$ and $C$ according to

[TABLE]

i.e. $O$ refers to the Ornstein-Uhlenbeck part of the dynamics, whereas $A$ and $B$ stand for the momentum and position updates, respectively.

One particular splitting scheme which has proven to be efficient is the $BAOAB$ scheme, (see LM (15) and references therein). The string of letters refers to the order in which the different parts are integrated, namely

[TABLE]

We note that many different discretisation schemes such as $ABOBA$, $OABAO$, etc. are viable, but that analytical and numerical evidence has shown that the $BAOAB$-ordering has particularly good properties to compute long-time ergodic averages with respect to $q$-dependent observables. Motivated by this, we introduce the following perturbed scheme, introducing additional Runge-Kutta integration steps between the $A$, $B$ and $O$ parts:

[TABLE]

where $RK_{4}(\Delta t,q_{0})$ refers to fourth order Runge-Kutta integration of the ODE

[TABLE]

up until time $\Delta t$. We remark that the $J_{2}$-perturbation is linear and can therefore be included in the $O$-part without much computational overhead. Clearly, other discretisation schemes are possible as well, for instance one could use a symplectic integrator for the ODE (93), noting that it is of Hamiltonian type. However, since $V$ as the Hamiltonian for (93) is not separable in general, such a symplectic integrator would have to be implcit. Moreover, (92c) and (92e) could be merged since (92e) commutes with (92d). In this paper, we content ourselves with the above scheme for our numerical experiments.

Remark 17.

The aformentioned schemes lead to an error in the approximation for $\pi(f)$, since the invariant measure $\pi$ is not preserved exactly by the numerical scheme. In practice, the $BAOAB$-scheme can therefore be accompanied by an accept-reject Metropolis step as in MWL (16), leading to an unbiased estimate of $\pi(f)$, albeit with an inflated variance. In this case, after every rejection the momentum variable has to be flipped ($p\mapsto-p$) in order to keep the correct invariant measure. We note here that our perturbed scheme can be ’Metropolized’ in a similar way by ’flipping the matrices $J_{1}$ and $J_{2}$ after every rejection ($J_{1}\mapsto-J_{1}$ and $J_{2}\mapsto-J_{2})$ and using an appropriate (volume-preserving and time-reversible) integrator for the dynamics given by (93). Implementations of this idea are the subject of ongoing work.

5.2 Diffusion Bridge Sampling

To numerically test our analytical results, we will apply the dynamics (9) to sample a measure on path space associated to a diffusion bridge. Specifically, consider the SDE

[TABLE]

with $X_{s}\in\mathbb{R}^{n}$, $\beta>0$ and the potential $U:\mathbb{R}^{n}\rightarrow\mathbb{R}$ obeying adequate growth and smoothness conditions (see HSV (07), Section 5 for precise statements). The law of the solution to this SDE conditioned on the events $X(0)=x_{-}$ and $X(s_{+})=x_{+}$ is a probability measure $\pi$ on $L^{2}([0,s_{+}],\mathbb{R}^{n})$ which poses a challenging and important sampling problem, especially if $U$ is multimodal. This setting has been used as a test case for sampling probability measures in high dimensions (see for example BPSSS (11) and OPPS (16)). For a more detailed introduction (including applications) see BS (09) and for a rigorous theoretical treatment the papers HSVW (05); HSV (07, 09); BS (09) .

In the case $U\equiv 0$, it can be shown that the law of the conditioned process is given by a Gaussian measure $\pi_{0}$ with mean zero and precision operator $\mathcal{S}=-\frac{\beta}{2}\Delta$ on the Sobolev space $H^{1}([0,s_{+}],\mathbb{R}^{d})$ equipped with appropriate boundary conditions. The general case can then be understood as a perturbation thereof: The measure $\pi$ is absolutely continuous with respect to $\pi_{0}$ with Radon-Nikodym derivative

[TABLE]

where

[TABLE]

and

[TABLE]

We will make the choice $x_{-}=x_{+}=0$, which is possible without loss of generality as explained in (BRSV, 08, Remark 3.1), leading to Dirichlet boundary conditions on $[0,s_{+}]$ for the precision operator $\mathcal{S}$. Furthermore, we choose $s_{+}=1$ and discretise the ensuing $s$-interval $[0,1]$ according to

[TABLE]

in an equidistant way with stespize $s_{j+1}-s_{j}\equiv\delta=\frac{1}{d+1}$. Functions on this grid are determined by the values $x(s_{1})=x_{1},\ldots,x(s_{n})=x_{n}$, recalling that $x(0)=x(1)=0$ by the Dirichlet boundary conditions. We discretise the functional $\Psi$ as

[TABLE]

such that its gradient is given by

[TABLE]

The discretised version $A$ of the Dirichlet-Laplacian $\Delta$ on $[0,1]$ is given by

[TABLE]

Following (94), the discretised target measure $\widehat{\pi}$ has the form

[TABLE]

with

[TABLE]

In the following we will consider the case $n=1$ with potential $U:\mathbb{R}\rightarrow\mathbb{R}$ given by $U(x)=\frac{1}{2}(x^{2}-1)^{2}$ and set $\beta=1$. To test our algorithm we adjust the parameters $M$, $\Gamma$, $J_{1}$ and $J_{2}$ according to the recommended choice in the Gaussian case,

[TABLE]

where we take $S=\frac{\beta}{2}\delta\cdot A$ as the precision operator of the Gaussian target. We will consider the linear observable $f_{1}(x)=l\cdot x$ with $l=(1,\ldots,1)$ and the quadratic observable $f_{2}(x)=|x|^{2}$. In a first experiment we adjust the perturbation $J_{1}$ (and via (95) also $J_{2}$) to the observable $f_{2}$ according to Algorithm 2. The dynamics (9) is integrated using the splitting scheme introduced in Section 5.1 with a stepsize of $\Delta t=10^{-4}$ over the time interval $[0,T]$ with $T=10^{2}$. Furthermore, we choose initial conditions $q_{0}=(1,\ldots,1)$, $p_{0}=(0,\ldots,0)$ and introduce a burn-in time $T_{0}=1$, i.e. we take the estimator to be

[TABLE]

We compute the variance of the above estimator from $N=500$ realisations and compare the results for different choices of the friction coefficient $\gamma$ and of the perturbation strength $\mu$.

The numerical experiments show that the perturbed dynamics generally outperform the unperturbed dynamics independently of the choice of $\mu$ and $\gamma$, both for linear and quadratic observables. One notable exception is the behaviour of the linear observable for small friction $\gamma=10^{-3}$ (see Figure 3(a)), where the asymptotic variance initially increases for small perturbation strengths $\mu$. However, this does not contradict our analytical results, since the small perturbation results from Section 4.1 generally require $\gamma$ to be sufficiently big (for example $\gamma\geq\sqrt{2}$ in Theorem 4.1). We remark here that the condition $\gamma\geq\sqrt{2}$, while necessary for the theoretical results from Section 4.1, is not a very advisable choice in practice (at least in this experiment), since Figures 3(b) and 4(b) clearly indicate that the optimal friction is around $\gamma\approx 10^{-1}$. Interestingly, the problem of choosing a suitable value for the friction coefficient coefficient $\gamma$ becomes mitigated by the introduction of the perturbation: While the performance of the unperturbed sampler depends quite sensitively on $\gamma$, the asymptotic variance of the perturbed dynamics is a lot more stable with respect to variations of $\gamma$.

In the regime of growing values of $\mu$, the experiments confirm the results from Section 4.3, i.e. the asymptotic variance approaches a limit that is smaller than the asymptotic variance of the unperturbed dynamics.

As a final remark we report our finding that the performance of the sampler for the linear observable is qualitatively independent of the coice of $J_{1}$ (as long as $J_{2}$ is adjusted according to (95)). This result is in alignment with Propostion 5 which predicts good properties of the sampler for antisymmetric observables. In contrast to this, a judicious choice of $J_{1}$ is critical for quadratic observables. In particular, applying Algorithm 2 significantly improves the performance of the perturbed sampler in comparison to choosing $J_{1}$ arbitrarily.

6 Outlook and Future Work

A new family of Langevin samplers was introduced in this paper. These new SDE samplers consist of perturbations of the underdamped Langevin dynamics (that is known to be ergodic with respect to the canonical measure), where auxiliary drift terms in the equations for both the position and the momentum are added, in a way that the perturbed family of dynamics is ergodic with respect to the same (canonical) distribution. These new Langevin samplers were studied in detail for Gaussian target distributions where it was shown, using tools from spectral theory for differential operators, that an appropriate choice of the perturbations in the equations for the position and momentum can improve the performance of the Langvin sampler, at least in terms of reducing the asymptotic variance. The performance of the perturbed Langevin sampler to non-Gaussian target densities was tested numerically on the problem of diffusion bridge sampling.

The work presented in this paper can be improved and extended in several directions. First, a rigorous analysis of the new family of Langevin samplers for non-Gaussian target densities is needed. The analytical tools developed in DLP (16) can be used as a starting point. Furthermore, the study of the actual computational cost and its minimization by an appropriate choice of the numerical scheme and of the perturbations in position and momentum would be of interest to practitioners. In addition, the analysis of our proposed samplers can be facilitated by using tools from symplectic and differential geometry. Finally, combining the new Langevin samplers with existing variance reduction techniques such as zero variance MCMC, preconditioning/Riemannian manifold MCMC can lead to sampling schemes that can be of interest to practitioners, in particular in molecular dynamics simulations. All these topics are currently under investigation.

Acknowledgments

AD was supported by the EPSRC under grant No. EP/J009636/1. NN is supported by EPSRC through a Roth Departmental Scholarship. GP is partially supported by the EPSRC under grants No. EP/J009636/1, EP/L024926/1, EP/L020564/1 and EP/L025159/1. Part of the work reported in this paper was done while NN and GP were visiting the Institut Henri Poincaré during the Trimester Program ”Stochastic Dynamics Out of Equilibrium”. The hospitality of the Institute and of the organizers of the program is greatly acknowledged.

Appendix A Estimates for the Bias and Variance

of Lemma 1.

Suppose that $(P_{t})_{t\geq 0}$ satisfies (27). Let $\pi_{0}$ be an initial distribution of $(X_{t})_{t\geq 0}$ such that $\pi_{0}\ll\pi$ and $h=\frac{d\pi_{0}}{d\pi}\in L^{2}(\pi)$. Slightly abusing notation, we denote by $\pi_{0}P_{t}$ the law of $X_{t}$ given $X_{0}\sim\pi$. Then

[TABLE]

where $P_{t}^{*}$ denotes the $L^{2}(\pi)$-adjoint of $P_{t}$. Since $f$ is assumed to be bounded, we immediately obtain

[TABLE]

and so, for $X_{0}\sim\pi_{0}$,

[TABLE]

as required.

of Lemma 2.

Given $f\in L^{2}(\pi)$, for fixed $T>0$,

[TABLE]

Then we have that $\chi_{T}\in\mathcal{D}(\mathcal{L})$ and $\mathcal{L}\chi_{T}=f-P_{T}f$, moreover

[TABLE]

so that $\{\chi_{T}\}_{T\geq 0}$ is a Cauchy sequence in $L^{2}(\pi)$ converging to $\chi=\int_{0}^{\infty}\left(\pi(f)-P_{t}f\right)\mathrm{d}t$. Since $\mathcal{L}$ is closed and

[TABLE]

in $L^{2}(\pi)$, it follows that $\chi\in\mathcal{D}(\mathcal{L})$ and $\mathcal{L}\chi=f-\pi(f)$. Moreover,

[TABLE]

where $K_{\lambda}=C\int_{0}^{\infty}e^{-\lambda t}\,\mathrm{d}t$. Since we assume that $f$ is smooth, the coefficients are smooth and $\operatorname{\mathcal{L}}$ is hypoelliptic, then $\mathcal{L}\chi=f-\pi(f)$ implies that $\chi\in C^{\infty}(\operatorname{\mathbb{R}}^{d})$, and thus we can apply Itô’s formula to $\chi(X_{t})$ to obtain:

[TABLE]

One can check that the conditions of (EK, 86, Theorem 7.1.4) hold. In particular, the following central limit theorem follows

[TABLE]

By Theorem 2.1, the generator $\operatorname{\mathcal{L}}$ has the form

[TABLE]

where $\nabla\cdot(\pi\gamma)=0$. It follows that

[TABLE]

First suppose that $X_{0}\sim\pi$. Then $(\chi(X_{t}))_{t\geq 0}$ is a stationary process, and so

[TABLE]

From which (29) follows. More generally, suppose that $X_{0}\sim\pi_{0}$, where $\pi_{0}(x)=h(x)\pi(x)$ for $h\in L^{2}(\pi)$. If $f\in L^{\infty}(\pi)$, then by (1),

[TABLE]

so that $\chi\in L^{\infty}(\pi)$. Therefore $\frac{1}{\sqrt{T}}(\chi(X_{0})-\chi(X_{T}))\xrightarrow{p}0$ as $T\rightarrow\infty$, and so (29) holds in this case, similarly.

Appendix B Proofs of Section 3

Proof.

of Lemma 3 We first note that $\operatorname{\mathcal{L}}$ in (37) can be written in the “sum of squares” form:

[TABLE]

where

[TABLE]

and

[TABLE]

Here $\{e_{k}\}_{k=1,\ldots,d}$ denotes the standard Euclidean basis and $\Gamma^{1/2}$ is the unique positive definite square root of the matrix $\Gamma$. The relevant commutators turn out to be

[TABLE]

Because $\Gamma$ has full rank on $\operatorname{\mathbb{R}}^{d}$, it follows that

[TABLE]

Since

[TABLE]

and $\operatorname{span}(\{\Gamma^{1/2}M^{-1}\nabla_{q}\,:\,k=1,\ldots,d\})=\operatorname{span}\{\partial_{q_{k}}:k=1,\ldots,d\}$, it follows that

[TABLE]

so the assumptions of Hörmander’s theorem hold.

B.1 The overdamped limit

The following is a technical lemma required for the proof of Proposition 1:

Lemma 9.

Assume the conditions from Proposition 1. Then for every $T>0$ there exists $C>0$ such that

[TABLE]

Proof.

Using variation of constants, we can write the second line of (39) as

[TABLE]

We then compute

[TABLE]

Clearly, the first term on the right hand side of (98) is bounded. For the second term, observe that

[TABLE]

since $V\in C^{1}(\mathbb{T}^{d})$ and therefore $\nabla_{q}V$ is bounded. By the basic matrix exponential estimate $\|e^{-t(\nu J_{2}+\Gamma)M^{-1}}\|\leq Ce^{-\omega t}$ for suitable $C$ and $\omega$, we see that (99) can further be bounded by

[TABLE]

so this term is bounded as well. The third term is bounded by the Burkholder–Davis–Gundy inequality and a similar argument to the one used for the second term applies. The cross terms can be bounded by the previous ones, using the Cauchy-Schwarz inequality and the elementary fact that $\sup(ab)\leq\sup a\cdot\sup b$ for $a,b>0$, so the result follows.

of Proposition 1.

Equations (39) can be written in integral form as

[TABLE]

and

[TABLE]

where the first line has been multiplied by the matrix $\nu J_{2}+\Gamma$. Combining both equations yields

[TABLE]

Now applying Lemma 9 gives the desired result, since the above equation differs from the integral version of (40) only by the term $\epsilon(\nu J_{2}+\Gamma)^{-1}(p_{t}^{\epsilon}-p_{0})$ which vanishes in the limit as $\epsilon\rightarrow 0$.

B.2 Hypocoercivity

The objective of this section is to prove that the perturbed dynamics (9) converges to equilibrium exponentially fast, i.e. that the associated semigroup $(P_{t})_{t\geq 0}$ satisfies the estimate (27). We we will be using the theory of hypocoercivity outlined in Vil (09) (see also the exposition in (Pav, 14, Section 6.2)). We provide a brief review of the theory of hypocoercivity.

Let $(\mathcal{H},\langle\cdot,\cdot\rangle)$ be a real separable Hilbert space and consider two unbounded operators $A$ and $B$ with domains $D(A)$ and $D(B)$ respectively, $B$ antisymmetric. Let $S\subset\mathcal{H}$ be a dense vectorspace such that $S\subset D(A)\cap D(B)$, i.e. the operations of $A$ and $B$ are authorised on $S$. The theory of hypocoercivity is concerned with equations of the form

[TABLE]

and the associated semigroup $(P_{t})_{t\geq 0}$ generated by $L=A^{*}A-B$. Let us also introduce the notation $K=\ker L$. With the choices $\mathcal{H}=L^{2}(\widehat{\pi})$, $A=\sigma\nabla_{p}$ and $B=M^{-1}p\cdot\nabla_{q}-\nabla_{q}V\cdot\nabla_{p}-\mu J_{1}\nabla_{q}V\cdot\nabla_{q}-\nu J_{2}M^{-1}p\cdot\nabla_{p},$ it turns out that $L$ is the (flat) $L^{2}(\mathbb{R}^{2d})$-adjoint of the generator $\mathcal{L}$ given in (37) and therefore equation (101) is the Fokker-Planck equation associated to the dynamics (9). In many situations of practical interest, the operator $A^{*}A$ is coercive only in certain directions of the state space, and therefore exponential return to equilibrium does not follow in general. In our case for instance, the noise acts only in the $p$-variables and therefore relaxation in the $q$-variables cannot be concluded a priori. However, intuitively speaking, the noise gets transported through the equations by the Hamiltonian part of the dynamics. This is what the theory of hypocoercivity makes precise. Under some conditions on the interactions between $A$ and $B$ (encoded in their iterated commutators), exponential return to equilibrium can be proved. To state the main abstract theorem, we need the following definitions:

Definition 1.

(Coercivity) Let $T$ be an unbounded operator on $\mathcal{H}$ with domain $D(T)$ and kernel $K$. Assume that there exists another Hilbert space $(\tilde{\mathcal{H}},\langle\cdot,\cdot\rangle_{\tilde{\mathcal{H}}})$, continuously and densely embedded in $K^{\perp}$. The operator is said to be $\lambda$-coercive if

[TABLE]

for all $h\in K^{\perp}\cap D(T)$.

Definition 2.

An operator $T$ on $\mathcal{H}$ is said to be relatively bounded with respect to the operators $T_{1},\ldots,T_{n}$ if the intersection of the domains $\cap D(T_{j})$ is contained in $D(T)$ and there exists a constant $\alpha>0$ such that

[TABLE]

holds for all $h\in D(T)$.

We can now proceed to the main result of the theory.

Theorem B.1.

(Vil, 09, Theorem 24)* Assume there exists $N\in\mathbb{N}$ and possibly unbounded operators*

[TABLE]

such that $C_{0}=A$,

[TABLE]

and for all $k=0,1,\ldots,N$

1. (a)

$[A,C_{k}]$* is relatively bounded with respect to $\{C_{j}\}_{0\leq j\leq k}$ and $\{C_{j}A\}_{0\leq j\leq k-1}$,* 2. (b)

$[C_{k},A^{*}]$* is relatively bounded with respect to $I$ and $\{C_{j}\}_{0\leq j\leq k}$ ,* 3. (c)

$R_{k}$* is relatively bounded with respect to $\{C_{j}\}_{0\leq j\leq k-1}$ and $\{C_{j}A\}_{0\leq j\leq k-1}$ and* 4. (d)

there are positive constants $\lambda_{i}$, $\Lambda_{i}$ such that $\lambda_{j}I\leq Z_{j}\leq\Lambda_{j}I$.

Furthermore, assume that $\sum_{j=0}^{N}C_{j}^{*}C_{j}$ is $\kappa$-coercive for some $\kappa>0$. Then, there exists $C\geq 0$ and $\lambda>0$ such that

[TABLE]

where $\mathcal{H}^{1}\subset\mathcal{H}$ is the subspace associated to the norm

[TABLE]

and $K=\ker(A^{*}A-B)$.

Remark 18.

Property (103) is called hypocoercivity of $L$ on $\mathcal{H}^{1}:=(K^{\perp},\|\cdot\|_{\mathcal{H}^{1}})$.

If the conditions of the above theorem hold, we also get a regularization result for the semigroup $e^{-tL}$ (see (Vil, 09, Theorem A.12)):

Theorem B.2.

Assume the setting and notation of Theorem B.1. Then there exists a constant $C>0$ such that for all $k=0,1,\ldots,N$ and $t\in(0,1]$ the following holds:

[TABLE]

of Theorem 3.1.

. We pove the claim by verifying the conditions of Theorem B.1. Recall that $C_{0}=A=\sigma\nabla_{p}$ and

[TABLE]

A quick calculation shows that

[TABLE]

so that indeed

[TABLE]

and

[TABLE]

We make the choice $N=1$ and calculate the commutator

[TABLE]

Let us now set $C_{1}=\sigma M^{-1}\nabla_{q}$, $Z_{1}=1$ and $R_{1}=\nu\sigma M^{-1}J_{2}\nabla_{p}$, such that (102) holds for $j=0$. Note that $[A,A]=0$111This is not true automatically, since $[A,A]$ stands for the array $([A_{j},A_{k}])_{jk}$., $[A,C_{1}]=0$ and $[A^{*},C_{1}]=0$. Furthermore, we have that

[TABLE]

We now compute

[TABLE]

and choose $R_{2}=[C_{1},B]$, $Z_{2}=1$ and recall that $C_{2}=0$ by assumption (of Theorem B.1). With those choices, assumptions (a)-(d) of Theorem B.1 are fulfilled. Indeed, assumption (a) holds trivially since all relevant commutators are zero. Assumption (b) follows from the fact that $[A,A^{*}]=\sigma M^{-1}\sigma$ is clearly bounded relative to $I$. To verify assumption (c), let us start with the case $k=1$. It is necessary to show that $R_{1}=\nu\sigma M^{-1}J_{2}\nabla_{p}$ is bounded relatively to $A=\sigma\nabla_{p}$ and $A^{2}$. This is obvious since the $p$-derivatives appearing in $R_{1}$ can be controlled by the $p$-derivatives appearing in $A$. For $k=2,$ a similar argument shows that $R_{2}=-\sigma M^{-1}\nabla^{2}V\nabla_{p}+\mu\sigma M^{-1}\nabla^{2}VJ_{1}\nabla_{q}$ is bounded relatively to $A=\sigma\nabla_{p}$ and $C_{1}=\sigma M^{-1}\nabla_{q}$ because of the assumption that $\nabla^{2}V$ is bounded. Note that it is crucial for the preceding arguments to assume that the matrices $\sigma$ and $M$ have full rank. Assumption (d) is trivially satisfied, since $Z_{1}$ and $Z_{2}$ are equal to the identity. It remains to show that

[TABLE]

is $\kappa$-coercive for some $\kappa>0$. It is straightforward to see that the kernel of $T$ consists of constant functions and therefore

[TABLE]

Hence, $\kappa$-coercivity of $T$ amounts to the functional inequality

[TABLE]

Since the transformation $\phi\mapsto\psi$, $\psi(q,p)=\phi(\sigma^{-1}Mq,\sigma^{-1}p)$ is bijective on $H^{1}(\mathbb{R}^{2d},\widehat{\pi})$, the above is equivalent to

[TABLE]

i.e. a Poincaré inequality for $\widehat{\pi}$. Since $\widehat{\pi}=\pi\otimes\mathcal{N}(0,M),$ coercivity of $T$ boils down to a Poincaré inequality for $\pi$ as in Assumption 1. This concludes the proof of the hypocoercive decay estimate (103). Clearly, the abstract $\mathcal{H}^{1}$-norm from (104) is equivalent to the Sobolev norm $H^{1}(\widehat{\pi})$, and therefore it follows that there exist constants $C\geq 0$ and $\lambda\geq 0$ such that

[TABLE]

for all $f\in H^{1}(\widehat{\pi})\setminus K$, where $K=\ker T$ consists of constant functions. Let us now lift this estimate to $L^{2}(\widehat{\pi})$. There exist a constant $\tilde{C}\geq 0$ such that

[TABLE]

Therefore, Theorem B.2 implies

[TABLE]

for $t=1$ and a possibly different constant $\tilde{C}$. Let us now assume that $t\geq 1$ and $f\in L^{2}(\widehat{\pi})\setminus K$. It holds that

[TABLE]

where the last inequality follows from (105). Now applying (107) and gathering constants results in

[TABLE]

Note that although we assumed $t\geq 1$, the above estimate also holds for $t\geq 0$ (although possibly with a different constant $C$) since $\|P_{t}\|_{L^{2}(\widehat{\pi})\rightarrow L^{2}(\widehat{\pi})}$ is bounded on $[0,1]$.

Appendix C Asymptotic Variance of Linear and Quadratic Observables in the Gaussian Case

We begin by deriving a formula for the asymptotic variance of observables of the form

[TABLE]

with $K\in\mathbb{R}_{sym}^{d\times d}$ and $l\in\mathbb{R}^{d}$. Note that the constant term is chosen such that $\widehat{\pi}(f)=0$. The following calculations are very much along the lines of (DLP, 16, Section 4). Since the Hessian of $V$ is bounded and the target measure $\pi$ is Gaussian, Assumption 1 is satisfied and exponential decay of the semigroup $(P_{t})_{t\geq 0}$ as in (27) follows by Theorem 3.1. According to Lemma 2, the asymptotic variance is then given by

[TABLE]

where $\chi$ is the solution to the Poisson equation

[TABLE]

Recall that

[TABLE]

is the generator as in (45), where for later convenience we have defined $A=B^{T}$, i.e.

[TABLE]

In the sequel we will solve (111) analytically. First, we introduce the notation

[TABLE]

and

[TABLE]

such that by slight abuse of notation $f$ is given by

[TABLE]

By uniqueness (up to a constant) of the solution to the Poisson equation (111) and linearity of $\mathcal{L}$, $g$ has to be a quadratic polynomial, so we can write

[TABLE]

where $C\in\mathbb{R}_{sym}^{2d\times 2d}$ and $D\in\mathbb{R}^{2d}$ (notice that $C$ can be chosen to be symmetrical since $x\cdot Cx$ does not depend on the antisymmetric part of $C$). Plugging this ansatz into (111) yields

[TABLE]

where

[TABLE]

denotes the trace of the momentum component of $C$. Comparing different powers of $x$, this leads to the conditions

[TABLE]

Note that (113c) will be satisfied eventually by existence and uniqueness of the solution to (111). Then, by the calculations in DLP (16), the asymptotic variance is given by

[TABLE]

of Proposition 2.

. According to (114) and (113a), the asymptotic variance satisfies

[TABLE]

where the matrix $C$ solves

[TABLE]

and $A$ is given as in (112). We will use the notation

[TABLE]

and the abbreviations $C(0):=C(0,0)$, $C^{\mu}(0):=\partial_{\mu}C|_{\mu,\nu=0}$ and $C^{\nu}(0):=\partial_{\nu}C|_{\mu,\nu=0}$. Let us first determine $C(0)$, i.e. the solution to the equation

[TABLE]

This leads to the following system of equations,

[TABLE]

Note that equations (116b) and (116c) are equivalent by taking the transpose. Plugging (116a) into (116e) yields

[TABLE]

Adding (116b) and (116c), together with (116a) and (117) leads to

[TABLE]

Solving (116b) we obtain,

[TABLE]

so that

[TABLE]

Taking the $\mu$-derivative of (115) and setting $\mu=\nu=0$ yields

[TABLE]

Notice that

[TABLE]

With computations similar to those in the derivation of (118) (or by simple substitution), equation (119) can be solved by

[TABLE]

We employ a similar strategy to determine $C^{\nu}(0)$: Taking the $\nu$-derivative in equation (115), setting $\mu=\nu=0$ and inserting $C(0)$ and $A(0)$ as in (118) and (112) leads to the equation

[TABLE]

which can be solved by

[TABLE]

Note that $\operatorname{Tr}(C\bar{K})=\operatorname{Tr}(C_{1}K)$, and so

[TABLE]

since clearly $\operatorname{Tr}([K,J],K)=\operatorname{Tr}(KJK)-\operatorname{Tr}(JK^{2})=0$. In the same way it follows that

[TABLE]

proving (47).

Taking the second $\mu$-derivative of (115) and setting $\mu=\nu=0$ yields

[TABLE]

employing the notation $C^{\mu\mu}(0)=\partial_{\mu}^{2}C|_{\mu,\nu=0}$ and noticing that $\partial_{\mu}^{2}A=0$. Using (120) we calculate

[TABLE]

As before, we make the ansatz

[TABLE]

leading to the equations

[TABLE]

Again, (122b) and (122c) are equivalent by taking the transpose. Plugging (122a) into (122d) and combing with (122b) or (122c) gives

[TABLE]

Now

[TABLE]

gives the first part of (48). We proceed in the same way to determine $C_{1}^{\nu\nu}(0)$. Analogously, we get

[TABLE]

Solving the resulting linear matrix system (similar to (122a)-(122d)) results in

[TABLE]

leading to

[TABLE]

To compute the cross term $C_{1}^{\mu\nu}(0)$ we take the mixed derivative $\partial_{\mu\nu}^{2}$ of (115) and set $\mu=\nu=0$ to arrive at

[TABLE]

Using $\eqref{eq:C^mu}$ and (121) we see that

[TABLE]

The ensuing linear matrix system yields the solution

[TABLE]

leading to

[TABLE]

This completes the proof.

Proof of Proposition 3.

By (113b) and (114) the function $\Theta$ satisfies

[TABLE]

Recall the following formula for blockwise inversion of matrices using the Schur complement:

[TABLE]

provided that $X$ and $U-VX^{-1}W$ are invertible. Using this, we obtain

[TABLE]

Taking derivatives, setting $\mu=\nu=0$ and using the fact that $J^{T}=-J$ leads to the desired result.

Lemma 10.

The following holds:

1. (a)

$\gamma-\frac{4}{\gamma^{3}}-\gamma^{3}-\frac{1}{\gamma}<0$* for $\gamma\in(0,\infty)$.* 2. (b)

Let $J=-J^{T}$ and $K=K^{T}$. Then $\operatorname{Tr}(JKJK)-\operatorname{Tr}(J^{2}K^{2})\geq 0$. Furthermore, equality holds if and only if $[J,K]=0.$

Proof.

To show (a) we note that $\gamma-\frac{4}{\gamma^{3}}-\gamma^{3}-\frac{1}{\gamma}<\gamma-\frac{4}{\gamma^{3}}-\gamma^{3}=\gamma(1-\frac{4}{\gamma^{4}}-\gamma^{2})$. The function $f(\gamma):=1-\frac{4}{\gamma^{4}}-\gamma^{2}$ has a unique global maximum on $(0,\infty)$ at $\gamma_{min}=8^{1/6}$ with $f(\gamma_{min})=-2$, so the result follows.

For (b) we note that $[J,K]^{T}=[J,K],$ and that $[J,K]^{2}$ is symmetric and nonnegative definite. We can write

[TABLE]

with $\lambda_{i}$ denoting the (real) eigenvalues of $[J,K]$. From this it follows that $\operatorname{Tr}([J,K]^{2})\geq 0$ with equality if and only if $[J,K]=0$. Now expand

[TABLE]

which implies the advertised claim.

Appendix D Orthogonal Transformation of Tracefree Symmetric Matrices into a Matrix with Zeros on the Diagonal

Given a symmetric matrix $K\in\mathbb{R}_{sym}^{d\times d}$ with $\operatorname{Tr}K=0$, we seek to find an orthogonal matrix $U\in O(\mathbb{R}^{d})$ such that $UKU^{T}$ has zeros on the diagonal. This is a crucial step in Algorithms 1 and 2 and has been addressed in various places in the literature (see for instance Kaz (88) or Bha (97), Chapter 2, Section 2). For the convenience of the reader, in the following we summarize an algorithm very similar to the one in Kaz (88).

Since $K$ is symmetric, there exists an orthogonal matrix $U_{0}\in O(\mathbb{R}^{d})$ such that $U_{0}KU_{0}^{T}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{d})$. Now the algorithm proceeds iteratively, orthogonally transforming this matrix into one with the first diagonal entry vanishing, then the first two diagonal entries vanishing, etc, until after $d$ steps we are left with a matrix with zeros on the diagonal. Starting with $\lambda_{1}$, assume that $\lambda_{1}\neq 0$ (otherwise proceed with $\lambda_{2}$). Since $\operatorname{Tr}(K)=\operatorname{Tr}(U_{0}KU_{0}^{T})=\sum\lambda_{i}=0$, there exists $\lambda_{j}$, $j\in\{2,\ldots,d\}$ such that $\lambda_{1}\lambda_{j}<0$ (i.e. $\lambda_{1}$ and $\lambda_{j}$ have opposing signs). We now apply a rotation in the $1j$-plane to transform the first diagonal entry into zero. More specifically, let

[TABLE]

with $\alpha=\arctan\sqrt{-\frac{\lambda_{1}}{\lambda_{j}}}.$ We then have $(U_{1}U_{0}KU_{0}^{T}U_{1}^{T})_{11}=0$. Now the same procedure can be applied to the second diagonal entry \text{\lambda}_{2}, leading to the matrix $U_{2}U_{1}U_{0}KU_{0}^{T}U_{1}^{T}U_{2}^{T}$ with

[TABLE]

Iterating this process, we obtain that $U_{d}\ldots U_{1}U_{0}KU_{0}^{T}U_{1}^{T}\ldots U_{d}^{T}$ has zeros on the diagonal, so $U_{d}\ldots U_{1}U_{0}\in O(\mathbb{R}^{d})$ is the required orthogonal transformation.