Variance Reduction Using Nonreversible Langevin Samplers

A standard approach to computing expectations with respect to a given target measure is to introduce an overdamped Langevin equation which is reversible with respect to the target distribution, and to approximate the expectation by a time-averaging estimator. As has been noted in recent papers [30, 37, 61, 72], introducing an appropriately chosen nonreversible component to the dynamics is beneficial, both in terms of reducing the asymptotic variance and of speeding up convergence to the target distribution. In this paper we present a detailed study of the dependence of the asymptotic variance on the deviation from reversibility. Our theoretical findings are supported by numerical simulations.


Motivation
In various applications arising in statistical mechanics, biochemistry, data science and machine learning [19,38,39,42], it is often necessary to compute expectations of an observable f with respect to a target probability distribution π(dx) on R d with density π(x) with respect to the Lebesgue measure, known up to the normalization constant. 1 When the dimension d is large, standard deterministic quadrature approaches become intractable, and one typically resorts to Markov-chain Monte Carlo (MCMC) methods [19,39,63]. In this approach, π( f ) is approximated by a long-time average of the form: where X t is a Markov process chosen to be ergodic with respect to the target distribution π. The Birkhoff-von Neumann Ergodic theorem [31,34,60] states that, for every observable f ∈ L 1 (π) we have If π possesses a smooth, strictly positive density, then a natural choice for X t is the overdamped Langevin dynamics where W t is a standard Brownian motion on R d . Assuming that (3) possesses a unique strong solution which is non-explosive, the process X t is ergodic, with unique invariant distribution π, such that (2) holds. Under additional assumptions on the distribution π and on the observable, this convergence result is accompanied by a central limit theorem which characterizes the asymptotic distribution of the fluctuations of π where σ 2 f is known as the asymptotic variance for the observable f . For the reversible process (3) started from stationarity (i.e. X 0 ∼ π), the Kipnis-Varadhan theorem [12,32] implies that (4) holds with asymptotic variance where ·, · π is the inner product in L 2 (π) and S is the infinitesimal generator of X t defined by S = ∇ log π(x)∇ · + .
Sampling methods based on Langevin diffusions have become increasingly popular due to their wide applicability and relative ease of implementation. In practice, a discretisation of (3) may be used, which in general, will not be ergodic with respect to the target distribution π. Thus, the discretisation is augmented with a Metropolis-Hastings accept/reject step [65] which guarantees that the resulting Markov chain remains reversible with respect to π. The resulting scheme is known as the Metropolis-Adjusted Langevin Algorithm (MALA), see [65]. The computational cost required to approximate π( f ) using π T ( f ) to a given level of accuracy depends on the target distribution π, the observable f and the process X t over which the long-time average is generated. Quite often, X t will exhibit some form of metastability [10,11,36,66]: the process X t will remain trapped for a long time exploring one mode, with transitions between different modes occurring over longer timescales. When the observable Fig. 1 Trajectories of a reversible MALA chain and a nonreversible Langevin sampler generating samples from a warped Gaussian distribution. The mass of the distribution is concentrated along a one-dimensional submanifold. Reversible samplers, such as MALA, perform a very slow exploration of the distribution, spending large amounts of time concentrated around a small region of the submanifold. On the other hand, nonreversible samplers are able to perform long "jumps" along the level-curves of the distribution, thus able to explore the distribution far more effectively (Color figure online) depends directly on the metastable degrees of freedom (i.e. the observable takes different values in different metastable regions), the asymptotic variance σ 2 f of the estimator π T ( f ) may be very large. As a result, more samples are required to obtain an estimate of π( f ) with the desired accuracy. A similar scenario arises when the mass of π is tightly concentrated along a low-dimensional submanifold of R d , as illustrated in Fig. 1. In this case, reversible dynamics, such as (3) will cause a very slow exploration of the support of π. As a result, π T ( f ) will exhibit very large variance for observables which vary strongly along the manifold.
As there are infinitely 2 many Markov processes with invariant distribution π, a natural question is whether such a process can be chosen to have optimal performance. The two standard optimality criteria that are commonly used are: (a) With respect to speeding up convergence to the target distribution. (b) With respect to minimizing the asymptotic variance.
These criteria can be used in order to introduce a partial ordering in the set of Markov chains or diffusion processes that are ergodic with respect to a given probability distribution [50,59]. From a practical perspective, the definite optimality criterion is that of minimizing the computational cost. We address this issue (at least partially) later in this paper.
Within the family of reversible samplers, much work has been done to derive samplers which exhibit good (if not optimal) computational performance. This has motivated a number of variants of MALA which exploit geometric features of π to explore the state space more effectively, including preconditioned MALA [64], Riemannian Manifold MALA [20] and Stochastic Newton Methods [45].

Nonreversible Langevin Dynamics
An MCMC scheme which departs from the assumption of reversible dynamics is Hamiltonian MCMC [53], which has proved successful in Bayesian inference. By augmenting the state space with a momentum variable, proposals for the next step of the Markov chain are generated by following Hamiltonian dynamics over a large, fixed time interval. The resulting nonreversible chain is able to make distant proposals. Various methods have been proposed which are related to this general idea of breaking nonreversibility by introducing an additional dimension to the state space and introducing dynamics which explore the enlarged space while still preserving the equilibrium distribution. In particular, the lifting method [15,27,71] is one such method applied to discrete state systems, where the Markov chain is "lifted" from the state space onto the space ×{1, −1}. The transition probabilities in each copy are modified to introduce transitions between the copies to preserve the invariant distribution but now promote the sampler to generate long strides or trajectories. Similar methods based on introducing nonreversibility into the dynamics of discrete state chains to speed up convergence have been applied with success in various applications [9,26,52,68,69]. These methods are also reminiscent of parallel tempering or replica exchange MCMC [51], which are aimed at efficiently sampling from multimodal distributions.
It is well documented that breaking detailed balance, i.e. considering a nonreversible diffusion process that is ergodic with respect to π, can help accelerate convergence to equilibrium. In fact, it has been proved that, among all diffusion processes with additive noise that are ergodic with respect to π, the reversible dynamics (3) has the slowest rate of convergence to equilibrium, measured in terms of the spectral gap of the generator in L 2 (R d ; π) =: L 2 (π), c.f. [37]. Adding a drift to (3) that is divergence-free with respect to π and that preserves the invariant measure of the dynamics will always accelerate convergence to equilibrium [28,29,37,54,61,72]. The optimal nonreversible perturbation can be identified and obtained in an algorithmic manner for diffusions with linear drift, whose invariant distribution is Gaussian [37]; see also [72].
The effect of nonreversibility in the dynamics on the asymptotic variance has also been studied. In [61] it is shown that small antisymmetric perturbations of the reversible dynamics always decrease the asymptotic variance, and more recently in [62] Friedlin-Wentzell theory is used to study the limit of infinitely strong antisymmetric perturbations. In [30] the authors use spectral theory for selfadjoint operators to study the effect of antisymmetric perturbations on diffusions on both R d and on compact manifolds and provide a general comparison result between reversible and nonreversible diffusions. This work is related with previous studies on the behavior of the spectral gap of the generator when the strength of the nonreversible perturbation is increased [13,18].

Objectives of This Paper
In this paper we present a detailed analytical and computational study of the effect of adding a nonreversible drift to the reversible overdamped Langevin dynamics (3) on the asymptotic variance σ 2 f . We will consider the nonreversible dynamics defined by [28,29] : where the smooth vector field γ is taken to be divergence-free with respect to the distribution π, ∇ · γ π = 0.
There are infinitely many vector fields that satisfy (7) and a general formula can be derived using Poincaré's lemma [29,62]. We can, for example, construct such vector fields by taking In (6) we have already normalized the various vector fields and we have introduced a parameter α which measures the strength of the deviation from reversible dynamics. The generator of the diffusion process X γ t can be decomposed into a symmetric and an antisymmetric part in L 2 (π), representing the reversible and irreversible parts of the dynamics, respectively [56,Sect. 4.6]: where S is given in (5) and A = γ · ∇. We are particularly interested in quantifying the reduction in the asymptotic variance obtained caused by breaking detailed balance. It is well known [32,33] that for square integrable observables, the asymptotic variance, denoted by σ 2 f (α) can be written in terms of the solution of the Poisson equation as Our first goal is to obtain an explicit formula for σ 2 f (α) in terms of S and A. Moreover, through a spectral decomposition of the operator (−S) −1 A we investigate the dependence of the asymptotic variance on the strength of the nonreversible perturbation. Based on earlier work by Golden and Papanicolaou [22], Majda et al. [3,44] and Bhattacharya et al. [7,8], in [58], an expression is obtained for the effective self-diffusion coefficient for a tagged particle whose microscopic behaviour is determined by a nonreversible Markov process. The basic idea is the introduction of the operator which is skewadjoint in the weighted Sobolev space where ·, · π denotes the L 2 (π) inner product; see Sect. 3. We follow a similar approach in this paper, from which a quite explicit formula for the asymptotic variance σ 2 f (α) is derived using the spectral theorem for selfadjoint operators. From this expression we can immediately deduce that the addition of a nonreversible perturbation reduces the asymptotic variance and we can also study the small and large α asymptotics of σ 2 f ; see Theorem 4. One of the results that we prove is that the large α behaviour of the asymptotic variance depends on the detailed properties of the vector field γ : when the nullspace of the antisymmetric part of the generator A = γ · ∇ consists of only constants in H 1 , then the asymptotic variance converges to 0. Indeed, in this case, in the |α| → ∞ limit, the limiting dynamics become deterministic, characterized by the Liouville operator γ · ∇. On the other hand, when the nullspace of A contains nontrivial functions in H 1 there will exist observables for which the asymptotic variance σ 2 f converges to a positive constant as |α| → ∞. The effect of the antisymmetric part on the long time dynamics of diffusion processes has been also studied extensively in the context of turbulent diffusion [43] and fluid mixing. The effect of an incompressible flow on the convergence of the solution to the advection-diffusion equation on a compact manifold to its mean value (i.e. when π ≡ 1) was first studied in [13].
In particular, the concept of a relaxation enhancing flow was introduced and it was shown that a divergence-free flow is relaxation enhancing if and only if the Liouville operator A = v · ∇ has no eigenfunctions in the Sobolev space H 1 . An equivalent formulation of this result is that an incompressible flow is relaxation enhancing if and only if it is weakly mixing. Examples of relaxation enhancing flows are given in [13]. This problem was studied further in [18], where it also mentioned that there are very few examples of relaxation enhancing flows. In [18] it is shown that the spectral gap of the advection-diffusion operator L = αv · ∇ + , for a divergence-free drift v and α ∈ R, remains bounded above in the limit as α → ±∞ by a negative constant if and only if the advection operator has an eigenfunction in H 1 . These results are reminiscent of the results we mention above on the necessary and sufficient condition to obtain a reduction of asymptotic variance in the limit α → ±∞.
Our analysis of the asymptotic variance σ 2 f , based on the careful study of the Poisson Eq. (10), enables us to study in detail the problem of finding the nonreversible perturbation giving rise to minimum asymptotic variance for diffusions with linear drift, i.e. diffusions whose invariant measure is Gaussian, over a large class of observables. Diffusions with linear drift were considered in [37] where the optimal nonreversible perturbation with respect to accelerating convergence was obtained. For linear and quadratic observables, we can give a complete solution to this problem, and construct nonreversible perturbations that provide a dramatic reduction in asymptotic variance. Moreover, we demonstrate that the conditions under which the variance is reduced are very different from those of maximising the spectral gap discussed in [37]. In particular, we show how a nonreversible perturbation can dramatically reduce the asymptotic variance for the estimator π T ( f ), even though no such improvement can be made on the rate of convergence to equilibrium.
Guided by our theoretical results, we can then study numerically the reduction in the asymptotic variance due to the addition of a nonreversible drift for some toy models from molecular dynamics. In particular, we study the problem of computing expectations of observables with respect to a warped Gaussian [23] in two dimensions, as well as a simple model for a dimer in a solvent [38]. The numerical experiments reported in this paper illustrate that a judicious choice of the nonreversible perturbation, dependent on the target distribution and the observable, can dramatically reduce the asymptotic variance.
To compute π T ( f ) numerically, we use an Euler-Maruyama discretisation of (6). The resulting discretisation error introduces an additional bias in the estimator for π( f ), see [47] for a comprehensive error analysis. This imposes additional constraints on the magnitude of the nonreversible drift, since increasing α arbitrarily will give rise to a large discretisation error which must be controlled by taking smaller timesteps. A natural question is whether the increase in the computational cost due to the necessity of taking smaller timesteps negates any benefits of the resulting variance reduction. To study this problem, we compare the computational cost of the unadjusted nonreversible Langevin sampler with the corresponding MALA scheme. 3 Our numerical results, together with the theoretical analysis for diffusions with linear drift, show that the nonreversible Langevin sampler can outperform the MALA algorithm, provided that the nonreversible perturbation is well-chosen. Finally, we consider a higher order numerical scheme for generating samples of (6), based on splitting the reversible and nonreversible dynamics. Numerically, we investigate the properties of this integrator, and demonstrate that its improved stability and discretisation error make it a good numerical scheme for computing the estimator π T ( f ) using a nonreversible diffusion.
The rest of the paper is organized as follows. In Sect. 2 we describe how the central limit theorem (4) arises from the solution of the Poisson equation associated with the generator of the dynamics. In Sect. 3 we analyse the asymptotic variance and formulate the problem of finding the optimal perturbation with respect to minimising σ 2 f , for a fixed observable, or over the space of square-integrable observables. Moreover, following the analysis described in [58], we derive a spectral representation for the asymptotic variance σ 2 f , in terms of the discrete spectrum of the operator G defined in (12), and recover estimates for the asymptotic variance for any value of α. In Sect. 4, we consider the case of Gaussian diffusions, which are amenable to explicit calculation to demonstrate the theory presented in this paper. In Sect. 5, we provide various numerical examples to complement the theoretical results. Finally, in Sect. 6 we describe the bias-variance tradeoff for nonreversible Langevin samplers, and explore their computational cost. Conclusions and discussion on further work are presented in Sect. 7.

The Central Limit Theorem and Estimates on the Asymptotic Variance via the Poisson Equation
In this section we make explicit some sufficient conditions under which the estimator where X γ t denotes the solution of (6), satisfies a central limit theorem of the form (4). The fundamental ingredient for proving such a central limit theorem is establishing the wellposedness of the Poisson equation for all bounded and continous functions f : R d → R, where L is defined by (9), and obtaining estimates on the growth of the unique solution φ. As described previously, we shall assume that π possesses a smooth, strictly positive density, also denoted by π(x), such that R d π(x) dx < ∞ and that the SDE (6) has a unique strong solution for all α ∈ R. Applying the results detailed in [21,49], we shall assume that the process X γ t possesses a Lyapunov function, which is sufficient to ensure the exponential ergodicity of X γ t , [46,48], as detailed in the following assumption. Assumption A (Foster-Lyapunov Criterion) There exists a function U : R d → R and constants c > 0 and b ∈ R such that π(U ) < ∞ and where 1 C is the indicator function over a petite set. For the definition of a petite set we refer the reader to [48]. For the generator L corresponding to the process (6), compact sets are always petite. As noted in [65], a sufficient condition on π for (6) to possess a Lyapunov function is the following. Assumption B The density π is bounded and for some 0 < β < 1: Moreover, if the nonreversible term is of the form γ (x) = J ∇ log π(x), for J ∈ R d×d antisymmetric, then this choice of U is a Lyapunov function for X γ t defined by (6) for all α ∈ R.
Proof For this choice of γ the generator of (6) has the form Thus, by assumption (15), there exists > 0 and M such that for |x| > M: and so Finally, we note that since π is bounded, then U is bounded above from zero uniformly. Thus, U can be rescaled to satisfy the condition U ≥ 1, as is required by the Foster-Lyapunov criterion.

Remark 1 Note that Assumption B holds trivially when π is a Gaussian distribution in
so that (14) holds. On the other hand, consider the case when when π is a Student's t- Then it is straightforward to check that (14) will not hold. Indeed, since |∇ log π(x)| → 0, by [65, Theorem 2.4] this process is not exponentially ergodic.
If condition (14) holds for X γ t , then the process will be exponentially ergodic. More specifically, the law of the process X γ t started from a point x ∈ R d will converge exponentially fast in the total variation norm to the equilibrium distribution π. In particular, denoting by (P γ t ) t≥0 the semigroup associated with the diffusion process (6), we have the following result. where p γ t (x, ·) denotes the law of X γ t given X γ 0 = x and ||·|| T V denotes the total variation norm. In particular, for any probability measure ν such that U ∈ L 1 (ν), For a central limit theorem to hold for the process X γ t , and thus for σ 2 f to be finite, it is necessary that the Poisson Eq. (10) is well-posed. The Foster-Lyapunov condition (14) is sufficient for this to hold.
The technique of using a Poisson equation to obtain a central limit theorem for an additive functional of a Markov process is widely known, [5,6,55]. The approach is based on the fact that, at least formally, we can decompose π t ( f ) − π( f ) into a martingale and a "remainder" term: Considering the rescaling √ t (π t ( f ) − π( f )), the martingale term √ t M t will converge in distribution to a Gaussian random variable with mean 0 and variance by the central limit theorem for martingales [25]. It remains to control the remainder term √ t R t . We distinguish between two cases: If X γ 0 ∼ π, then since φ ∈ L 2 (π), √ t R t converges to 0 in L 2 (π) and the result follows. In the more general case we must resort to a "propagation of chaos" argument (c.f. [12, Sect. 8]), i.e. apply Theorem 1 to show that as r → ∞, for all continuous bounded functions H . The result then follows by decomposing and applying the propagation of chaos argument to the second term. The conclusion is summarized in the following result, which provides a central limit theorem for X γ t starting from an arbitrary initial distribution ν.

Theorem 3 [21, Thm 4.4] If Assumption A holds for Lyapunov function U , then for any f such that f
In the remainder of this paper we shall study the dependence of φ, and thus σ 2 f on the choice of non-reversible perturbation γ . We note that (16) is precisely the Dirichlet form associated with the dynamics L evaluated at the solution φ of the Poisson Eq. (13).

Mathematical Setting
We present a few definitions from [33] that will be useful in the sequel. Let L 2 0 (π) denote the set of L 2 (π)-integrable functions with zero mean, with corresponding inner product ·, · π and norm ||·|| π . Consider the operator S given by (5) densely defined on L 2 0 (π). For k ∈ N, given the family of seminorms we define the function spaces It follows that ||·|| k is a norm on H k . For k ≥ 1, the norm · k satisfies the parallelogram identity and, consequently, the completion of H k with respect to the norm · k , which is denoted by H k , is a Hilbert space. The inner product ·, · k in H k is defined through polarization. It is easy to check that, for f, g ∈ D(L), We note that H 0 = L 2 0 (π). A careful analysis of the function space H k is presented in [33]. The operator S is symmetric with respect to π and can be extended to a selfadjoint operator on L 2 0 (π), which is also denoted by S, with domain D(S) = H 2 . We shall make the following assumption on π which is required, in addition to Assumption B, to ensure that L possesses a spectral gap in L 2 0 (π).
In the following lemma we establish a number of fundamental properties relating to the spectrum of L. Using these results, we then establish the well-posedness of the Poisson Eq. (13) for any f ∈ L 2 0 (π).
Lemma 2 Suppose that Assumptions B and C hold. Then the embedding H 1 ⊂ L 2 0 (π) is compact. For all α ∈ R, the operator L satisfies the following Poincaré inequality in L 2 0 (π): where λ is a positive constant, independent of the nonreversible perturbation. Moreover, for for some δ ∈ (0, 1) and M 1 > 0. Applying [40,Theorem 8.5.3], relations (19) and (21) imply the compactness of the embedding H 1 ⊂ L 2 0 (π). As a result, following [40, Theorem 8.6.1], this is sufficient for the Poincaré inequality to hold for S, i.e. λ ||g|| 2 π ≤ ||g|| 2 1 = g, (−S)g π , g ∈ H 1 . from which (20) follows by the antisymmetry of A in L 2 0 (π). The existence of this Poincaré inequality implies that the semigroup P γ t associated with (6) converges exponentially fast to equilibrium, that is, for all f ∈ L 2 0 (π): By (22) it follows that φ ∈ L 2 0 (π) and from (23) we have −Lφ = f . Moreover, we have the bound Throughout this section, we shall assume that Assumptions B and C hold, and moreover, for simplicity we shall make the following additional assumption: Assumption D The nonreversible perturbation γ is smooth and bounded in L ∞ . We believe that this assumption could be relaxed, see in particular [30] for more general assumptions under which the following results should hold. We stick here to a simple presentation. In practice, this assumption is not very stringent. Indeed, suppose that there exists a smooth function ψ : then γ satisfies (7) since using the fact that J is antisymmetric. Moreover, it is clear that γ (x) is bounded and smooth, thus satisfying Assumption D. In particular, if V has compact level sets (which holds for example if V is a nonnegative polynomial function), then condition (24) is satisfied by chosing ψ to be a smooth non-negative function with compact support.
We note that γ of the form (25) leaves the potential V invariant under the flowż t = γ (z t ), i.e. V (z t ) is constant for all t ≥ 0. Thus, for large |α|, the flow αγ will result in rapid exploration of the level surfaces of V , but the motion of X t between level surfaces is entirely due to the reversible dynamics of the process. In particular, for potentials with energy barriers, the transition time for X γ t to cross a barrier will still satisfy the same Arrenhius law as the corresponding reversible process. Other choices of flow γ are possible. For example, one could alternatively consider a skew-symmetric matrix function J (x) as detailed in [41]. The corresponding flow would then be defined by It is straightforward to check that ∇ ·(γ π) = 0, using the fact that J (x) is skew-symmetric. If additionally, the matrix function J is smooth with bounded derivative and compact support, then γ satisfies Assumption D. As detailed in [41] one can further generalise this choice of dynamics by additionally introducing a space dependent diffusion tensor, and an appropriate correction of the drift to maintain ergodicity with respect to π. We do not consider this choice of dynamics in this paper, noting that most of the presented results can be readily generalized to this scenario. Given that Assumption D holds for all function u ∈ H 1 , then we have (7), so that the operator A :

An Expression for the Asymptotic Variance
The main objective of this paper is to study the effect of the nonreversible perturbation γ on the asymptotic variance σ 2 f , with the aim of choosing γ so that σ 2 f is minimized. Integrating (16) by parts we can express σ 2 f in terms of the Dirichlet form for L as follows: where the last line follows from the fact that A is antisymmetric in L 2 0 (π). Starting from (27) we can obtain a quite explicit characterisation of σ 2 f . Indeed, given f ∈ L 2 0 (π), we can rewrite the last line of (27) as where [·] S denotes the symmetric part of the operator. Using the fact that −L is invertible on L 2 0 (π) for all α ∈ R by Lemma 2, the following result yields an expression for σ 2 f (α) in terms of S, A and α.
where A * = −A denotes the adjoint of A in L 2 0 (π). In particular, for all f ∈ L 2 0 (π), Proof Since −S is positive, the operator Q = (−S) − 1 2 : L 2 0 (π) → H 1 can be defined using functional calculus. Consider the operator C := QA. We can write −S + α 2 A * (−S) −1 A = −S + α 2 C * C, which can be shown to be closed in L 2 0 (π) with domain H 2 . Moreover, this operator has nullspace {0}, so that the inverse is also densely defined on L 2 0 (π). To show that (29) holds, we expand the left hand side to get as required. The inequality (30) is then simply a consequence of the fact that where denotes the partial ordering between selfadjoint operators on L 2 0 (π). Thus, the asymptotic variance is never increased by introducing a nonreversible perturbation, for all f ∈ L 2 0 (π). This had already been noted in [61] where an expression for the the asymptotic was derived as the curvature of the rate function of the empirical measure, and also in [30] using an approach similar to that above. Expression (28) provides us with a formula for σ 2 f in terms of a symmetric quadratic form which is explicit in terms of A and S.

Quantitative Estimates for the Asymptotic Variance
In this section we derive quantitative versions of (29) and (30), using techniques developed in [58] for the analysis of the Green-Kubo formula, which is itself based on earlier work on the estimation of the eddy diffusivity in turbulent diffusion [2,8,44,57].
Following the approach of [22,58], to quantify the effect of the antisymmetric perturbation αA on the asymptotic variance, we define the operator G = (−S) −1 A : H 1 → H 1 . First we note that we can rewrite (29) as follows: and therefore, from (28) and (29): From the boundedness of the nonreversible perturbation we have the following properties:

Lemma 4 Suppose that Assumption D holds, then the operator
so that G is antisymmetric. Since G is also bounded, it follows that G is skewadjoint on H 1 .
Moreover, under appropriate assumptions on the target distribution π, one can show that the operator G is compact.

Lemma 5 Suppose that Assumptions B, C and D hold, then the operator G is compact on
Proof By Lemma 2, the embedding H 1 ⊂ L 2 0 (π) is compact, and it follows immediately that H 2 ⊂ H 1 is a compact embedding. Moreover, since γ is bounded we have, wheref N ∈ N andf n = f , e n 1 . Consequently, using this spectral decomposition, we can rewrite (32) as The conclusion of the above computation is summarized in the following result.
Theorem 4 Suppose that Assumptions B, C and D hold. Let f ∈ L 2 0 (π) and α ∈ R, then the asymptotic variance σ 2 f (α) corresponding to X γ t is given by (34). In particular, we obtain the following limiting values for the asymptotic variance: and, lim From (34) we obtain the following bounds on the asymptotic variance The problem of choosing the optimal nonreversible perturbation in (6) to minimize the asymptotic variance over all observables in L 2 0 (π) can be expressed as the following min-max problem: where, for some constant M > 0, A M denotes the set of admissible nonreversible drifts, typically, This is equivalent to finding γ ∈ A M which solves this max-min problem where σ [M] denotes the spectrum of the operator M on L 2 0 (π). It is important to make the distinction between this problem and that considered in [37] for finding the optimal spectral gap, namely finding γ ∈ A M such that From (37) and (38) it is evident that the decrease in asymptotic variance, and the increase in spectral gap arising from a nonreversible perturbation are due to very different mechanisms. For the latter problem we see that the increase of spectral gap arises from the nonnormality of S +A. In addition, since the operator A * (−S) −1 A is nonnegative in H 1 , a nonreversible perturbation cannot increase the variance. This is in contrast with the problem of maximising the speed of convergence to equilibrium, as was considered in [37], where increasing the strength of the nonreversible perturbation could result in a decrease of the speed of convergence.
We note however that a nonreversible perturbation γ which solves the min-max problem (37) need not be a good candidate for reducing the variance of π T ( f ) for a given fixed observable f ∈ L 2 0 (π). Indeed, unless N = {0} for all g ∈ H 1 , there will always be an observable for which the nonreversible perturbation does not reduce the variance, since if g ∈ N is nonzero, f = (−S)g is nonzero and such that σ 2 f (α) = σ 2 f (0) for any α. Thus, from a practical point of view, it makes more sense to consider the problem choosing γ to minimise the asymptotic variance of π T ( f ) for a particular observable. To this end, for a fixed f ∈ L 2 0 (π), we can identify two distinct cases: More generally, consider f ∈ L 2 (π) so that f −π( f ) ∈ L 2 0 (π). Assuming we can increase α arbitrarily, and neglecting any computational issues arising from the resulting discretisation error (which will be discussed in Sect. 6), the problem of minimising the asymptotic variance of f reduces to finding γ such that holds in H 1 . Checking this condition requires the solution of an elliptic boundary value problem, thus this condition is not of practical use. Nonetheless, we can derive some intuition from (39). Clearly, if we can choose γ so that N = {0} in H 1 , the nonreversible perturbation will be optimal for all observables f ∈ L 2 (π). In general, it might not be possible to find γ that satisfies this condition. For the nonreversible drift considered in [37], namely γ = J ∇V , with J = −J , N will always be nontrivial; indeed, in this case H • V ∈ N for all functions H such that H • V ∈ H 1 . In this case, it is always possible to choose an observable f such that γ will not be optimal for f , in the sense that σ 2 f (α) is nonzero in the limit α → ∞. We also remark that these asymptotic results, in particular the distinction between these two cases is reminiscent of similar results that have been obtained in the context of turbulent diffusion [44,57].
An analogous classification of the asymptotic behaviour of the spectral gap of the operator (9) in the limit of large α is considered in [18] on a compact manifold. Indeed, in [18,Theorem 1] it is determined that the spectral gap is finite in the limit of α → ±∞ if and only if A has a nonconstant eigenfunction in H 1 . However, one should note that, the asymptotic variance σ 2 f (α) may converge to 0 as α → ±∞ even when the spectral gap is finite, see also [61, Example 2.9] for a counterexample.

A Two Dimensional Example
In this section we present a simple example on R 2 . Consider the problem of calculating π( f ) with respect to the Gaussian distribution π(x) = 1 Z e −|x| 2 /2 . The reversible overdamped Langevin equation corresponding to π is given by the OU process: We introduce an antisymmetric perturbation αγ (X t ) where γ (x) is the flow given by The infinitesimal generator for the perturbed nonreversible dynamics is In polar coordinates, this is given by 2π). As |α| → +∞, we expect the deterministic flow to move increasingly along the level curves of the distribution π(x). Given an observable f , any variance in π T ( f ) arising from the variation of f along the level curves should vanish as |α| → ∞, leaving only the variance contributed by the variation of f between level curves. We make this precise with a particular example. Consider the observable f (x 1 , x 2 ) = 2x 2 1 , expressible in polar coordinates as f (r, θ) = 2r 2 cos 2 (θ ) = r 2 (1 + cos(2θ)) .
Noting that π( f ) = 2, it is straightforward to check that the Poisson equation −Lφ = f − π( f ), has mean zero solution given by The asymptotic variance can be evaluated directly as follows (cos(2θ) − α sin(2θ)) 1 + α 2 r 2 (1 + cos(2θ)) − 2 re −r 2 /2 dr dθ Therefore, as |α| → ∞, the asymptotic variance converges to 4. We note that, in this case, where r 2 /2 − 1 is perpendicular to r 2 /2 cos(2θ) in H 1 . Since the nullspace of A consists of all H 1 functions which depend only on r , we havef N = r 2 /2 − 1. It follows that which agrees with the conclusions of Theorem 4. As |α| → ∞, the motion in the θ direction is averaged out. Indeed, in this limit, 2 f N 2 1 corresponds to the asymptotic variance of the observable T 0 r 2 t dt, where r t is the following 1D reversible process, For more general observables f (r, θ), we expect maximum variance reduction as |α| → ∞ when f varies strongly with respect to θ . In the other extreme, we expect zero improvement when f depends only on r . This intuition is formalised in the following section, in particular in Proposition 1 and the subsequent bound (51). For more general potentials, using a flow field of the form γ (x) = J ∇ log π(x), J = −J , the mechanism for reducing the asymptotic variance is analogous: the large antisymmetric drift gives rise to fast deterministic mixing along the level curves of the potential, while the reversible dynamics induce slow diffusive motion along the gradient of the potential. When the nullspace of A is trivial in H 1 , the fast deterministic flow is ergodic, so that, for α large, the antisymmetric component will cause a rapid exploration of the entire state space. Consequently, the asymptotic variance converges to 0 as α → ∞. On the other hand, if A has a nontrivial nullspace, the antisymmetric perturbation is no longer ergodic, and the state space can be decomposed into components such that the rapid flow behaves ergodically in each individual component. In the limit of large α, X t becomes a fast-slow system, with rapid exploration within the ergodic components coupled to a slow diffusion between components. Very recently, Rey-Bellet and Spiliopoulos [62] have applied Freidlin-Wenzell theory to rigorously analyse this case in the large α limit for a large class of potentials.

Nonreversible Perturbations of Gaussian Diffusions
For the case when the target distribution is Gaussian, the SDE (6) for X γ t is linear. In this case, we can obtain an explicit analytical expression for the asymptotic variance for a large class of observables f . Indeed, consider the nonsymmetric Ornstein-Uhlenbeck process in where J is an antisymmetric matrix, α > 0, and W t is a standard d-dimensional Brownian motion. The stationary distribution π(x) is N (0, I ), independent of α and J . Although this system does not fall under the framework of Theorem 4 we are still able to obtain analogous conditions for a reduction in the asymptotic variance. The objective of this section is to mirror the results for speeding up convergence to equilibrium of X t that were derived in [37] to the case of minimizing the asymptotic variance. In particular, following arguments similar to [58,Sect. 4.2], an explicit formula for the asymptotic variance will be derived, from which an optimal J can be chosen, in a manner similar to [37]. We note that for the process (41), the optimal nonreversible perturbation obtained in [37] does not provide any increase to the rate of convergence to equilibrium, since all eigenvalues of the covariance matrix of the Gaussian stationary distribution are the same. Nonetheless, in this section we show that for certain observables, the asymptotic variance of π T ( f ) can be dramatically decreased.

Explicit Formula for the Asymptotic Variance
We shall assume that the observable f is a quadratic functional of the form where M ∈ R d×d is a symmetric positive definite matrix, l ∈ R d and k is a constant, chosen so that f (x) is centered with respect to π(x): k = − Tr M. Consider the Poisson Eq. (10): where L is the infinitesimal generator of (41) given by L = −(I + α J )x · ∇ + . For such observables we can solve the Poisson Eq. (43) analytically and obtain a closed-form formula for the asymptotic variance for the observable f .

Proposition 1 Let A = (I + α J ) = (I − α J ). The unique mean zero solution of the Poisson Eq. (43) is given by
and Moreover, the asymptotic variance is given by

the Frobenius norm of A, and [A, B] = AB − B A is the commutator of A and B. In particular,
Proof Clearly φ(x) must also be a quadratic function of x, so we make the ansatz Plugging this into (43) we obtain Comparing equal powers of x we have for all x ∈ R d . Since x is arbitrary, it follows that Equation (46a) is a Lyapunov equation, which is well-posed as M is positive definite and spec(A) ⊂ {λ ∈ C : Reλ > 0}. Indeed, for C given by both (46a) and (46c) are satisfied. The asymptotic variance σ 2 f (α) is then given by (see (27)), Using the fact that C and M are symmetric, Therefore, Since A = I − α J , where J is antisymmetric, we have that, for all l ∈ R d : Therefore, we obtain the following expression for σ 2 f (α): Writing the first term as follows: Since J is skew-symmetric, the matrix exponential e α J s is a rotation matrix. Thus, the matrix e α J s Me −α J s has the same eigenvalues as M and M . From [4, III. 6.14] we have where U is a d × d orthogonal matrix with the first N columns spanning N . Then where l N is the projection of l onto N . Thus, for quadratic observables, from (49), the asymptotic variance has the following lower bound in the limit of large α: In particular, for diffusions with linear drift, the asymptotic variance cannot be decreased arbitrarily. The worst case scenario is when f (x) = m |x| 2 − d , for some m > 0, for which the asymptotic variance will not be decreased, for any antisymmetric matrix J and α ∈ R, analogous to the situation which occurs in [37] for maximising the spectral gap, when all the eigenvalues of the covariance matrix are equal.

Finding the Optimal Perturbation
We first focus on the case when the observable is linear, i.e. f (x) = l · x, similar to that considered in [58,Sect. 4.2], so that σ 2 f = 2l · (I + α 2 J J ) −1 l. Suppose that we fix α ∈ R, and wish to choose J which minimizes the asymptotic variance subject to ||J || F = 1. In this case, we have the equality which is optimal if l is orthogonal to N . Thus, the best we can do is to choose J such that l is an eigenvector of J J with maximal eigenvalue. This can be done by choosing a unit vector ω ∈ R d orthogonal to l and setting wherel = l/|l| Then, is a projector onto {l, ω} and ||J || F = 1, where ||·|| F is the Frobenius norm with respect to the Euclidean basis. In this case We note in addition that, when minimising the asymptotic variance, there are many antisymmetric matrices J which give the minimal asymptotic variance. As an example, let d = 3, consider the observables f (x) = l (i) · x for l (1) = (0, 1, 1) / √ 2 and l (2) = (1, 0, 1) / √ 2 and l (3) = (1, −1, 1) / √ 3, respectively. We choose J to be In this case The decay of σ 2 f (α) as α → ∞ is strongly dependent on the nullspace of J , given by In Fig. 2a we plot the asymptotic variance for these observables. For f (x) = l (1) · x, see that as α → ∞ the asymptotic variance converges to 0. This is due to the fact that l (1) is perpendicular to the nullspace of J . Thus, for this observable the matrix J given by (52) is an optimal perturbation. For f (x) = l (2) · x, since l (2) is not orthogonal to ζ , as α → ∞, Finally, for f (x) = l (3) · x we observe that the asymptotic variance remains constant at 2, since l (3) ∈ N , the nonreversible perturbation has no effect on this observable. We now focus on the case when l = 0 so that f (x) = x · M x − Tr(M). In this case, it is not clear how to construct an optimal J , however the bound (51) suggests a good candidate for J . Suppose that M has eigenvalues λ 1 ≤ λ 2 ≤ · · · ≤ λ d with corresponding eigenvectors e 1 , . . . , e d . Suppose that d is even. Let be a partition of 1, . . . d into disjoint pairs. Define the antisymmetric matrix J by In this case e α J s can be decomposed into a product of rotations between the pairs of eigenvectors: where R v,w (θ ) is an anticlockwise rotation of angle θ in the {v, w} plane. Then For this choice of J , the asymptotic variance is given by Arguing by contradiction, this is minimized by choosing i k = k, and j k = (d − k + 1) for k = 1, . . . , d/2, in which case Clearly, however these bounds are only tight when M = m I, for m > 0.

Numerical Experiments
In this section we present three numerical examples illustrating the effects of the antisymmetric perturbation on the asymptotic variance σ 2 f . We will consider target measures defined by Gibbs distributions of the form: where V (x) is a given smooth, confining potential with finite unknown normalisation constant Z . We will also assume that the divergence-free vector field γ (x) is given by so that the drift of (6) will be given by The symmetric and antisymmetric parts of the generator in L 2 (T d ; π(x)dx) are, respectively, Provided that V ∈ H 1 (π) the nullspace of the antisymmetric part of the generator is always nonempty. Indeed, the antisymmetry of A implies that AV = 0. Thus, for any observable f such that π( f ) = 0 and V π the conclusion of Theorem 4 implies the asymptotic variance σ 2 f (α) will converge to a nonzero constant as α → ∞.

Periodic Distribution
In the first example we consider a two-dimensional target density given by (54) where the periodic potential V (x) is given by and β −1 = 0.1. Our objective is to calculate the expectation π( f ) of the observable In Fig. 3a we plot the asymptotic variance over α ∈ [0, 12]. The stepsize t for the Euler-Maruyama discretisation of (6) is chosen to be 10 −3 . For each α, M = 10 3 independent realisations of the Markov process are run, starting from X 0 = (0, 0), each for T = 10 5 time units to ensure that the process is close to stationarity. We observe that increasing α from 0 to 10 decreases the asymptotic variance by approximately two orders of magnitude. In Fig. 3b we plot the value of the estimator along with the confidence intervals estimated over 10 3 independent realisations. For this particular example, it appears that increasing the magnitude of the nonreversible perturbation does not appear to give rise to any noticeable increase in bias, while it significantly decreases the asymptotic variance, giving rise to a dramatic increase in performance of the estimator π T ( f ).

Warped Gaussian Distribution
As a second numerical example we consider computing an observable with respect to a two-dimensional warped Gaussian distribution [23], defined by (54) with (a) (b) Fig. 3 The asymptotic variance and value of the estimator π T ( f ), where T = 10 5 , for the target Gibbs distribution with potential defined by (56) and observable as in (57). The plot was generated from 10 3 independent realisations of an Euler-Maruyama discretisation of the nonreversible diffusion (6) The parameter b > 0 is chosen to be b = 0.05. The potental V (x) is plotted in Fig. 4. Our objective is to compute π( f ) where the observable f is given by: The nonreversible perturbation is chosen to be: In Fig. 5a we plot the asymptotic variance of the estimator π T ( f ), approximated from 10 3 independent realisations of the process, each run for 10 8 timesteps of size 10 −3 so that the total time is T = 10 5 . Each realisation was started at (0, 0). As expected, we observe a significant decrease in asymptotic variance as the magnitude of the nonreversible perturbation is increased. Increasing α from 0 to 10 gives rise to a decrease in variance by a factor of 80. In Fig. 5b we plot the value of the estimator π T averaged over the 10 3 realiations, with corresponding confidence intervals. As α is increased, the bias arising from the discretisation error also increases. To mitigate this bias one could decrease the stepsize, thus requiring more steps to generate π T ( f ) or otherwise introduce a Metropolis-Hastings accept-reject step, which will be discussed in Sect. 5.3. The tradeoff between bias and variance is considered in more detail in Sect. 6.

Introducing a Metropolis-Hastings Accept-Reject Step
When sampling from a distribution π it is natural to introduce a Metropolis-Hastings (MH) accept-reject step, rather than use the Markov chain obtained directly from an Euler-Maruyama discretization of the SDE (6). Given a current state X (n) , a next state is proposed according to the Euler-Maruyama discretisation of (6): The proposed state is then accepted i.e.X (n+1) :=X with probability where p(·, x) is the proposal density corresponding to (60). By introducing this accept-reject mechanism, the resulting chain X (n) is guaranteed to have a stationary distribution which is exactly π, independently of t. Thus, introducing a Metropolis-Hastings accept-reject step allows for far larger stepsizes to be used, while still preserving the correct invariant distribution of the chain, eliminating any bias arising from the discretisation of the SDE, making it beneficial both in terms of computational performance and stability. A natural question is whether an MH chain using a proposal distribution based on the SDE (6) with antisymmetric drift will inherit the superior mixing properties of the nonreversible diffusion process. As the MH algorithm works by enforcing the detailed balance of the chain X (n) with respect to the distribution π, we expect that any benefits of the antisymmetric drift term will be negated when introducing this accept-reject step. To test this, we repeat the numerical experiment of Sect. 5.2 using the MH algorithm using the Euler discretisation of (6) as a proposal scheme, for various values of α. The effect of introducing this acceptreject step to the nonreversible diffusion is evident from Fig. 5. While the accept-reject step removes any bias due to discretisation error, as is evident from Fig. 5b, the asymptotic variance actually increases as α increases. This is due to the fact that for large α, proposals are more likely to be rejected as they are far away from the current state.

Dimer in a Solvent
We now test the effect of adding a nonreversible perturbation to a test model from molecular dynamics, as described in [38,Sect. 1.3.2.4]. We consider a system composed of N particles P 1 , . . . , P N in a two-dimensional periodic box of side length L. Particles P 1 and P 2 are assumed to form a dimer pair, in a solvent comprising the particles P 3 , . . . , P N . The solvent particles interact through a truncated Lennard-Jones potential: where r is the distance between two particles, and σ are two positive parameters, and r 0 = 2 1 6 σ . The interaction potential between the dimer pair is given by a double-well potential where h and w are two positive parameters. The total energy of the system is given by with q = (q 1 , . . . , q N ) ∈ (LT) d N , where q i denotes the position of particle P i . The potential V S has two energy minima at r = r 0 (corresponding to the compact state), and at r = r 0 +2w (corresponding to the stretched state). The energy barrier separating the two states is h. Define ξ(q) to be This reaction coordinate describes the transition from the compact state, ξ(q) = 0 to the stretched state ξ(q) = 1. The standard overdamped reversible dynamics to sample from the stationary distribution π(q) ∝ exp(−βV (q)) is given by where β is the inverse temperature. Our objective is to compute the average reaction coordinate E π ξ(q). We introduce an antisymmetric drift term as follows where J ∈ R 2N ×2N is an antisymmetric matrix. We consider two types of nonreversible perturbations, determined by antisymmetric matrices, J 1 and J 2 . The first matrix J 1 is the block-circulant matrix defined by: where I 2 and O 2 denote the 2 × 2 identity and zero matrix, respectively. The second matrix is given by where R is the following rotation matrix on R 4×4 : and O 4 is the 4 × 4 zero matrix. The effect, in this case, is to apply the antisymmetric transformation only on the first two coordinates, which correspond to the positions of the particles composing the dimer.
Using an Euler-Maruyama discretisation of (63), samples are generated for N = 16 particles, with parameter values β = 1.0, σ = 1.0, = 1.0, h = 1.0, w = 0, and α ∈ {0, 2, 4, 6, 8, 10}. For each value of α, a single realisation, starting from 0 is simulated for 10 10 timesteps of size t = 10 −5 , so that the total running time is T = 10 5 . The asymptotic variance is approximated from a single realisation of the process using a batchmeans estimator (see [1,IV .5]). In Fig. 6 we plot the value of the generated estimator, along with the asymptotic variance for different values of α and for J 1 and J 2 . The first perturbation provides marginally lower asymptotic variance for this particular observable, giving rise to a 66 % decrease, as opposed to a 50 % decrease for J 2 , increasing α from 0 to 10. This decrease in asymptotic variance, comes at the cost of increased computation time due to having to reduce the stepsize to compensate for the discretisation error and numerical instability caused by the large drift. While this may appear as a negative result, since the antisymmetric drift terms were chosen arbitrarily, no claim is made about optimality. An important issue is whether the increase in computational cost outweighs the benefits. This issue will be studied more carefully in Sect. 6.

The Computational Cost of Nonreversible Langevin Samplers
As observed Sect. 5, while increasing α is guaranteed to decrease the asymptotic variance σ 2 f (α) for the estimator π T ( f ), this will also give rise to an increase in the discretisation error arising from the particular discretisation being used. Moreover, as α increases, the SDE Fig. 6 Value of the estimator π T (ξ ) for T = 10 5 with corresponding 95 % confidence intervals, for the dimer model, generated from a single realisation of (63). The right plot shows the asymptotic variance of the estimator over varying α, with antisymmetric drifts defined by J 1 ∇V (x) and J 2 ∇V (x), respectively (Color figure online) (6) becomes more and more stiff, to the extent that the discretisation becomes numerically unstable unless the stepsize is chosen to be accordingly small. As a result, any discretisation {X (n) } N n=1 will require smaller timesteps to guarantee that the stationary distribution of X (n) is sufficiently close to π(x). This tradeoff between computational cost and asymptotic variance of the estimator must be taken into consideration when comparing reversible to nonreversible diffusions.
A rigorous error analysis of the long time average estimator was carried out in [47]. In this paper careful estimates of the mean square error were derived for discretisations of a general overdamped Langevin diffusion on the unit torus T d , see also [70]. In particular, in [47,Theorem 5.2] it is shown that the mean squared error can be bounded as follows where C is a positive constant independent of t and N , which depends on the coefficients of the SDE and the observable f . This estimate makes explicit the tradeoff between discretisation error and sampling error. For a fixed computational budget N , the right hand side of (64) is minimized when t ∝ (N ) − 1 3 . For an SDE of the form (6), we expect that the constant C will increase with α. Identifying the correct scaling of the error with respect to α is an interesting problem that we intend to study.
To obtain a clearer idea of the bias variance tradeoff we compute the mean-square error for the Euler-Maruyama discretisation for two particular examples. In Fig. 7a, we consider the warped Gaussian distribution defined by (58) and the observable f (x) = |x| 2 . A value for π( f ) is obtained by integrating R d f (x)π(dx) numerically, using a globally adaptive quadrature scheme to obtain an approximation with error less than 10 −12 . In Fig. 7a we plot the relative mean-squared-error defined by Err N , t [ f ]/π( f ) 2 for an Euler-Maruyama discretisation of (6), for timestep t in the interval [2 −5 , 1]. The total number of timesteps is kept fixed at N = 10 6 . For each value of α, the mean square error is approximated over an ensemble of 256 independent realisations. Missing points indicate finite time blowup of the discretized diffusion. The dashed line denotes the MSE generated from the corresponding MALA sampler, namely an Euler-Maruyama discretisation of the reversible diffusion with an added Metropolis-Hastings accept-reject step. We note that both the Euler-Maruyama discretisation and the MALA sampler require one evaluation of the gradient term ∇ log π per timestep, so that comparing an ergodic average obtained from 10 6 steps of each scheme is fair. A trade-off between discretisation error and variance is evident from Fig. 7a, and is consistent with the error estimate (64). We observe that the nonreversible Langevin sampler outperforms the reversible Langevin sampler by an order of magnitude, with the lowest MSE attained when α = 10. As α is increased beyond this point, the discretisation error is balanced by the decrease in variance, and we observe no further gain in performance. Nonetheless, despite the fact that the MALA scheme has no bias, the nonreversible sampler, with α = 5, outperforms MALA (in terms of MSE) by a significant factor of 8.8.
We repeat this numerical experiment for the target distribution π given by a standard Gaussian distribution in R d and observable f (x) = x 2 + x 3 . In this case π( f ) is exacty 0. We use the linear diffusion specified by (41) where the antisymmetric matrix J is given in (52), which is optimal for this observable. We plot the (absolute) MSE for the estimator π T ( f ) in Fig. 8a. While the smallest MSE is attained by the nonreversible Langevin sampler, when α = 25, the increase in performance is only marginal. This is due to the fact that increasingly smaller timesteps must be taken to ensure that the EM approximation does not blow up. Indeed, the α = 25 sampler would not converge to a finite value for t greater than 10 −3 , while the reversible sampler (α = 0) and the MALA scheme were accurate even for timesteps of order 1.
From both examples it is clear that managing the numerical stability and discretisation error of the skew-symmetric drift term is essential for any practical implementation of the nonreversible sampling scheme. This suggests that using higher order and/or more stable numerical integrators to compute long time averages would be beneficial to eliminate the bias arising from discretisation, as well as permit the use of larger timesteps. Naturally, such schemes would require multiple evaluations of ∇ log π at each step, thus it is possible that the additional computational cost offsets any performance gain. In the remainder of this section, we will perform the same numerical experiments using an integrator based on a Strang splitting [35,67] of the stochastic reversible and deterministic nonreversible dynamics. The reversible part will be simulated using a standard MALA scheme and the nonreversible flow using an appropriate higher-order integrator. Indeed, denote by r,t the evolution of the reversible SDE (3) from time 0 to t, and let n,t (x) denote the flow map corresponding to the ODE:ż (t) = γ (z(t)), z(0) = x.
We shall consider an integrator based on the following map from time t to t + t: between the results obtained by MALA and the splitting scheme, we note that while a careful implementation of MALA requires only one evaluation of ∇ log π per timestep, the splitting scheme requires six evaluations of ∇ log π per timestep (naively a single timestep would require two evaluations for each reversible substep and four evaluations for the nonreversible substep, however we can reuse two evaluations of ∇ log π between the steps). Thus, we shall compare the MSE obtained from trajectories of 10 6 timesteps of the nonreversible sampler with 6 · 10 6 timesteps of the corresponding MALA scheme, for stepsizes ranging from 10 −5 to 1. The results for the warped Gaussian distribution in R 2 and standard Gaussian in R 3 are plotted in Figs. 7b and 8b, respectively. Note that we omit the α = 0 case since, in this case, the splitting scheme reduces to standard the MALA scheme. We observe that with this splitting scheme, the nonreversible sampler outperforms MALA by a factor of 13 for the warped Gaussian model, and by a factor of 20 for the standard Gaussian model. The benefits of the splitting scheme appear to be twofold: firstly the integrator is more stable, in both models, the long time simulation of X γ t did not blow up, even for large values of α, and for t = 0.1. Moreover, compared to the corresponding Euler-Maruyama discretisation, the MSE is consistently an order of magnitude less.
While this splitting scheme is only a first step into properly investigating appropriate integrators for nonreversible Langevin schemes, the above numerical experiments demonstrate clearly that there is a significant benefit in doing so, which motivates future investigation.

Conclusions and Further Work
In this paper we have presented a detailed analytical and numerical study of the effect of nonreversible perturbations to Langevin samplers. In particular, we have focused on the effect on the asymptotic variance of adding a nonreversible drift to the overdamped Langevin dynamics. Our theoretical analysis, presented for diffusions with periodic coefficients and for diffusions with linear drift for which a complete analytical study can be performed, and our numerical investigations on toy models clearly show that a judicious choice of the nonreversible drift can lead to a substantial reduction in the asymptotic variance. On the other hand, as observed from the dimer model example in Sect. 5.4, an arbitrary choice of nonreversible drift will not always give rise to significant improvement. We have also presented a careful study of the computational cost of the algorithm based on a nonreversible Langevin sampler, in which the competing effects of reducing the asymptotic variance and of increasing the stiffness due to the addition of a nonreversible drift are monitored. The main conclusions that can be drawn from our numerical experiments is that a nonreversible Langevin sampler with close-to-the-optimal choice of the nonreversible drift significantly outperforms the (reversible) Metropolis-Hastings sampler.
There are many open problems that need to be studied further: (1) The effect of using degenerate, hypoelliptic diffusions for sampling from a given distribution.
(2) Combining the use of nonreversible Langevin samplers with standard variance reduction techniques such as the zero variance reduction MCMC methodology [14]. (3) Optimizing Langevin samplers within the class of reversible diffusions. (4) The development of nonreversible Metropolis-Hastings algorithms based on the above techniques, possibly related to approach described in [9]. (5) The development and analysis of numerical schemes specifically designed to simulate nonreversible Langevin diffusions.
All these topics are currently under investigation.