Using Perturbed Underdamped Langevin Dynamics to Efficiently Sample from Probability Distributions

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.


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 π on R d , i.e. we wish to evaluate integrals of the form: (1) 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 [19,32,52], where a Markov process (X t ) t≥0 is constructed which is ergodic with respect to the probability measure π. Then, defining the long-time average for f ∈ L 1 (π), the ergodic theorem guarantees almost sure convergence of the long-time average π T ( f ) to π( 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≥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≥0 of the following stochastic differential equation (SDE): where V = − log π is the potential associated with the smooth positive density π. Under appropriate assumptions on V , i.e. on the measure π(dx), the process (X t ) t≥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≥0 = (q t , p t ) t≥0 defined on the extended space (phase space) R d × R d by the following pair of coupled SDEs: where the mass and friction tensors M and Γ are assumed to be symmetric positive definite matrices. It is well-known [36,46] that (q t , p t ) t≥0 is ergodic with respect to the measure π := π ⊗ N (0, M), having density with respect to the Lebesgue measure on R 2d given by where Z is a normalization constant. Note that π has marginal π with respect to p and thus for functions f ∈ L 1 (π), we have that 1 t t 0 f (q t ) dt → π( 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 [12].
As there is a lot of freedom in choosing the dynamics in (2), see the discussion in Sect. 2, it is desirable to choose the diffusion process (X t ) t≥0 in such a way that π T ( f ) can provide a good estimation of π( 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 Sect. 5 and [14,Sect. 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≥0 so that this rate of convergence is maximised. Conditions on the potential V which guarantee exponential convergence to equilibrium, both in L 2 (π) and in relative entropy can be found in [7,39,54]. In the case when the target measure π 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 [38,43,44]. 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 [1,2].
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 π and the observable f , the estimator π T ( f ) satisfies a central limit theorem (CLT) [31], that is, where σ 2 f < ∞ is the asymptotic variance of the estimator π T ( f ). The asymptotic variance characterises the magnitude of fluctuations of π T ( f ) around π( f ). Consequently, another natural objective is to choose the process (X t ) t≥0 such that σ 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 [31] − 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 [14], see also [23], and it will also be used in this paper.
Other measures of performance have also been considered. For example, in [50,51], performance of the estimator is quantified in terms of the rate functional of the ensemble measure 1 t t 0 δ X (t) (dx). See also [28] 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 [18] (see the discussion in Sect. 2) and the nonreversible Langevin samplers [20,21]. As a particular example of the general framework that was introduced in [18], we mention the preconditioned overdamped Langevin dynamics dX t = −P∇V (X t ) dt + √ 2P dW t , that was presented in [4]. There, the long-time behaviour of as well as the asymptotic variance of the corresponding estimator π 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 [6].
The nonreversible overdamped Langevin dynamics where the vector field γ satisfies ∇ · (πγ ) = 0 is ergodic (but not reversible) with respect to the target measure π for all choices of the divergence-free vector field γ . The asymptotic behaviour of this process was considered for Gaussian diffusions in [20], where the rate of convergence of the covariance to equilibrium was quantified in terms of the choice of γ . This work was extended to the case of non-Gaussian target densities, and consequently for nonlinear SDEs of the form (7) in [21]. The problem of constructing the optimal nonreversible perturbation, in terms of the L 2 (π) spectral gap for Gaussian target densities was studied in [34] see also [55]. Optimal nonreversible perturbations with respect to miniziming the asymptotic variance were studied in [14,23]. 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 (7)], the nonreversible Langevin sampler (7) is never worse that its reversible counterpart (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 γ (x). The goal of this paper is to extend the analysis presented in [14,34] by introducing the following modification of the standard underdamped Langevin dynamics: where M, Γ ∈ R d×d are constant strictly positive definite matrices, μ and ν are scalar constants and J 1 , J 2 ∈ R d×d are constant skew-symmetric matrices. As demonstrated in Sect. 2, the process defined by (8) will be ergodic with respect to the Gibbs measure π 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 (8) 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 (8) in order to study the problem of diffusion bridge sampling.
The rest of the paper is organized as follows. In Sect. 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 Sect. 3 we study qualitative properties of the perturbed underdamped Langevin dynamics (8) including exponentially fast convergence to equilibrium and the overdamped limit. In Sect. 4 we study in detail the performance of the Langevin sampler (8) for the case of Gaussian target distributions. In Sect. 5 we introduce a numerical scheme for simulating the perturbed dynamics (8) 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.

Background and Preliminaries
In this section we consider estimators of the form (2) where (X t ) t≥0 is a diffusion process given by the solution of the following Itô SDE: with drift coefficient a : R d → R d and diffusion coefficient b : R d → R d×m both having smooth components, and where (W t ) t≥0 is a standard R m -valued Brownian motion. Associated with (9) is the infinitesimal generator L formally given by where Σ = bb , ∇∇ f denotes the Hessian of the function f and : denotes the Frobenius inner product. In general, Σ is nonnegative definite, and could possibly be degenerate. In particular, the infinitesimal generator (10) need not be uniformly elliptic. To ensure that the corresponding semigroup exhibits sufficient smoothing behaviour, we shall require that the process (9) is hypoelliptic in the sense of Hörmander. If this condition holds, then irreducibility of the process (X t ) t≥0 will be an immediate consequence of the existence of a strictly positive invariant distribution π(x)dx, see [30]. Suppose that (X t ) t≥0 is nonexplosive. It follows from the hypoellipticity assumption that the process (X t ) t≥0 possesses a smooth transition density p(t, x, y) which is defined for all t ≥ 0 and x, y ∈ R d , [5,Theorem VII.5.6]. The associated strongly continuous Markov semigroup (P t ) t≥0 is defined by P t f (x) = R d p(t, x, y) f (y) dy. Suppose that (P t ) t≥0 is invariant with respect to the target measure π, i.e. R d P t f (x)π(dx) = R d f (x)π(dx) for t ≥ 0 and all bounded continuous functions f . Then (P t ) t≥0 can be extended to a positivity preserving contraction semigroup on L 2 (π) which is strongly continuous. Moreover, the infinitesimal generator corresponding to (P t ) t≥0 is given by an extension of (L, C 2 c (R d )), also denoted by L.
Due to hypoellipticity and invariance with respect to (P t ) t≥0 , the probability measure π on R d has a smooth density with respect to the Lebesgue measure. If this density is strictly positive, it follows that π is necessarily the unique invariant distribution. Slightly abusing the notation, we will denote both the measure and its density by π. Furthermore, we will denote by L 2 (π) be the Hilbert space of π-square integrable functions equipped with inner product ·, · L 2 (π ) and norm · L 2 (π ) .

A General Characterisation of Ergodic Diffusions
A natural question is what conditions on the coefficients a and b of (9) are required to ensure that (X t ) t≥0 is invariant with respect to the distribution π(x) dx. The following result provides a necessary and sufficient condition for a diffusion process to be invariant with respect to a given target distribution.
Theorem 1 Consider a diffusion process (X t ) t≥0 on R d defined by the unique, nonexplosive solution to the Itô SDE (9) with drift a ∈ C 1 (R d ; R d ) and diffusion coefficient b ∈ C 1 (R d ; R d×m ). Then (X t ) t≥0 is invariant with respect to π if and only if a = Σ∇ log π + ∇ · Σ + γ, where Σ = bb and γ : R D → R D is a continuously differentiable vector field satisfying ∇ · (πγ ) = 0. (12) If additionally γ ∈ L 1 (π), then there exists a skew-symmetric matrix function C : R d → R d×d such that γ = 1 π ∇ · (πC) . In this case the infinitesimal generator can be written as an L 2 (π)-extension of The proof of the first part of this result can be found in [46,Chap. 4]; similar versions of this characterisation can be found in [54] and [21]. For the existence of the skew-symmetric matrix C see, e.g., [16,Sec.4,Prop. 1]. See also [37].
Remark 1 If (11) holds and L is hypoelliptic it follows immediately that (X t ) t≥0 is ergodic with unique invariant distribution π(x) dx (see [30]).
More generally, we can consider Itô diffusions in an extended phase space: where (W t ) t≥0 is a standard Brownian motion in R N , N ≥ d. This is a Markov process with generator where Σ(z) = σ σ T (z). We will consider dynamics (Z t ) t≥0 that is ergodic with respect to There are various well-known choices of dynamics which are invariant (and indeed ergodic) with respect to the target distribution π(x)dx.
1. Choosing b = I and γ = 0 we immediately recover the overdamped Langevin dynamics (3). 2. Choosing b = I , and γ = 0 such that (12) holds gives rise to the nonreversible overdamped equation defined by (7). As it satisfies the conditions of Theorem 1, it is ergodic with respect to π. In particular choosing γ (x) = J ∇V (x) for a constant skew-symmetric matrix J we obtain which has been studied in previous works. 3. Given a target density π > 0 on R d , if we consider the augmented target density π on R 2d given in (5), then choosing where M and Γ are positive definite symmetric matrices, the conditions of Theorem 1 are satisfied for the target density π. The resulting dynamics (q t , p t ) t≥0 is determined by the underdamped Langevin equation (4). It is straightforward to verify that the generator is hypoelliptic, [35, Sec 2.2.3.1], and thus (q t , p t ) t≥0 is ergodic. 4. More generally, consider the augmented target density π on R 2d as above, and choose where μ and ν are scalar constants and J 1 , J 2 ∈ R d×d are constant skew-symmetric matrices. With this choice we recover the perturbed Langevin dynamics (8). It is straightforward to check that (17) satisfies the invariance condition (12), and thus Theorem 1 guarantees that (8) is invariant with respect to π. 5. In a similar fashion, one can introduce an augmented target density on R (m+2)d , with where λ i ∈ R and α i > 0, for i = 1, . . . , m. The resulting process (9) is given by where (W j t ) t≥0, j=1,...,m are independent R d -valued Brownian motions. This process is ergodic with unique invariant distribution π, and under appropriate conditions on V , converges exponentially fast to equilibrium in relative entropy [42]. Equation (18) is a Markovian representation of a generalised Langevin equation of the form where N (t) is a mean-zero stationary Gaussian process with autocorrelation function Then choosing b = I D×D and γ = 0 we obtain the dynamics

Comparison Criteria
For a fixed observable f , a natural measure of accuracy of the estimator π T ( f ) = t −1 t 0 f (X s ) ds is the mean square error (MSE) defined by where E x denotes the expectation conditioned on the process Here μ( f, T ) measures the bias of the estimator π T ( f ) and σ 2 ( f, T ) measures the variance of fluctuations of π T ( f ) around the mean. The speed of convergence to equilibrium of the process (X t ) t≥0 will control both the bias term μ( f, T ) and the variance σ 2 ( f, T ). To make this claim more precise, suppose that the semigroup (P t ) t≥0 associated with (X t ) t≥0 decays exponentially fast in L 2 (π), i.e. there exist constants λ > 0 and C ≥ 1 such that Remark 2 If (21) holds with C = 1, this estimate is equivalent to − L having a spectral gap in L 2 (π). 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 [54].
The following lemma characterises the decay of the bias μ( f, T ) as T → ∞ in terms of λ and C. The proof can be found in [41].

Lemma 1
Let (X t ) t≥0 be the unique, non-explosive solution of (9), such that X 0 ∼ π 0 π and dπ 0 dπ ∈ L 2 (π), where dπ 0 dπ denotes the Radon-Nikodym derivative of π 0 with respect to π. Suppose that the process is ergodic with respect to π such that the Markov semigroup (P t ) t≥0 satisfies (21). Then for f ∈ L ∞ (π), The study of the long time behaviour of the variance σ 2 ( f, T ) involves deriving a central limit theorem for the additive functional t 0 f (X t )−π( f ) dt. As discussed in [13], we reduce this problem to proving well-posedness of the Poisson equation The only complications in this approach arise from the fact that the generator L need not be symmetric in L 2 (π) 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. Again, the proof can be found in [41].

Lemma 2
Let (X t ) t≥0 be the unique, non-explosive solution of (9) with smooth drift and diffusion coefficients, such that the corresponding infinitesimal generator is hypoelliptic. Syppose that (X t ) t≥0 is ergodic with respect to π and moreover, (P t ) t≥0 decays exponentially fast in L 2 (π) as in (21). Then for all f ∈ L 2 (π), there exists a unique mean zero solution χ to the Poisson equation (22).
where σ 2 f is the asymptotic variance defined by 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 ∈ L 2 (π) satisfying π( f ) = 0, simplifying expressions (22) and (23). The corresponding subspace of L 2 (π) will be denoted by L 2 0 (π). If the exponential decay estimate (21) is satisfied, then Lemma 2 shows that − L is invertible on L 2 0 (π), so we can express the asymptoptic variance as We note that the constants C and λ appearing in the exponential decay estimate (21) also control the speed of convergence of σ 2 ( f, T ) to zero. Indeed, it is straightforward to show that if (21) is satisfied, then the solution χ of (22) satisfies Lemmas 1 and 2 would suggest that choosing the coefficients Σ and γ to optimize the constants C and λ in (26) would be an effective means of improving the performance of the estimator π 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 [20,21,34], maximising the speed of convergence to equilibrium is a delicate task. As the leading order term in M SE( f, T ), it is typically sufficient to focus specifically on the asymptotic variance σ 2 f and study how the parameters of the SDE (9) can be chosen to minimise σ 2 f . This study was undertaken in [14] for processes of the form (7).

Perturbation of Underdamped Langevin Dynamics
The primary objective of this work is to compare the performances of the perturbed underdamped Langevin dynamics (8) and the unperturbed dynamics (4) according to the criteria outlined in Sect. 2.3 and to find suitable choices for the matrices J 1 , J 2 , M and Γ that improve the performance of the sampler. We begin our investigations of (8) by establishing ergodicity and exponentially fast return to equilibrium, and by studying the overdamped limit of (8).
As the latter turns out to be nonreversible and therefore in principle superior to the usual overdamped limit (3), e.g. [21], 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) = 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: Under the above choices (27), we show that the large perturbation limit lim μ→∞ σ 2 f exists and is finite and we provide an explicit expression for it (see Corollary 4). 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 (γ > √ 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 · K q + l · q + C: see Theorem 3. In fact, we conjecture that this statement is true for arbitrary observables f ∈ L 2 (π), but we have not been able to prove this. The dynamics (8) [used in conjunction with the conditions (27)] proves to be especially effective when the observable is antisymmetric (i.e. when it is invariant under the substitution q → −q) or when it has a significant antisymmetric part. In particular, in Proposition 3 we show that under certain conditions on the spectrum of J 1 , for any antisymmetric observable f ∈ L 2 (π) it holds that lim μ→∞ σ 2 f = 0. Numerical experiments and analysis show that departing significantly from (27) in fact possibly decreases the performance of the sampler. This is in stark contrast to (7), where it is not possible to increase the asymptotic variance by any perturbation. For that reason, until now it seems practical to use (8) 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 π 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 (27).

Remark 3
In [34, Rem. 3] another modification of (4) was suggested (albeit with the simplifications Γ = γ · I and M = I ): J again denoting an antisymmetric matrix. However, under the change of variables p → (1 + J )p the above equations transform into Since any observable f depends only on q (the p-variables are merely auxiliary), the estimator π 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, (28) reduces to the underdamped Langevin dynamics (4) and does not represent an independent approach to sampling. Suitable choices of M and Γ will be discussed in Sect. 4.5.

Properties of Perturbed Underdamped Langevin Dynamics
In this section we study some of the properties of the perturbed underdamped dynamics (8).
First, note that its generator is given by decomposed into the perturbation L pert and the unperturbed operator L 0 , which can be further split into the Hamiltonian part L ham and the thermostat (Ornstein-Uhlenbeck) part L therm , see [35,36,46].
Proof The proof consists of verifying the conditions of Hörmander's Theorem for the generator (30) and can be found in [41].
An immediate corollary of this result and of Theorem 1 is that the perturbed underdamped Langevin process (8) is ergodic with unique invariant distribution π given by (5).
As explained in Sect. 2.3, the exponential decay estimate (21) is crucial for our approach, as in particular it guarantees the well-posedness of the Poisson equation (22). From now on, we will therefore make the following assumption on the potential V, required to prove exponential decay in L 2 (π): Assumption 1 Assume that the Hessian of V is bounded and that the target measure π(dq) = 1 Z e −V dq satisfies a Poincare inequality, i.e. there exists a constant ρ > 0 such that Sufficient conditions on the potential so that Poincaré's inequality holds, e.g. the Bakry-Emery criterion, are presented in [7].
Theorem 2 Under Assumption 1 there exist constants C ≥ 1 and λ > 0 such that the semigroup (P t ) t≥0 generated by L satisfies exponential decay in L 2 (π) as in (21).
Proof The proof uses the machinery of hypocoercivity developed in [54] and can be found in [41]. Using the framework of [15], we conjecture that the assumption on the boundedness of the Hessian of V can be substantially weakened and more quantitative decay estimates (in particular with respect to μ and ν) can be obtained. This approach has recently been successfully applied to equilibrium and nonequilibirum Langevin dynamics, see [27,53]. We leave this work track for future study.

The Overdamped Limit
In this section we develop a connection between the perturbed underdamped Langevin dynamics (8) and the nonreversible overdamped Langevin dynamics (7). The analysis is very similar to the one presented in [35, Sect. 2.2.2] and we will be brief. For convenience in this section we will perform the analysis on the d-dimensional torus T d ∼ = (R/Z) d , i.e. we will assume q ∈ T d . Consider the following scaling of (8): valid for the small mass/small momentum regime M → 2 M, p t → p t . Equivalently, those modifications can be obtained from subsituting Γ → −1 Γ and t → −1 t, and so in the limit as → 0 the dynamics (32) describes the limit of large friction with rescaled time. It turns out that as → 0, the dynamics (32) converges to the limiting SDE The following proposition makes this statement precise.
Proof The proof follows standard arguments (see for instance [46]) and can be found in [41]. By a more refined analysis, it is also possible to get information on the rate of convergence; see, e.g. [48,49].

Remark 4
The overdamped limit (33) respects the invariant distribution, in the sense that it is ergodic with respect to π(dq) = 1 Z e −V dq.
The limiting SDE (33) is nonreversible due to the term −μJ 1 ∇ q V (q t )dt and also because the matrix (ν J 2 + Γ ) −1 is in general neither symmetric nor antisymmetric. This result, together with the fact that nonreversible perturbations of overdamped Langevin dynamics of the form (7) are by now well-known to have improved performance properties, motivates further investigation of the dynamics (8).

Sampling from a Gaussian Distribution
In this section we study in detail the performance of the Langevin sampler (8) 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.

Unit Covariance: Small Perturbations
In our study of the dynamics given by (8) we first consider the simple case when V (q) = 1 2 |q| 2 , i.e. the task of sampling from a Gaussian measure with unit covariance. We will assume M = I , Γ = γ 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 μ and ν). Our first result concerns the asymptotic variance for linear and quadratic observables for small perturbations of equal strength (μ = ν). For sufficiently strong damping (γ > √ 2) always leads to an improvement in asymptotic variance under the nondegeneracy conditions [J, K ] = 0 and l / ∈ ker J :

Theorem 3 Consider the dynamics
with γ > √ 2 and an observable of the form f (q) = q · K q + l · q + C, where K ∈ R d×d sym , l ∈ R d and C ∈ R. If at least one of the conditions [J, K ] = 0 and l / ∈ ker J is satisfied, then the asymptotic variance of the unperturbed sampler is at a local maximum independently of K and J (and γ , as long as γ > √ 2), i.e.
Proof The dynamics (34) are of Ornstein-Uhlenbeck type, i.e. we can write with X = (q, p) T , and (W t ) t≥0 denoting a standard Wiener process on R 2d . The generator of (35) is then given by According to Lemma 2, the asymptotic variance can be expressed as By calculations similar to those in [14,Sect. 4], χ is given by using the notationsK The asymptotic variance is then given by Taking derivatives of 38 and solving the ensuing matrix equations, it is possible to obtain explicit expressions for ∂ μ C| μ=0 , ∂ 2 μ C| μ=0 , ∂ μ D| μ=0 and ∂ 2 μ C| μ=0 as detailed in [41]. We obtain Numerical examples show that the conditions γ > √ 2 and μ = ν are indeed necessary for the conclusions of Theorem 3 to hold (an explanation in terms of the spectrum of the generator will be provided in Sect. 4.2). In particular, an unfortunate choice of the perturbations will actually increase the asymptotic variance of the dynamics.
Let us illustrate this by plotting the asymptotic variance as a function of the perturbation strength μ (see Fig. 1), making the choices d = 2, l = (1, 1) T , The asymptotic variance has been computed according to (38) and (40). Going beyond the results of this section, the graphs give an impression of the behaviour of the asymptotic variace for large values of μ, discussed further in Sect. 4.3. Figure 1a, b, c show the asymptotic variance associated with the quadratic observable f (q) = q · K q. In accordance with Theorem 3, the asymptotic variance is at a local maximum at zero perturbation in the case μ = ν (see Fig. 1a). For increasing perturbation strength, the graph shows monotone decay for μ → ∞ (this limiting behaviour will be explored analytically in Sect. 4.3). If the condition μ = ν is only approximately satisfied (Fig. 1b), 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 μ. If the perturbations are opposed (μ = −ν), it is possible for certain observables that the unperturbed dynamics represents a global minimum. Such a case is observed in Fig. 1c. In Fig. 1d, e the observable f (q) = l · q is considered. If the damping is sufficiently strong (γ > √ 2), the unperturbed dynamics is at a local maximum of the asymptotic variance (Fig. 1d). Furthermore, the asymptotic variance approaches zero as μ → ∞ (for a theoretical explanation see again Sect. 4.3). The graph in Fig. 1e shows that the assumption of γ not being too small cannot be dropped from Theorem 3. Even in this case though the example shows decay of the asymptotic variance for large values of μ.

Exponential Decay Rate
Let us denote by λ * the optimal exponential decay rate in (21), i.e. λ * = sup{λ > 0 | There exists C ≥ 1 such that (21) holds}. (42) Note that λ * is well-defined and positive by Theorem 2. We also define the spectral bound of the generator L by In [38] it is proven that the Ornstein-Uhlenbeck semigroup (P t ) t≥0 considered in this section is differentiable (see Proposition 2.1). In this case (see Corollary 3.12 of [17]), it is known that the exponential decay rate and the spectral bound coincide, i.e. λ * = s(L), whereas in general only λ * ≤ s(L) holds. In this section we will therefore analyse the spectral properties of the generator L in the Gaussian case. In particular, this leads to some intuition of why choosing equal perturbations (μ = ν) is crucial for the performance of the sampler.
In [38] (see also [43]), it was proven that the spectrum of L as in (36) in L 2 ( π) is given by Note that σ (L) only depends on the drift matrix B. In the case where μ = ν, the spectrum of B can be computed explicitly.

Lemma 4 Assume μ = ν. Then the spectrum of B is given by
Proof We will compute σ B − γ 2 I and then use the identity σ ( where I is understood to denote the identity matrix of appropriate dimension. The above quantity is zero if and only if Together with (4.2), the claim follows.
Using formula (45), in Fig. 2a we show a sketch of the spectrum σ (−L) for the case of equal perturbations (μ = ν) with the convenient choices n = 1 and γ = 2. Of course, the eigenvalue at 0 is associated to the invariant measure since L † π = 0. The arrows indicate the movement of the eigenvalues as the perturbation μ increases in accordance with Lemma 4. Clearly, the spectral bound of L is not affected by the perturbation. Note that the eigenvalues on the real axis stay invariant under the perturbation. The subspace of L 2 0 ( π) associated to those will turn out to be crucial for the characterisation of the limiting asymptotic variance as μ → ∞ (see Remark 10).
To illustrate the suboptimal properties of the perturbed dynamics when the perturbations are not equal, we plot the spectrum of the drift matrix σ (B) in the case when the dynamics is only perturbed by the term ν J 2 pdt (i.e. μ = 0) for n = 2, γ = 2 and (see Fig. 2b). Note that the full spectrum σ (−L) can be inferred from (44). For ν = 0 we have that the spectrum σ (B) only consists of the (degenerate) eigenvalue 1. For increasing ν, the figure shows that the degenerate eigenvalue splits up into four eigenvalues, two of which get closer to the imaginary axis as ν increases, leading to a smaller spectral bound and therefore to a decrease in the speed of convergence to equilibrium. Figure 2a, b give an intuitive explanation of why the fine-tuning of the perturbation strengths is crucial. We close this subsection by providing autocorrelation plots (see Fig. 3) for the linear observable considered in Fig. 1d (with a friction coefficient of γ = 2.5). It is well-known that the asymptotic variance is given by the integrated autocorrelation function (see e.g. Proposition IV 1.3 in [3]), 2 Notice that γ 2 2 − 1 is understood to be a complex number for γ < 2. Comparing Fig. 3a, b yields additional insight into the mechanics of the variance reduction: the increase of the imaginary part of the eigenvalues of L (as indicated in Fig. 2a) leads to oscillations in the autocorrelation function and therefore to cancellations in (47). A similar effect has already been observed in [50] for the nonreversible overdamped Langevin dynamics (15).

Unit Covariance: Large Perturbations
In the previous subsection we observed that for the particular perturbation J 1 = J 2 and μ = ν [see equation (34)] the perturbed Langevin dynamics demonstrated an improvement in performance for μ in a neighbourhood of 0, when the observable is linear or quadratic.
Recall that this dynamics is ergodic with respect to a standard Gaussian measureπ on R 2d with marginal π with respect to the q-variable. As before, we shall consider only observables that do not depend on p. Moreover, we assume without loss of generality that π( f ) = 0. For such observables we will write f ∈ L 2 0 (π) and consider the canonical embedding L 2 0 (π) ⊂ L 2 (π). We emphasize that L 2 0 (π) consists of functions that only depend on q, whereas functions in L 2 (π) may depend on both q and p.
In this subsection will analyse the asymptotic variance for large values of μ. The infinitesimal generator of (34) can be written as (48) where we have introduced the notation L pert = μA. In the sequel, the adjoint of an operator B in L 2 ( π) will be denoted by B * . In the rest of this section we will make repeated use of the Hermite polynomials invoking the notation x = (q, p) ∈ R 2d . For m ∈ N 0 define the Hilbert spaces The following result (Theorem 4) holds for operators of the form (36) providing an orthogonal decomposition of L 2 ( π) into invariant subspaces. The drift and diffusion matrices B and Q are assumed to be such that L is the generator of an ergodic stochastic process (see [2, Definition 2.1] for precise conditions).

Theorem 4 [2, Sect. 5]. The following holds:
(a) The space L 2 ( π) has a decomposition into mutually orthogonal subspaces: (c) The spectrum of L has the following decomposition: Remark 6 Note that by the ergodicity of the dynamics, ker L consists of constant functions and so ker L = H 0 . Therefore, L 2 0 ( π) has the decomposition Our first main result of this section is an expression for the asymptotic variance in terms of the unperturbed operator L 0 and the perturbation A: Proposition 2 Let f ∈ L 2 0 (π) (so in particular f = f (q)). Then the associated asymptotic variance is given by Remark 7 The proof of the preceding Proposition will show that L 2 0 + μ 2 A * A is invertible on L 2 0 ( π) and that ( To prove Proposition 2 we will make use of the generator with reversed perturbation L − = L 0 − μA and the momentum flip operator Clearly, P 2 = I and P * = P. Further properties of L 0 , A and the auxiliary operators L − and P are gathered in the following lemma: Lemma 5 For all φ, ψ ∈ C ∞ (R 2d ) ∩ L 2 ( π) the following holds: (a) The generator L 0 is symmetric in L 2 ( π) with respect to P: The perturbation A is skewadjoint in L 2 ( π): (c) The operators L 0 and A commute: (e) L and L − commute, [L, L − ]φ = 0, and the following relation holds: (f) The operators L, L 0 , L − , A and P leave the Hermite spaces H m invariant.

Remark 8
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 ).

Proof of Lemma 5
The statement (a) is well-known and its proof can be found in [35, Sect. 2.2.3.1] for instance. The claim (b) follows by noting that the flow vector field b(q, p) = (−J q, −J p) associated to A is divergence-free with respect to π, i.e. ∇·( π b) = 0. Therefore, A is the generator of a strongly continuous unitary semigroup on L 2 ( π) and hence skewadjoint by Stone's Theorem. The claims (c), (d) and (e) follow by direct computations which can be found in [41]. To prove (f) first notice that L, L 0 and L − are of the form (36) and therefore leave the spaces H m invariant by Theorem 4. It follows immediately that also A leaves those spaces invariant. The fact that P leaves the spaces H m invariant follows directly by inspection of (49) and (51).

Now we proceed with the proof of Proposition 2:
Proof of Proposition 2 Since the potential V is quadratic, Assumption 1 clearly holds and thus Lemma 2 ensures that L and L − are invertible on L 2 0 ( π) with and analogously for L −1 − . In particular, the asymptotic variance can be written as Due to the respresentation (53) and Theorem 4, the inverses of L and L − leave the Hermite spaces H m invariant. We will prove the claim from Proposition 2 under the assumption that P f = f which includes the case f = f (q). For the following calculations we will assume f ∈ H m for fixed m ≥ 1. Combining statement (f) with (a) and (e) of Lemma 5 (and noting that H m ⊂ C ∞ (R 2d ) ∩ L 2 ( π)) we see that when restricted to H m . Therefore, the following calculations are justified: where in the second line we have used the assumption P f = f and in the third line the properties P 2 = I , P * = P and Eq. (54). Since L and L − commute on H m according to Lemma 5(e),(f) we can write for the restrictions on H m , using L+L − = 2L 0 . We also have LL − = (L 0 +μA)(L 0 −μA) = L 2 0 + μ 2 A * A, since L 0 and A commute. We thus arrive at the formula Now since ( , it follows that the operator −L 0 (L 2 0 + μ 2 A * A) −1 is bounded. We can therefore extend formula (55) to the whole of L 2 ( π) by continuity, using the fact that L 2 0 ( π) = m≥1 H m .
Applying Proposition 2 we can analyse the behaviour of σ 2 f in the limit of large perturbation strength μ → ∞. To this end, we introduce the orthogonal decomposition where J q ·∇ q is understood as an unbounded operator acting on L 2 0 (π), obtained as the smallest closed extension of J q · ∇ q acting on C ∞ c (R d ). In particular, ker(J q · ∇ q ) is a closed linear subspace of L 2 0 (π). Let Π denote the L 2 0 (π)-orthogonal projection onto ker(J q · ∇ q ). We will write σ 2 f (μ) 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 Π f for the unperturbed dynamics.

Remark 9
Note that the fact that the limit exists and is finite is nontrivial. In particular, as Fig. 1b, c demonstrate, it is often the case that lim μ→∞ σ 2 f (μ) = ∞ if the condition μ = ν is not satisfied.

Remark 10
The decomposition (56) can be interpreted in terms of the spectrum σ (L) as follows: First observe that for functions f that only depend on q, f ∈ ker(J q · ∇ q ) is equivalent to f ∈ ker A. Let us denote byσ = μ∈R σ (L 0 + μA) the part of σ (L 0 ) that is not affected by the perturbation and bȳ the corresponding subspace. Then it is straightforward to see that ker(A) =Ē. 3 In Fig. 2a, σ has been highlighted by diamonds.
Proof of Theorem 5 Note that L 0 and A * 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 A * A is nonnegative self-adjoint, there exists an orthogonal decomposition L 2 0 (π) = i W i into eigenspaces of the operator −L 0 (L 2 0 + μ 2 A * A) −1 , the decomposition W i being finer then H m in the sense that every W i is a subspace of some H m . Moreover, (50) can be written as where f = i f i and f i ∈ W i . Let us assume now without loss of generality that W 0 = ker A * A, so in particular λ 0 = 0. Then clearly showing the equality in the claim. It remains to show that σ 2 Π f (0) ≤ σ 2 f (0). To see this, we write Note that since we only consider observables that do not depend on p, Π f ∈ ker(J q · ∇ q ) and (1 − Π) f ∈ i≥1 W i . Since L 0 commutes with A, it follows that (−L 0 ) −1 leaves both W 0 and i≥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 5 it follows that in the limit as μ → ∞, the asymptotic variance σ 2 f (μ) is not decreased by the perturbation if f ∈ ker(J q · ∇ q ). In fact, this result also holds true non-asymptotically, i.e. observables in ker(J q ·∇ q ) are not affected at all by the perturbation: Proof From f ∈ ker(J q · ∇ q ) it follows immediately that f ∈ ker A * A. Then the claim follows from the expression (57).
Example 1 Recall the case of observables of the form f (q) = q · K q + l · q + C with K ∈ R d×d sym , l ∈ R d and C ∈ R from Sect. 4.1. If [J, K ] = 0 and l ∈ ker J , then f ∈ ker(J q · ∇ q ) as From the preceding lemma it follows that σ 2 showing that the assumption in Theorem 3 does not exclude nontrivial cases.
The following result shows that the dynamics (34) is particularly effective for antisymmetric observables (at least in the limit of large perturbations): and assume that ker J = {0}. Furthermore, assume that the eigenvalues of J are rationally independent, i.e.

Proof of Proposition 3
The claim would immediately follow from f ∈ ker(J q ·∇) ⊥ according to Theorem 5, 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 2 that L is invertible on L 2 0 ( π) and its inverse leaves the Hermite spaces H m invariant. Consequently, the asymptotic variance of an observable f ∈ L 2 0 ( π) can be written as where Π m : L 2 0 ( π) → H m denotes the orthogonal projection onto H m . From (49) it is clear that g a is symmetric for |α| even and antisymmetric for |α| odd. Therefore, from f being antisymmetric it follows that f ∈ m≥1,m odd H m . In view of (45), ((c)) and (58) the spectrum of L |H m can be written as with appropriate real constants C α,γ ∈ R that depend on α and γ , but not on μ. For |α| = 2d j=1 α j = m odd, we have that Indeed, assume to the contrary that the above expression is zero. Then it would follow that α j = α j+d for all j = 1, . . . , d by rational independence of λ 1 , . . . , λ d and |m| would have to be even. From (60) and (61) it is clear that where B(0, r ) denotes the ball of radius r centered at the origin in C. Consequently, the spectral radius of (−L| H m ) −1 and hence (−L| H m ) −1 itself converges to zero as μ → ∞.
The result then follows from (59).

Remark 11
The idea of the preceding proof can be explained using Fig. 2a and Remark 10.
The eigenvalues in the fixed spectrumĒ (on the real axis, highlighted by diamonds) correspond to Hermite polynomials of even order. The independence condition on the eigenvalues of J prevents cancellations that would lead to fixed eigenvalues associated to Hermite polynomials of odd order. Therefore, antisymmetric observables are orthogonal toĒ = ker A.
The following corollary gives a version of the converse of Proposition 3 and provides further intuition into the mechanics of the variance reduction achieved by the perturbation.
Proof According to Theorem 5, lim μ→∞ σ 2 f (μ) = 0 implies σ 2 Π f (0) = 0. We can write and recall from the proof of Proposition 2 that (−L 0 ) −1 and (−L * 0 ) −1 leave the Hermite spaces H m invariant. Therefore ker (−L 0 ) −1 + (−L * 0 ) −1 = 0 in L 2 0 ( π), and in particular Hence, there exists a sequence (φ n ) n ∈ C ∞ c (R d ) such that J q · ∇φ n → f in L 2 (π). Taking a subsequence if necessary, we can assume that the convergence is pointwise π-almost everywhere and that the sequence is pointwise bounded by a function in L 1 (π). Since J is antisymmetric, we have that J q · ∇φ n = ∇ · (φ n J q). Now Gauss's theorem yields where n denotes the outward normal to the sphere ∂ B(0, r ). This quantity is zero due to the orthogonality of J q and n, and so the result follows from Lebesgue's dominated convergence theorem.

Optimal Choices of J for Quadratic Observables
Assume f ∈ L 2 0 (π) is given by f (q) = q · K q + l · q − Tr K , with K ∈ R d×d sym and l ∈ R d (note that the constant term is chosen such that π( f ) = 0). Our objective is to choose J in such a way that lim μ→∞ σ 2 f (μ) becomes as small as possible. To stress the dependence on the choice of J , we introduce the notation σ 2 f (μ, J ). Also, we denote the orthogonal projection onto (ker J ) ⊥ by Π ⊥ ker J . Lemma 7 (Zero variance limit for linear observables). Assume K = 0 and Π ⊥ ker J l = 0. Then lim Proof According to Theorem 5, we have to show that Π f = 0, where Π is the L 2 (π)orthogonal projection onto ker(J q · ∇). Let us thus use (63) and prove that f ∈ im(J q · ∇). Indeed, since Π ⊥ ker J l = 0, by Fredholm's alternative there exists u ∈ R d such that J u = l. Now define φ ∈ L 2 0 (π) by φ(q) = −u · q, leading to f = J q · ∇φ, 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 − Tr K d · I and the trace-part the following holds: (a) There exists an antisymmetric matrix J such that lim μ→∞ σ 2 f 0 (μ, J ) = 0, and there is an algorithmic way (see Algorithm 1) to compute an appropriate J in terms of K . (b) The trace-part is not effected by the perturbation, i.e. σ 2 f 1 (μ, J ) = σ 2 f 1 (0) for all μ ∈ R. Proof To prove the first claim, according to Theorem 5 it is sufficient to show that f 0 ∈ ker(J q · ∇) ⊥ = im(J q · ∇). Let us consider the function φ(q) = q · Aq, with A ∈ R d×d sym . It holds that J q · ∇φ = 2q · (J T Aq) = q · [A, J ]q. The task of finding an antisymmetric matrix J such that lim μ→∞ σ 2 f 0 (μ, J ) = 0 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 ∈ O(R d ) such that U K 0 U T has zero entries on the diagonal, and that U can be obtained in an algorithmic manner (see for example [29] or [22,Chap. 2,Sect. 2,Problem 3]) Assume thus that such a matrix U ∈ O(R d ) has been found and choose real numbers a 1 , . . . , a d ∈ R such that a i = a j if i = j. We now setĀ = diag(a 1 , . . . , a n ), and Observe that since U K 0 U T is symmetric,J is antisymmetric. A short calculation shows that [Ā,J ] = U K 0 U T . We can thus define A = U TĀ U and J = U TJ U to obtain [A, J ] = K 0 . Therefore, the J constructed in this way indeed satisfies (4.4). For the second claim, note that f 1 ∈ ker(J q · ∇), since J q · ∇ q · Tr K d q = 2 Tr K d q · J q = 0 due to 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 , . . . , a d ∈ R in its proof. However, it is asymptotically optimal: Corollary 2 In the setting of Lemma 8 the following holds: Proof The claim follows immediately since f 1 ∈ ker(J q · ∇) for arbitrary antisymmetric J as shown in (4.4), 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 ∈ R d×d sym , determine an optimal antisymmetric perturbation J as follows: . . d such that a i = a j for i = j and set Remark 12 In [14], the authors consider the task of finding optimal perturbations J for the nonreversible overdamped Langevin dynamics given in (15). In the Gaussian case this optimization problem turns out be equivalent to the one considered in this section. Indeed, equation (39) of [14] can be rephrased as f ∈ ker(J q · ∇) ⊥ . Therefore, Algorithm 1 and its generalization Algorithm 2 (described in Sect. 4.5) can be used without modifications to find optimal perturbations of overdamped Langevin dynamics.

Gaussians with Arbitrary Covariance and Preconditioning
In this section we extend the results of the preceding sections to the case when the target measure π is given by a Gaussian with arbitrary covariance, i.e. V (q) = 1 2 q · Sq with S ∈ R d×d sym symmetric and positive definite. The dynamics (8) then takes the form The key observation is now that the choices M = S and Γ = γ S together with the transformation q = S 1/2 q and p = S −1/2 p lead to the dynamics which is of the form (34) if J 1 and J 2 obey the condition S J 1 S = J 2 . Clearly the dynamics (66) is ergodic with respect to a Gaussian measure with unit covariance, in the following denoted by π. The connection between the asymptotic variances associated to (65) and (66) is as follows: For an observable f ∈ L 2 in which case it is given by Naturally, it is reasonable to choose m (i) in such a way that the exponential rate (λ (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 (69) quite stiff, requiring a very small timestep Δt in a discretisation scheme. The choice of c will therefore need to strike a balance between those two competing effects. The constraint (71) then implies Γ = 2cS. By a coordinate transformation, the preceding argument also applies if S, M and Γ are diagonal in the same basis, and of course M and Γ 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 Γ nondiagonally with respect to S (although only by a small margin). A clearer understanding of this is a topic of further investigation.

Numerical Scheme
In this section we introduce a splitting scheme for simulating the perturbed underdamped Langevin dynamics given by Eq. (8). In the unpertubed case, i.e. when J 1 = J 2 = 0, the B AO AB scheme (see [33] and references therein) has proven to be efficient for computing 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: where RK 4 (Δt, q 0 ) refers to fourth order Runge-Kutta integration of the ODĖ up until time Δt. We remark that the J 2 -perturbation is linear and can therefore be included in the O-part without much computational overhead. We emphasize the fact that many different splitting schemes could be investigated: although the BAOAB-scheme works well for unperturbed Langevin dynamics, it is not clear whether this remains true for the perturbed dynamics. Moreover, the perturbations introduced by J 1 and J 2 can be added in various places. Other discretisation schemes for the ODE (73) could be useful as well, for instance one could use a symplectic integrator, using the Hamiltonian structure of (73). However, since V as the Hamiltonian for (73) is not separable in general, such a symplectic integrator would have to be implicit. Note that (72c) and (72e) could be merged since (72e) commutes with (72d). In this paper, we content ourselves with the above scheme for our numerical experiments. Investigation of optimal numerical schemes for perturbed Langevin dynamics is an interesting problem for further research.

Remark 15
The aformentioned schemes lead to an error in the approximation for π( f ), since the invariant measure π is not preserved exactly by the numerical scheme. In practice, the B AO AB-scheme can therefore be accompanied by an accept-reject Metropolis step as in [40], leading to an unbiased estimate of π( f ), albeit with an inflated variance. In this case, after every rejection the momentum variable has to be flipped ( p → −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 → −J 1 and J 2 → −J 2 ) and using an appropriate (volume-preserving and time-reversible) integrator for the dynamics given by (73). Implementations of this idea are the subject of ongoing work. See [47] for a similar approach to nonreversible overdamped Langevin dynamics.

Diffusion Bridge Sampling
To numerically test our analytical results, we will apply the dynamics (8) to sample a measure on path space associated to a diffusion bridge. Specifically, consider the SDE dX s = −∇U (X s )ds + 2β −1 dW s , with X s ∈ R n , β > 0 and the potential U : R n → R obeying adequate growth and smoothness conditions (see [24], Sect. 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 π on L 2 ([0, s + ], 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 [9] and [45]). For a more detailed introduction (including applications) see [11] and for a rigorous theoretical treatment the papers [11,[24][25][26].
In the case U ≡ 0, it can be shown that the law of the conditioned process is given by a Gaussian measure π 0 with mean zero and precision operator S = − β 2 Δ on the Sobolev space H 1 ([0, s + ], R d ) equipped with appropriate boundary conditions. The general case can then be understood as a perturbation thereof: The measure π is absolutely continuous with respect to π 0 with Radon-Nikodym derivative where We will make the choice x − = x + = 0, which is possible without loss of generality as explained in [ in an equidistant way with stespize s j+1 − s j ≡ δ = 1 d+1 . Functions on this grid are determined by the values x(s 1 ) = x 1 , . . . , x(s n ) = x n , recalling that x(0) = x(1) = 0 by the Dirichlet boundary conditions. We discretise the functional Ψ as such that its gradient is given by We denote by A δ the discretised Dirichlet Laplacian on [0, 1] with stepsize δ. Following (74), the discretised target measure π has the form In the following we will consider the case n = 1 with potential U : R → R given by U (x) = 1 2 (x 2 − 1) 2 and set β = 1. To test our algorithm we adjust the parameters M, Γ , J 1 and J 2 according to the recommended choice in the Gaussian case, (27), where we take S = β 2 δ · A δ as the precision operator of the Gaussian target. We will consider the linear observable f 1 (x) = l · x with l = (1, . . . , 1) and the quadratic observable f 2 (x) = |x| 2 . In a first experiment we adjust the perturbations J 1 and J 2 to the observable f 2 according to Algorithm 2. The dynamics (8) is integrated using the splitting scheme introduced in Sect. 5.1 with a stepsize of Δt = 10 −4 over the time interval [0, T ] with T = 10 2 . Furthermore, we choose initial conditions q 0 = (1, . . . , 1), p 0 = (0, . . . , 0) and introduce a burn-in time T 0 = 1, i.e. we take the estimator to beπ( f ) ≈ 1 T −T 0 T T 0 f (q t )dt. We compute the variance of the above estimator from N = 500 realisations and compare the results for different choices of the friction coefficient γ and of the perturbation strength μ.
The numerical experiments show that the perturbed dynamics generally outperform the unperturbed dynamics independently of the choice of μ and γ , both for linear and quadratic observables. One notable exception is the behaviour of the linear observable for small friction γ = 10 −3 (see Fig. 4a), where the asymptotic variance initially increases for small perturbation strengths μ. This does not contradict our analytical results, since the small perturbation results Theorem 3 and Corollary 3 are only valid if γ > √ 2. We remark here that the condition γ > √ 2, while necessary for the theoretical results from Sect. 4.1, is not a very advisable choice in practice (at least in this experiment), since Figs. 4b and 5b clearly indicate that the optimal friction is around γ ≈ 10 −1 . Interestingly, the problem of choosing a suitable value for the friction coefficient coefficient γ becomes mitigated by the introduction of the perturbation: While the performance of the unperturbed sampler depends quite sensitively on γ , the asymptotic variance of the perturbed dynamics is a lot more stable with respect to variations of γ . A somewhat surprising phenomenon is that the standard deviation σ associated to the linear observable decays in the range γ ∈ [10, 100] for the unperturbed sampler (see Fig. 4b). We confirmed this behaviour by further numerical experiments and remark that as the target measureπ is fairly complicated, convexity of the function γ → σ should not be expected. In the regime of growing values of μ, the experiments confirm the results from Sect. 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 choice of J 1 [as long as J 2 is adjusted according to (27)]. This result is in alignment with Propostion 3 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.

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 [14] 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.