Romain Poncet
email: romain.poncet@cmap.polytechnique.fr

GENERALIZED AND HYBRID METROPOLIS-HASTINGS OVERDAMPED LANGEVIN ALGORITHMS

Keywords: Non-reversible diffusions, Langevin samplers, Markov Chain Monte Carlo, Metropolis-Hastings, variance reduction, lifting method, MCMC

It has been shown that the nonreversible overdamped Langevin dynamics enjoy better convergence properties in terms of spectral gap and asymptotic variance than the reversible one ([12, 13, 16, 25, 20, 21, 8]). In this article we propose a variance reduction method for the Metropolis-Hastings Adjusted Langevin Algorithm (MALA) that makes use of the good behaviour of the these nonreversible dynamics. It consists in constructing a nonreversible Markov chain (with respect to the target invariant measure) by using a Generalized Metropolis-Hastings adjustment on a lifted state space. We present two variations of this method and we discuss the importance of a well-chosen proposal distribution in terms of average rejection probability. We conclude with numerical experimentations to compare our algorithms with the MALA, and show variance reduction of several order of magnitude in some favourable toy cases.

Introduction

This article proposes a new class of MCMC algorithms whose objective is to compute expectations

π(f ) := E π (f ) = R d f (x)π(dx), (1.1)
for a given observable f , with respect to a probability measure π(dx) absolutely continuous, with respect to the Lebesgue measure, with density π(x) = e -U (x) . We suppose, as it is the case in many practical situations, that π is only known up to a multiplicative constant.

Many techniques have been developed to solve this problem. Deterministic quadratures can be very efficient at low dimension. Yet, in the high dimensional case, these methods tend to become inefficient, and MCMC methods can be used instead. The basic idea is to construct an ergodic Markov chain with respect to π, and to approximate π(f ) by the time average of this Markov chain. There are infinitely many ways to construct such a discrete time process. The general idea is to use an approximate time discretization of a time-continuous process known to be ergodic with respect to π. Generally, we cannot expect this discrete time process to be ergodic with respect to π. Thus one can use a Metropolis-Hastings acceptance-rejection step that ensures the detailed balance, and thus makes the chain reversible and ergodic with respect to π. In the case of the Euler-Maruyama discretization of the overdamped Langevin dynamics,

dX t = ∇ log π(X) + √ 2dW t , (1.2)
with (W t ) t≥0 a standard Brownian motion in R d , the method is called Metropolis Adjusted Langevin Algorithm (MALA, [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]).

Yet, it has been noticed in several contexts that departing from the reversibility can improve the performances of MCMC methods. This article aims to propose a generalization of the standard MALA that can be able construct nonreversible Markov chains that can outperform classical MALA.

1.1. Nonreversible dynamics. On the continous time setting, analysis have been carried out to compare the convergence properties of some time-continuous dynamics that are ergodic with respect to π [START_REF] Hwang | Accelerating gaussian diffusions[END_REF][START_REF]Accelerating diffusions[END_REF][START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Wu | Attaining the optimal gaussian diffusion acceleration[END_REF][START_REF] Rey | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF][START_REF] Rey | Variance reduction for irreversible langevin samplers and diffusion on graphs[END_REF][START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF], based on two kinds of optimality criterion: the speed of convergence toward equilibrium, measured in terms of spectral gap in L 2 (π), and asymptotic variance for the time averages. Obviously from a computational point of view an increase of the spectral gap enables to reduce the burnin, and a reduction of the asymptotic variance leads to a decrease of the computational complexity of the corresponding MCMC method. These analysis compare, for different vector fields γ, the overdamped Langevin dynamics given by,

dX t = ∇ log π(X t )dt + γ(X t )dt + √ 2dW t . (1.3)
Under the condition of non explosion and that the vector field γ is taken such that ∇ • (γπ) = 0, this dynamic is ergodic with respect to π. Such vector fields can be constructed easily: for any skew-symmetric matrix J, the vector field γ defined by γ(x) = J∇ log π(x) satisfies this divergence-free equation. Moreover, under this hypothesis, the following equation,

dX t = γ(X t )dt, (1.4)
conserves the energy U , which justifies the Hamiltonian denomination of the term γ. Moreover, this dynamic is time reversible if and only if γ = 0, and in this case, the detailed balance is satisfied (that is to say that the generator of the diffusion (1.3) is self-adjoint in L 2 (π)). It is well known that among all vector fields γ satisfying the non explosion condition and such that ∇ • (γπ) = 0, the dynamics given by (1.2) in the reversible case (γ = 0) has the worse rate of convergence in terms of spectral gap in L 2 (π) [START_REF] Hwang | Accelerating gaussian diffusions[END_REF][START_REF]Accelerating diffusions[END_REF][START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF]. Recent work has been done to construct divergence free (with respect to π) perturbations of the drift that achieve optimal convergence properties in the Gaussian case [START_REF] Lelièvre | Optimal non-reversible linear drift for the convergence to equilibrium of a diffusion[END_REF][START_REF] Wu | Attaining the optimal gaussian diffusion acceleration[END_REF]. Recent works also show that breaking the non reversibility with such divergence free perturbations on the drift also leads to improvement on the asymptotic variance. It is shown in [START_REF] Rey | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF] that the asymptotic variance decreases under the addition of the irreversible drift. Moreover, is has been shown in [START_REF] Rey | Variance reduction for irreversible langevin samplers and diffusion on graphs[END_REF] that the asymptotic variance is monotonically decreasing with respect to the growth of the drift, and the limiting behavior for infinitely strong drifts is characterized. More recently, in [START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF] the authors investigate the dependence of the asymptotic variance on the strength of the nonreversible perturbation.

On the discrete time setting, classical methods that depart from reversible sampling consist in hybrid (Hamiltonian) MCMC [START_REF] Simon Duane | Hybrid monte carlo[END_REF][START_REF] Lelièvre | Free energy computations : a mathematical perspective[END_REF] and generalized hybrid MCMC methods [START_REF] Kennedy | Cost of the generalised hybrid monte carlo algorithm for free field theory[END_REF]. In the former method, the drift direction is chosen isotropically at each time step and long time Hamiltonian integration is then carried out in this direction. The latter can be seen as a generalization of the former that brings some inertia in the direction of the Hamiltonian dynamics. Another class of nonreversible samplers is composed by lifting methods. They are designed in the discrete state space case, to construct a Markov chain that satisfies some skew detailed balance [START_REF] Diaconis | Analysis of a nonreversible markov chain sampler[END_REF][START_REF] Chen | Lifting markov chains to speed up mixing[END_REF][START_REF] Hukushima | An irreversible markov-chain monte carlo method with skew detailed balance conditions[END_REF][START_REF] Vucelja | Lifting-a nonreversible markov chain monte carlo algorithm[END_REF]. They consist in increasing the state space to take into account a privileged drift direction that is explored more efficiently. More recently, Bierkens proposed an extension of the classical Metropolis-Hastings algorithm to generate unbiased nonreversible Markov chain [START_REF] Bierkens | Non-reversible metropolis-hastings[END_REF]. This is achieved by modifying the acceptance probability to depart from detailed balance. Eventually, a recent and quite different approach has been proposed in [START_REF] Bierkens | The zig-zag process and super-efficient sampling for bayesian analysis of big data[END_REF] in the big data settings to circumvent the poor scalability of standard MCMC methods. The authors construct a continuous time piecewise deterministic Markov process. It is a constant velocity model where the velocity direction switches at random times with a rate depending on the target distribution. 1.2. Outline. We propose in this article a bias-free algorithm similar to MALA that aims to exploit the asymptotic variance reduction of the nonreversible time-continuous process. The idea is to construct a Markov chain with invariant measure π, by discretizing an equation of the form (1.3), instead of equation (1.2), enhanced with an acceptance-rejection step. The main difficulty consists in unbiasing the unadjusted chain. Indeed, it is not worth considering the use of a standard Metropolis-Hastings acceptance probability since it is designed to impose detailed balance with respect to the target distribution π, and thus to define a reversible Markov chain with respect to π. It would lead to a poor average acceptance ratio. An elegant way would be to construct an adequate acceptance probability with respect to this proposal, to ensure a high average acceptance ratio. In the setting of [START_REF] Bierkens | Non-reversible metropolis-hastings[END_REF], it would consist in finding a good vorticity kernel. Yet, we are not able to exactly do this. Instead, we propose a class of lifted algorithms that rely on these unadjusted chains. More precisely the first algorithm is a generalized Metropolis-Hastings algorithm in an enhanced state space, and the second one can be seen as the analogous of the generalized hybrid Monte Carlo method for the overdamped Langevin equation.

In section 2 we present the first algorithm (generalized MALA). We discuss in 2.1 how its performances are closely related to the choice of the transition kernel of the unadjusted chain. In section 2.2 we prove geometric convergence of the Markov chain constructed with this algorithm under some hypotheses, that ensure the existence of a central limit theorem. Then, in section 3 we propose a modification of this algorithm (the generalized hybrid MALA). In section 4 we present numerical comparisons of these algorithms with respect to classical MALA, which is followed by concluding remarks.

Generalized MALA

In this section, we construct a nonreversible Markov chain, ergodic with respect to a target distribution π known up to a normalizing constant. The algorithm is similar to MALA in the sense that it constructs a Markov chain from the discretization of an overdamped Langevin dynamic, augmented with an acceptance-rejection step that makes it ergodic with respect to π. The difference is that we construct a Markov chain on the discretization of a nonreversible Langevin equation to try to benefit from the smaller asymptotic variance of this kind of Markov processes, than the reversible ones. The main issue is then to choose a right acceptance probability that preserves the good ergodic properties of the underlying Markov process.

To state our algorithm, we slightly modify Equation (1.3). For ξ ∈ R, we consider the diffusions,

dX t = ∇ log π(X t )dt + ξγ(X t )dt + √ 2dW t , (2.1)
with divergence-free condition ∇ • (γπ) = 0. This way, ξ specifies the direction and the intensity of the nonreversibility. We denote now by Q ξ a proposal kernel that correspond to some discretization of the diffusions (2.1) with parameter ξ. We propose the following algorithm that we call Generalized MALA (GMALA), Algorithm 2.1 (Generalized MALA). Let h > 0, (x 0 , ξ 0 ) ∈ R d × R be an initial point and an initial direction. Iterate on n ≥ 0.

(1) Sample y n+1 according to Q ξ n (x n , dy).

(2) Accept the move with probability

A ξ n (x n , y n+1 ) = 1 ∧ π(y n+1 )Q -ξ n (y n+1 , x n ) π(x n )Q ξ n (x n , y n+1 ) . (2.2)
and set (x n+1 , ξ n+1 ) = (y n+1 , ξ n ); otherwise set (x n+1 , ξ n+1 ) = (x n , -ξ n ).

The important part of this algorithm is that the direction ξ n of the Hamiltoninan exploration must be inverted at each rejection to ensure its unbiasedness. A good choice of Q ξ is given in the next section.

Unbiasedness is obvious since this algorithm is actually built as a Generalized Metropolis-Hastings algorithm on the increased state space E = R d × {-ξ 0 , ξ 0 }. To simplify notations, we denote by x ξ all element (x, ξ) ∈ E. We set S the involutive transformation defined for all element x ξ ∈ E by S(x ξ ) = x -ξ . We extend the definition of π on E by π(x ξ ) = π(x) 2 for x ξ ∈ E. Obviously π is unchanged by S. Then, the algorithm constructs a Markov chain with transition kernel density P given by,

P (x ξ , y η ) = Q ξ (x, y)A ξ (x, y)1 ξ (η) + δ S(x ξ ) (y η ) 1 - R d Q ξ (x, z)A ξ (x, z)dz ,
where 1 xi denotes the characteristic function and δ S(x ξ ) a Dirac delta function, that satisfies the following skew detailed balance, ∀x ξ , y η ∈ E, π(x ξ )P (x ξ , y η ) = π(y η )P (S(y η ), S(x ξ )).

(2.3)

Heuristically, we hope that the discretization of the time-continuous process (1.3) specified by Q ξ benefits from the same good behavior in terms of asymptotic variance reduction. Since the Markov chain is constructed as parts of the discretization of the time-continuous dynamics between the rejections, then the closer to one is the average acceptance probability, the longer are these parts, and the more we can hope to benefit from this good behavior. Thus, to give a hint about the relevance of this algorithm, we compute in section 2.1 these average acceptance probabilities, and we show that they are of the same order of those for MALA for some well-chosen discretization. Moreover, we numerically show in section 4 that it can outperform MALA by several order of magnitude in terms of asymptotic variance.

Yet, before showing these results, we present a heuristic justification about this algorithm. In the reversible case, the time-continous equation satisfies the detailed balance with respect to π, ∀x, y ∈ R d , ∀h > 0, π(x)P h (x, y) = π(y)P h (y, x), where P h (x, y) denotes the density of the transition kernel to go from x to y after a time h for the dynamics (1.2). Moreover, MALA imposes that the Markov chain also satisfies the detailed balance with respect to π. In the nonreversible case, the time continuous process (2.1) satisfies a skew detailed balance as stated by the following lemma, Lemma 2.1. For all x, y ∈ R d , for all ξ ∈ R and for all h > 0, the following relation holds,

π(x)P ξ h (x, y) = π(y)P -ξ h (y, x), (2.4)
where P ξ h (x, dy) is the transition probability measure of the h-skeleton of the process (X ξ t ) t≥0 , solution of Equation (2.1).

The analogous of the reversible case would be to construct a π-invariant Markov chain that satisfies the same kind of skew detailed balance. Because the skew detailed balance (2.4) can be seen as a detailed balance on the enhanced state space E up to the transformation S,

∀x, y ∈ R d , ∀ξ ∈ R, ∀h > 0 π(x)P h (x ξ , y ξ ) = π(y)P h (S(y ξ ), S(x ξ )), (2.5)
where P h (x ξ , y ξ ) = P ξ h (x, y), the Generalized Metropolis-Hastings method given by algorithm 2.1 is actually the classical way to construct such a Markov chain. Nevertheless, the main difference between the time-continous and the discrete time dynamics is the fact that the latter requires some direction switching of the nonreversible component of the dynamics. As stated previously, the idea of lifting the state space to construct a nonreversible chain has been used in the discrete state space setting. Yet, the idea is quite different. With classical lifting methods, the goal is to switch between several directions with well-chosen probabilities, to quickly explore the state space in all of these directions. In our case, we do not aim to switch between several nonreversible directions. Yet, we are forced to do so at each rejection, and we have no choice but to reverse the current nonreversible directions. That is to say that we control neither the probability of switching nor the direction.

Proof of Lemma 2.1. We denote by L ξ the generator of diffusion (2.1). One can show that,

L ξ = S + ξA,
where S = ∇ log π(x)∇•, and A = γ •∇. Then, one can show that S and A are respectively the symmetric and the antisymmetric parts of L 1 , with respect to L 2 (π). Moreover, one can show that for all ξ ∈ R, L ξ and L -ξ are adjoint from one another. This amounts to show that D((L -ξ ) * ) ⊂ D(L ξ ), which follows from the injectivity of the operator (L -ξ ) * -iI, and the surjectivity of the operator L ξ -iI. It follows that the semigroups e L ξ t and e L -ξ t are adjoint from one another (Corollary 10.6 [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]). Then, for all bounded functions f and g,

(R d ) 2 f (x)g(y)P ξ h (x, dy)π(dx) = R d f (x)e L ξ h g(x)π(dx), = R d e (L ξ ) * h f (x)g(x)π(dx), = R d e L -ξ h f (x)g(x)π(dx), = (R d ) 2 f (y)g(x)P -ξ h (x, dy)π(dx).
2.1. Choice of the proposition kernel. The method presented above can be tuned with the choice of the proposition kernel Q ξ . It should be chosen such that it approximates the law of the transition probability of the h-skeleton of the process (X ξ t ) t≥0 solution of equation (2.1). The basic idea is to define Q ξ as the density of the law of an approximate discretization of equation (2.1). Doing so, we can hope that the skew-detailed balance for the unadjusted chain would be almost satisfied (since it is satisfied for the time-continous dynamics), that is to say we can hope to benefit from an acceptance ratio close to one. Moreover, we recall that the choice of the proposal kernel is the only degree of freedom that enables to tune the rate of the direction switching. Formally, the closer the proposal kernel is from the transition kernel of the h-skeleton of the continuous process, the higher the acceptance ratio is (and thus the less frequent is the direction switching).

The basic idea would be to propose according to the Maruyama-Euler approximation of equation (2.1), as it is done for MALA. We denote by Q ξ 1 this kernel. That is to say Q ξ 1 (x, dy) is given by the law of y, solution of

y = x -h∇U (x) -hξγ(x) + √ 2hχ, (2.6)
where χ a standard normal deviate. Then, Q ξ 1 is given by

Q ξ 1 (x, dy) = 1 (4πh) d/2 exp 1 4h ||y -(x -h∇U (x) -hξγ(x))|| 2 dy. (2.7)
Sadly, even though the simplicity of this proposal is appealing, this proposition kernel leads to an average rejection rate of order h when ξ = 0 and γ is non-linear as stated in Proposition 2. It is significantly worse than MALA that enjoys an average rejection rate of order h 3/2 (see [START_REF] Bou-Rabee | Pathwise accuracy and ergodicity of metropolized integrators for SDEs[END_REF]). As shown later by the numerical simulations, this bad rejection rate is not a pure theoretical problem: it forbids to use large discretization steps h, and thus the method only generates highly correlated samples.

To overcome this problem, we propose to implicit the resolution of the nonreversible term in equation (2.1) with a centered point discretization. More precisely, we propose a move y ξ from x ξ , such that y would be distributed according to the solution of

Φ hξ x (y) = x -h∇U (x) + √ 2hχ, (2.8)
where χ a standard normal deviate, and Φ hξ x is the function defined by,

Φ hξ x : R d → R d , y → y + hξγ x + y 2 .
(2.9)

Existence of such a proposal kernel, denoted by Q ξ 2 , is ensured under the hypothesis that U is twice differentiable with uniformly bounded second derivative.

Assumption 1. We suppose that U is twice differentiable with uniformly bounded second derivative. Proposition 1. Under Assumption 1, there exists h 0 > 0, such that for all h < h 0 , for all x ξ ∈ E, for all χ ∈ R d , there exists a unique y ∈ R d , solution of equation (2.8), and if χ denotes a standard normal deviate, the function Y ξ defined almost surely as the solution of equation (2.8) is a well-defined random variable, with law Q ξ 2 (x, dy), given by,

Q ξ 2 (x, dy) = 1 (4πh) d/2 | Jac φ ξh x (y)| exp 1 4h ||φ hξ x (y) -(x -h∇ log π(x))|| 2 dy. (2.10) Moreover x -y = √ 2hχ + O h ||∇U (x)|| + h 3/2 ||χ|| ,
and there exists C > 0, independent of x, χ, and h, such that,

||∇U (y)|| ≤ (1 + Ch) ||∇U (x)|| + C √ 2h ||χ|| .
Proof of Proposition 1. The first point is a corollary of Lemma 2.2 below. The second point can be proven by making use of the fact that ∇U is Lipschitz.

The following lemma is used to prove Proposition 1. We can note that the h 0 given by this lemma depends on the Lipschitz coefficient of ∇U . Lemma 2.2. Under Assumption 1, there exists h 0 > 0, such that for all h < h 0 , for all x ξ ∈ E, the function

Φ hξ x : R d → R d , defined by (2.9) is a C 1 -diffeomorphism. Proof of Lemma 2.2. The proof that Φ hξ
x is bijective relies on Picard fixed point theorem. For any fixed

z ∈ R d , we define Ψ on R d by Ψ(y) = z -hξJ∇U x + y 2 .
Then, since ∇U is supposed to be Lipschitz, then for small enough h this application is a contraction mapping from R d to R d . The fact that Φ hξ x is everywhere differentiable is clear, and the fact that its inverse is everywhere differentiable as well comes from the fact that the determinant of the Jacobian of Φ hξ

x is strictly positive for sufficiently small h.

Sampling from the this proposal kernel can be done by a fixed point method when no analytical solution is available since the proof of Lemma 2.2 uses a Picard fixed point argument. This transition kernel leads to a rejection rate of order h 3/2 (see Proposition 2), which is en improvement from Q ξ 1 . The main advantage of this kernel is that the computation of | Jac φ ξh x (y)| can be avoided in the case where γ is defined by γ(x) = J∇ log π(x), with J a skew symmetric matrix. Indeed, only the ratio

| Jac φ ξh x (y)|/| Jac φ -ξh y (x)
| is required to compute the acceptance probability A(x, y), and this ratio is equal to 1 if γ is of gradient type, as stated by the following lemma. Lemma 2.3. For all x, y ∈ R d , and for all h > 0,

| Jac φ ξh x (y)| = | Jac φ -ξh y (x)|
Proof of lemma 2.3. This result uses the more general fact that for any skew-symmetric matrix A and any symmetric matrix S, the matrices Id + AS and Id -AS have same determinant. This statement is equivalent to say that χ AS = χ -AS , where we denote by χ M the characteristic polynomial of any square matrix M . This last statement is true since for any square matrices A and S (of the same size) χ AS = χ SA . Then using the transposition, χ AS = χ A t S t . Eventually using the fact that A t = -A and S t = S, it comes χ AS = χ -AS . The result follows from the fact that the matrix Jac φ ξh x (y) is of the form Id + ξAS.

We denote by α 1 h,ξ and α 2 h,ξ respectively the acceptance probability for proposal kernels Q ξ 1 and Q ξ 2 (defined respectively by (2.7) and (2.10)). The following proposition provides an upper bound on the moments of the rejection probability.

Proposition 2. Suppose that U is three times differentiable with bounded second and third derivatives.

Then for all l ≥ 1, there exists C(l) > 0 and h 0 > 0 such that for all positive h < h 0 , and for all x ∈ R d ,

E (1 -α 1 h,ξ (x, Y 1 h,ξ )) 2l ≤ C(l)(1 + ||∇U (x)|| 4l )h 2l E (1 -α 2 h,ξ (x, Y 2 h,ξ )) 2l ≤ C(l)(1 + ||∇U (x)|| 4l )h 3l
Proof of Proposition 2. To deal with both results at once, we define for all x ∈ R d , for all ξ ∈ {-1, 1}, and for all θ ∈ {0, 1}, the random variable Y x ξ ,θ that satisfies the following implicit equation almost surely,

Y x ξ ,θ = x -h∇U (x) -hξJ θ∇U x + Y x ξ ,θ 2 + (1 -θ)∇U (x) + √ 2hχ,
where χ is a standard normal deviate in R d . Well-posedness of Y x ξ ,θ is given by Proposition 1. We set R ξ θ the associated proposal kernel. We get

R ξ 0 = Q ξ 1 and R ξ 1 = Q ξ 2 .
We define the Metropolis-Hastings ratio r(x ξ , y ξ ) for proposing y ξ from x ξ with kernel R ξ θ by,

r(x ξ , y ξ ) = π(y -ξ )R θ (y -ξ , x -ξ ) π(x ξ )R θ (x ξ , y ξ ) ,
where we set χ ∈ R d such that,

y = x -h∇U (x) -hξJ θ∇U x + y 2 + (1 -θ)∇U (x) + √ 2hχ,
Then, a straightforward computation gives, log(r(x ξ , y ξ )) =U (x) -U (y) + y -x, ∇U (x)

+ 1 2 y -x, ∇U (y) -∇U (x) - ξ 2 (1 -θ) y -x, J(∇U (y) -∇U (x)) + h 4 ||∇U (x)|| 2 -||∇U (y)|| 2 + h 2 ξθ J∇U x + y 2 , ∇U (x) + ∇U (y) + h 2 ξ 2 θ(1 -θ) J∇U x + y 2 , J(∇U (x) -∇U (y)) + h 4 (1 -θ) 2 ξ 2 ||J∇U (x)|| 2 -||J∇U (y)|| 2 .
A Taylor expansion of these terms, making use of Proposition 1 and noticing that θ(1

-θ) = 0 for θ ∈ {0, 1} leads to log(r(x ξ , y ξ )) = -ξh(1 -θ) χ, JD 2 U (x) • χ + O(h 3/2 (||∇U (x)|| 2 + ||χ|| 2 + ||χ|| 3 )) =O(h(1 -θ) ||χ|| 2 + h 3/2 (||∇U (x)|| 2 + ||χ|| 2 + ||χ|| 3 )).
Moreover,

(1 -1 ∧ r(x ξ , y ξ )) 2l =O(log(r(x ξ , y ξ )) 2l ) =O(h 2l (1 -θ) 2l ||χ|| 4l + h 3l (||∇U (x)|| 4l + ||χ|| 4l + ||χ|| 6l )),
and thus there exists C(l) > 0 such that,

E (1 -1 ∧ r(x ξ , Y x ξ ,θ )) 2l ≤ C(l)(1 + ||∇U (x)|| 4l )h 3l + C(l)h 2l (1 -θ) 2l
The previous method is quite efficient since no computation of the Hessian of log π is required. Nevertheless, global Lipschitzness of ∇U is required to justify the method. We propose a last kernel, denoted by Q ξ 3 , that does not require this hypothesis but still require γ to be of gradient type and the computation of the Hessian of U . More precisely, Q ξ 3 (x, dy) is the law of y, solution of

Id + hξ 2 J Hess(U )(x) (y -x) = -h(Id + ξJ)∇U (x) + √ 2hχ, (2.11)
where χ a standard normal deviate. Then, Q ξ 3 is given by

Q ξ 3 (x, dy) = 1 (4πh) d/2 det M ξ (x) exp 1 4h ||M ξ (x)(y -x) + h(Id + ξJ)∇U (x)|| 2 dy, (2.12)
with M ξ (x) = Id + hξ 2 J Hess U (x). This transition kernel offers also a rejection rate a 3 h,ξ of order h 3/2 . Proposition 3. Suppose that U is three times differentiable with bounded second and third derivatives. Then for all l ≥ 1, there exists C(l) > 0 and h 0 > 0 such that for all positive h < h 0 , and for all x ∈ R d ,

E (1 -α 3 h,ξ (x, Y 2 h,ξ )) 2l ≤ C(l)(1 + ||∇U (x)|| 4l )h 3l
Proof of Proposition 3. The proof involves the same arguments as in Proposition 2 and is left to the reader.

Even though Kernel Q ξ 2 enables to circumvent the bad acceptance ratio of the explicit proposal given by Q ξ 1 , it raises some difficulties. First, it requires the potential U to be globally Lipschitz, which is quite restrictive. To use the GMALA method in the non globally Lipschitz case, one can resort to importance sampling (to get back to a globally Lipschitz setting), or use a globally Lipschitz truncation of the potential to build the proposition kernel. The same idea is used in MALTA [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]. This last trick might lead to high average rejection ratio in the truncated regions of the potential. Moreover, the computation of the proposed move requires the resolution of a nonlinearly implicit equation, that can be solved by a fixed point method, which is more costly than the Euler-Maruyama discretization used by MALA. Specific methods of preconditioning should be considered to accelerate the convergence of the fixed point.

It would be interesting to propose a proposition kernel based on an explicit scheme, that would achieve a better acceptance ratio than Q ξ 1 , with a lower computational cost than Q ξ 2 . The question of decreasing the Metropolis-Hastings rejection rate has been recently studied in [START_REF] Fathi | Improving dynamical properties of stabilized discretizations of overdamped langevin dynamics[END_REF]. The authors propose a proposition kernel constructed from an explicit scheme, which is a correction (at order h 3/2 ) of the standard Euler-Maruyama proposal. Nevertheless, in our case this approach would not enable us to construct a proposition kernel based on an explicit scheme with a high average acceptance ratio. Work has also been done in this direction with the Metropolis-Hastings algorithm with delayed rejection, proposed in [START_REF] Tierney | Some adaptive monte carlo methods for bayesian inference[END_REF], and more recently [START_REF] Banterle | Accelerating metropolis-hastings algorithms by delayed acceptance[END_REF]. Yet, it is unclear whether this approach could be use efficiently in our case.

Convergence of GMALA.

In this section, we only treat the case in which GMALA is used with Q ξ 2 . To obtain a central limit theorem for the Markov chain built with GMALA, it is convenient to prove the geometric convergence in total variation norm toward the target measure π. Because we require ∇U to be globally Lipschitz to ensure well-posedness of GMALA with transition kernel Q ξ 2 , MALA would be likely to benefit from geometric convergence in this case. Indeed, it is well-known that MALA can be geometrically ergodic when the tails of π are heavier than Gaussian, under suitable hypotheses (see Theorem 4.1 [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]). For strictly lighter tails, we cannot hope for geometric convergence (see Theorem 4.2 [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]). This part is devoted to showing that under some hypotheses, GMALA can benefit from geometric convergence.

Classically for MALA, the convergence follows on from the aperiodicity and the irreducibility of the chain, since by construction the chain is positive with invariant measure π. Under these conditions, the geometric convergence can be proven by exhibiting a suitable Foster-Lyapunov function V such that,

lim |x|→∞ P V (x) V (x) < 1. (2.13)
In our case, the situation is slightly different. The Markov chain built with GMALA is still aperiodic and phi-irreducible. This is a simple consequence of the surjectivity of the application Φ hξ x in Lemma 2.2. To be able to prove a drift condition such as (2.13), it is usually required to be able to show that the acceptance probability for a proposed move starting from x ξ does not vanish in expectation for large x, which is not always true in our setting. More precisely, we can only say that the maximum of the two average acceptance probabilities for the proposed moves starting from x ξ and x -ξ does not vanishes for large x. Then one strategy could be to choose a Foster-Lyapunov function V that decreases in expectation when ξ is changed to -ξ after a rejection. An other strategy, that we present in the following and that seems to be more natural and more easily generalizable to potential U that does not satisfy the specific hypotheses we use in this section. We show that the odd and the even subsequences of the Markov chain converge geometrically quickly to the target measure π by showing a drift condition on P 2 : the transition probability kernel of the two-steps Markov chain defined by P 2 (x, dy) = E P (x, dz)P (z, dy), ∀x ∈ E, under the following assumption.

Assumption 2. We suppose that U is three times differentiable, and that,

(1) eigenvalues of D 2 U (x) are uniformly upper bounded and lower bounded away from 0, for x outside of a ball centered in 0, (2) the product D 3 U (x) ||∇U (x)|| is bounded for x outside of a ball centered in 0. Proposition 4. Suppose assumption 2. There exists s > 0 and h 0 > 0 such that for all h ≤ h 0 , for ξ ∈ {-1, 1},

lim ||x||→+∞ P 2 V s (x ξ ) V s (x ξ ) = 0, (2.14)
where the Foster-Lyapunov function V is defined by V s (x ξ ) = exp(sU (x)).

Remark 2.1. In order to prove the drift condition (2.13), one can use the Foster-Lyapunov V defined by,

V s (x ξ ) = exp sU (x) + s ξh 2 ∇U (x), D 2 U (x)J∇U (x) |ξh 2 ∇U (x), D 2 U (x)J∇U (x) | ,
for s small enough. Yet, this proof is left to the reader, but uses the arguments developed in the following.

Assumption 2 is not meant to be sharp. We are not striving for optimality here, but instead we aim to propose a simple criterion which is likely to be satisfied in smooth cases.

The first step to prove equation (2.14), is to show that the proposed move decreases the Lyapunov function V s in expectation.

Lemma 2.4. Under assumption 1, there exists s 0 > 0 such that for all s < s 0 , there exists C 1 , C 2 > 0 such that for all h < h 0 given by Proposition 1, and for all x ∈ R d large enough (depending on h),

E [V s (Y h,ξ )] ≤ V s (x)e -s 3h(1-C 1 h) 4 ||∇U (x)|| 2 +sC2 ,
where Y h,ξ is defined in Proposition 1. Thus, for h small enough,

E [V s (Y h,ξ )] ≤ V s (x)e -s h 2 ||∇U (x)|| 2 +sC2 ,
and for x large enough,

E [V s (Y h,ξ )] ≤ V s (x)e -s h 4 ||∇U (x)|| 2 .
Proof of Lemma 2.4. For all y ∈ R d , we set χ ∈ R d such that,

y = x -h∇U (x) -hξJ θ∇U x + y 2 + (1 -θ)∇U (x) + √ 2hχ,
A Taylor expansion of y in h yields,

U (y) = U (x) + ∇U (x) • (y -x) + O(||y -x|| 2 ) = U (x) -h ||∇U (x)|| 2 + √ 2h∇U (x) • χ + hξ∇U (x) • J∇U ( x + y 2 ) + O(||x -y|| 2 ).
Noticing that

hξ∇U (x) • J∇U ( x + y 2 ) = hξ∇U (x) • J(∇U ( x + y 2 ) -∇U (x)),
the Cauchy-Schwarz and the triangular inequality yield,

hξ∇U (x) • J∇U ( x + y 2 ) ≤ Ch ||∇U (x)|| ||x -y|| .
Thus, by Young inequality and Proposition 1, for h ≤ 1,

U (y) ≤ U (x) - 3h(1 -Ch) 4 ||∇U (x)|| 2 + O(||χ|| 2 ),
The conclusion holds for small enough s such that e sO(||G|| 2 ) is integrable, where G is a standard normal deviate in R d .

Classically ( [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]), proofs of geometric convergence of MALA require to show that the average acceptance ratio of the proposed move from x, does not vanishes when x is large. Namely that there exists ε > 0 such that I(x) = {y : α(x, y) ≤ 0} asymptotically has q-measure 0, where α is the acceptance probability. Our case is sightly different since it is possible that the average acceptance ratio of the proposed move from x ξ ∈ E vanishes when ||x|| → +∞, which leads to a rejection and to the switching ξ ← -ξ with probability close to one. Yet the average acceptance probability of the next proposed move from x -ξ is close to one. This behavior is described by the following Lemma. Lemma 2.5. Under Assumption 2, there exists h 0 > 0 such that for all h < h 0 and for all ε > 0, there exists C(h, ε) > 0 such that for all

x ξ ∈ E such that ||x|| ≥ C(h, ε), ξ ∇U (x), D 2 U (x) • J∇U (x) ≥ 0 =⇒ P α 2 h,ξ (x ξ , y ξ ) = 1 ≥ 1 -ε,
where y is defined by Equation (2.8).

Proof of Lemma 2.5. We recall that α 2 h,ξ is defined for all x ξ , y ξ ∈ E by α 2 h,ξ (x ξ , y ξ ) = 1 ∧ e r(x ξ ,y ξ ) , with,

r(x ξ , y ξ ) = U (x) -U (y) + √ 2hχ, ∇U (x) + 1 2 √ 2hχ, ∇U (y) -∇U (x) -h ||∇U (x)|| 2 -h ∇U (y) -∇U (x), ∇U (x) - h 4 ||∇U (x) -∇U (y)|| 2 .
The proof is done by computing a Taylor expansion in h of this quantity to evaluate its sign in the asymptotic case ||x|| → +∞. We denote by • the matrix vector product, and by : and . . . respectively the double and triple dot products.

U (x) -U (y) + √ 2hχ, ∇U (x) -h ||∇U (x)|| 2 = hξ∇U (x) • J ∇U x + y 2 -∇U (x) - 1 2 D 2 U (x) : (y -x) 2 + O(||D 3 U (x) . . .(y -x) 3 ||) = h 2 ξ∇U (x) • J D 2 U (x) • (y -x) + O(D 3 U (x) : (y -x) 2 ) - 1 2 D 2 U (x) : (y -x) 2 + O(||D 3 U (x) . . .(y -x) 3 ||) 1 2 √ 2hχ, ∇U (y) -∇U (x) = 1 2 √ 2hχ, D 2 U (x) • (y -x) + O(D 3 U (x) : (y -x) 2 ) , -h ∇U (y) -∇U (x), ∇U (x) = -h D 2 U (x) • (y -x) + O(D 3 U (x) : (y -x) 2 ), ∇U (x) - h 4 ||∇U (x) -∇U (y)|| 2 = hO(||x -y|| 2 ).
Collecting these terms, and using Proposition 1 and Assumption 2, it comes for h ≤ 1

r(x ξ , y ξ ) = - h 2 ∇U (x), D 2 U (x) • (y -x) + O(h 3 ||∇U (x)|| 2 ) + O(||χ|| 2 + ||χ|| 3 ),
where the big O are independent of h. and eventually,

r(x ξ , y ξ ) = h 2 2 ∇U (x), D 2 U (x) • ∇U (x) + ξ ∇U (x), D 2 U (x) • J∇U (x) + O(h 3 ||∇U (x)|| 2 ) + O(||χ|| 2 + ||χ|| 3 ),
Then, using the definite positivity of D 2 U (x) and the fact that ||∇U (x)|| is coercive, for all h small enough and for all ε > 0, there exists C(h, ε) > 0 such that for all x ∈ R d such that ||x|| ≥ C(h, ε),

ξ ∇U (x), D 2 U (x) • J∇U (x) ≥ 0 =⇒ P (r(x ξ , y ξ ) ≥ 0) ≥ 1 -ε.
Proposition 4 can now be obtained as a consequence of Lemmas 2.5 and 2.4.

Proof of proposition 4. Because of the acceptance-rejection step, we get for all x ξ ∈ E,

P 2 V (x) = R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ ) R d q(y ξ , z ξ )α 2 h,ξ (y ξ , z ξ )V (z ξ )dzdy + R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ )V (y -ξ ) 1 - R d q(y ξ , z ξ )α 2 h,ξ (y ξ , z ξ )dz dy + 1 - R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ )dy R d q(x -ξ , y -ξ )α 2 h,ξ (x -ξ , y -ξ )V (y -ξ )dy + V (x -ξ ) 1 - R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ )dy 1 - R d q(x -ξ , y -ξ )α 2 h,ξ (x -ξ , y -ξ )dy .
Simply using α 2 h,ξ ≤ 1 leads to,

P 2 V (x) ≤ R d q(x ξ , y ξ ) R d q(y ξ , z ξ )V (z ξ )dzdy + R d q(x ξ , y ξ )V (y -ξ ) 1 - R d q(y ξ , z ξ )α 2 h,ξ (y ξ , z ξ )dz dy + 1 - R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ )dy R d q(x -ξ , y -ξ )V (y -ξ )dy + V (x -ξ ) 1 - R d q(x ξ , y ξ )α 2 h,ξ (x ξ , y ξ )dy 1 - R d q(x -ξ , y -ξ )α 2 h,ξ (x -ξ , y -ξ )dy .
And eventually Lemmas 2.5 and 2.4 gives,

P 2 V (x) V (x) ≤ e -s h 2 ||∇U (x)|| 2 +2sC2 + e -s h 2 ||∇U (x)|| 2 +sC2 + (1 -E α 2 h,ξ (x ξ , Y h,ξ ) )e -s h 2 ||∇U (x)|| 2 +sC2 + (1 -E α 2 h,ξ (x ξ , Y h,ξ ) )(1 -E α 2 h,ξ (x -ξ , Y h,-ξ ) ).
And thus,

lim ||x||→+∞ P 2 V (x) V (x) = 0 3. Generalized hybrid MALA
We propose in this part a second algorithm. It is based on a splitting method by solving by turns the reversible and the non reversible parts of equation (2.1) with measure preserving schemes. The idea is then to integrate the pure reversible equation dx = -∇U (x)dt + dW t .

(3.1) by using MALA, and to integrate the pure Hamiltonian equation

dx = -ξJ∇U (x)dt. (3.2)
by an hybrid Monte-Carlo method based on a suitable Hamiltonian integrator. This method presents some theoretical advantages on Generalized MALA. Before presenting them, we begin with explaining the algorithm.

We denote by Ψ ξ t the flow of the Hamiltonian equation (3.2) on the time interval [0, t], and by Φ ξ h a numerical integrator for (3.2) on a time step h. We precise later necessary conditions on this integrator that ensure unbiasedness of the algorithm. MALA is used to sample x n+1/2 from x n , with time-step h. (2) Integration of the non reversible part (3.2):

(a) Compute

x n+1 = Φ ξ h (x n+1/2 ). (b) Set x n+1 = x n+1 with probability β h,ξ (x n+1/2 ) = min 1, exp(U (x n+1/2 ) -U ( x n+1 )) min 1, exp(U (x n+1/2 ) -U (Φ ξ h (x n+1/2 ))) .
Otherwise set x n+1 = x n+1/2 and ξ ← -ξ.

Similarly to the Hybrid Monte-Carlo algorithm, the first step enables to explore the state space across the iso-potential lines, whereas the second step enables to explore it along the iso-potential lines.

To ensure unbiasedness of the second step, the integrator Φ ξ h must satisfy the following properties,

Φ ξ h = (Φ -ξ h ) -1 , (3.3)
det Jac Φ ξ h = 1. (

These properties are classical for hybrid Monte-Carlo methods (see Chapter 2. [START_REF] Lelièvre | Free energy computations : a mathematical perspective[END_REF]).

Lemma 3.1. Under conditions (3.3) and (3.4), the second step leaves the measure π invariant.

Proof of Lemma 3.1. The proof can be found in [START_REF] Lelièvre | Free energy computations : a mathematical perspective[END_REF]. It consists in seeing this step as a generalized Metropolis-Hastings step, with proposal Q(x ξ , y η ) = δ Φ ξ h (x) δ ξ (η), and symmetric operator S(x ξ ) = x -ξ . Then, it is enough to show that the Metropolis-Hastings ratio r defined by the following Radon-Nikodym derivative

r(x ξ , y η ) = Q(S(y η ), S(dx ξ ))π(dy η ) Q(x ξ , dy η )π(dx ξ ) , is equal to exp (log π(y) -log π(x)).
For example, for all x ξ ∈ E the centered point integrator Φ ξ h (x) defined by the solution y of the following equation

y = x -hξJ∇ log π x + y 2 , (3.5)
satisfies these properties, under some assumptions that ensure well-posedness.

Lemma 3.2. Under Assumption 1, there exists h 0 > 0 such that for all positive h < h 0 , there exists a unique solution to equation 3.5, and the integrator Φ ξ h is well-defined and satisfies equations (3.3) and (3.4).

Proof of Lemma 3.2. The proof follows on from Lemma 2.2.

The main benefit of this algorithm with respect to GMALA, is the better average acceptance ratio of the Hybrid step with the centered point integrator, which is of order O(h 3 ) instead of O(h 3/2 ) for GMALA, which may enable to reduce the rate of switching directions. Lemma 3.3. Suppose Assumption 2. Then, for all l ≥ 1, there exists h 0 > 0 such that for all positive h < h 0 , and for all x ∈ R d ,

(1 -β h,ξ (x)) 2l ≤ C ||∇U (x)|| 2l h 6l .
Proof of Lemma 3.3. For x ∈ R d , we set y = Φ ξ h (x). A Taylor expansion of U (y) yields,

|U (y) -U (x)| = O(h 3 D 3 U (x) ||∇U (x)|| 3 ).
The centered point integrator is an example of integrator for the Hamiltonian dynamics, that can be used as soon as the potential U is globally Lipschitz. Actually, similarly to GMALA in the non globally Lispchitz case, one can construct an approximate integrator Φ ξ h by using a truncation of ∇U to ensure its global Lipschitzness. Yet, this trick may lead to high rejection rates in the areas where ∇U is truncated. Again, other strategies can be used like importance sampling or even a change of variable in the integrand (1.1). The GHMALA algorithm is especially interesting when we are able to integrate efficiently the Hamiltonian dynamics. We show later on the numerical experimentations two examples of specific integrators that enable to improve significantly the performance of GHMALA with respect to GMALA with proposal kernel Q ξ 2 . We can recall for example, the case of separated Hamiltonian dynamics that can be integrated with explicit volume preserving schemes which are time-reversible and symmetric ( [START_REF] Qin | Volume-preserving schemes and numerical experiments[END_REF]).

The authors propose in [START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF] a similar splitted scheme where the hybrid step is replaced by fourth-order Runge-Kutta method. Even though this choice leads to a biased estimator, it enabled to get rid of the centered point integrator, that may be more costly than the Runge-Kutta method in terms of computation time.

Numerical experimentations

4.1. Anisotropic distribution. In this section, we test our nonreversible MALA algorithm by computing an observable f with respect to a two dimensional anisotropic distribution. More precisely we want to estimate E [f (X)] with f ((x 1 , x 2 )) = x 2 1 1 x1>15 , and X ∼ π(x) with π(x) ∝ e -V (x) , where

V ((x 1 , x 2 )) = x 2 1 1 + 50x 2 1 + x 2 2 .
Such a distribution is more stretched out in the x 1 direction rather than the x 2 . We expect our lifted algorithm to be more favorable than MALA when the anisotropy is strong. Indeed, MALA tends to perform a slower exploration of the state space in the x 1 direction with respect to x 2 the direction. The lifted algorithm is supposed to correct this problem since the Hamiltonian dynamics should lead to a fast exploration of the iso-potential lines, which are stretched in the direction x 1 . We choose as the skew-symmetric matrix J defined by,

J = α 0 1 -1 0 , (4.1)
with α a real parameter.

To compare MALA with GMALA and GHMALA, different optimal time steps parameters should be used. There is a tradeoff between achieving high average acceptance ratios (obtained with small time steps) and small correlation between the successive samples (obtained with large time-steps). More precisely MALA usually gives its best results with a large h that ensures a significant average rejection ratio. On the contrary, GMALA and GHMALA require much smaller time steps to avoid the regress that happens with the direction switching at each rejection. Thus, MALA should be used with a higher time step than the two others. Moreover, as GMALA with proposal kernel Q ξ 1 leads to a worse average acceptance ratio than Q ξ 2 , the former requires a smaller time-step than the latter. Figure 1 shows the average acceptance ratio with respect to the time step h for GMALA with proposal kernels Q ξ 1 and Q ξ 2 . We can verify that the average rejection ratio of GMALA with Q ξ 1 is indeed much worse than MALA and is of order h. In practice, this limitation leads to very correlated successive samples, and thus bad asymptotic variances for the estimators based on such Markov chain. In the following, we always consider GMALA with proposal kernel Q ξ 2 instead of Q ξ 1 . Figure 2 shows the acceptance ratio with respect to the time step h for those three algorithms. GH-MALA is composed by two steps: the first one consists of MALA and the second one consists of a Hybrid iteration. Both steps are composed by an acceptanceection step. The average acceptance rate of the first step does not depend on α, and is denoted by GHMALA MALA in the legend. Moreover, it is the same as a plain MALA. The average acceptance rate of the second step is denoted by GHMALA Hybrid. We also observe that the rejection ratio for GMALA with Q ξ 2 is close to MALA's and scales as h 3/2 This is due to the fact that the non reversibility contributes at order h 2 in the expression of the average rejection ratio. We can also verify the upper bound on the rejection probability of the hybrid step, given by Lemma 3.3. We compare now the asymptotic variance of the estimators build with MALA, GMALA and GHMALA. To do so, we compute 1000 independent estimators composed by the time average of 10 5 samples. We plot in Figure 3 the empirical relative variance of these estimators, We observe that GMALA performs much better than MALA by a factor 20. We should precise that a reduction in the variance does not necessarily mean a reduction in computing time since one iteration of GMALA is more costly than MALA since it requires a fixed point iteration.

The question of choosing α is quite natural. In the case of the time-continuous process, it is known that the decrease in asymptotic variance is monotonic with the intensity α, which suggests to use the algorithms with large α [START_REF] Rey | Irreversible langevin samplers and variance reduction: a large deviations approach[END_REF][START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF][START_REF] Hwang | Variance reduction for diffusions[END_REF]]. Yet, well-posedness of the proposal kernel Q ξ 2 requires the product αh to be small enough to ensure convergence of the fixed point. From a computational point of view, the cost of the method is proportional to the number of Picard iterations, which scales like -log ρ, where ρ denotes the contraction ratio (that scales as αh). Then, the choice of parameters α and h should take into account these two effects. 4.2. Warped Gaussian distribution. This example deals with a non quadratic potential. Both GMALA and GHMALA can be adapted to this case. The simple idea is to build a proposal kernel with a truncation of ∇U , to make it globally Lipschitz. This is the same idea as the one used for MALTA ( [START_REF] Gareth | Exponential convergence of langevin distributions and their discrete approximations[END_REF]). Moreover, it is also possible to choose a more efficient integrator than the centered point integrator defined by Equation (3.5).

To illustrate these two methods, we test now our algorithms in the case of a two-dimensional warped Gaussian distribution. This toy case has been introduced in [START_REF] Haario | An adaptive metropolis algorithm[END_REF] and used as a benchmark for variance reduction methods based on nonreversible Langevin samplers in [START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF]. More precisely, we aim to estimate E [f (X)] with the observable f and X distributed with π, defined by,

f (x) = x 2 1 + x 2 2 , π(x) ∝ e -V (x) , with V (x) = x 2 1 100 + x 2 + x 2 1 20 - 5 2
.

We define the skew symmetric matrix J by (4.1). It appeared in [START_REF] Duncan | Variance reduction using nonreversible langevin samplers[END_REF] that the nonreversible Langevin dynamics enables to reduce the asymptotic variance by several orders of magnitude.

We propose to implement GMALA using a truncated drift to make it globally Lipschitz to ensure the well-posedness of the method, and to implement GHMALA using a specific integrator of the Hamiltonian dynamics. We show that GHMALA performs better than GMALA and MALA in this casee. More precisely, the integrator is defined as the centered point integrator, after the symplectic change of variable ψ defined for all (x 1 , x 2 ) ∈ R 2 by ψ(x 1 , x 2 ) = (x 1 , x 2 +

x 2 1 20 -5).

Typically, this integrator enables to solve this dynamics for larger time steps than GMALA with proposal kernel Q ξ 2 , and thus reduces the asymptotic variance. Figure 4 displays the asymptotic variance for the estimators build with these algorithms. It is computed as the empirical variance of 2000 independent estimators constructed as the time average of 10 5 iterations of the algorithms. We can observe that GMALA and GHMALA performs similarly for small time steps h. Yet, for larger time-steps, it is not possible to define the proposal kernel for GMALA, whereas it is still the case for GHMALA. Eventually, we achieve a variance reduction of about a factor 500 with GHMALA and 60 with GMALA, compared with classical MALA. 4.3. Quartic Gaussian distribution. This toy case aims to present a particular case where GHMALA can be used without implicit integrator, and in a non globally Lipschitz case, which may enable to reduce the computational cost by avoiding the fixed point iteration. Again, we aim to estimate E [f (X)] where X is distributed with π, and where, f (x) = x 2 1 + x 2 2 , π(x) ∝ e -V (x) , with V (x) =

x 2 1 100 + x 4 2 .

We define the skew-symmetric matrix J by (4.1). In this case, the Hamiltonian dynamics defines then a separable system and volume-preserving explicit methods can be used (see [START_REF] Qin | Volume-preserving schemes and numerical experiments[END_REF]). More precisely, we define Φ ξ h for all x = (x 1 , x 2 ) ∈ R 2 by Φ ξ h (x) = (y 1 , y 2 ), where, We plot in figure 5 the asymptotic variance for the time average estimator build with GHMALA in the same way as for the warped Gaussian distribution. We can observe that the decrease in variance between MALA and GHMALA for small h is around 280. Nevertheless, for larger h, the explicit integration is not accurate enough, which leads to an increase in the asymptotic variance for larger time steps. Eventually, the smallest asymptotic variance of the time average estimator of GHMALA is around 50 times lower than the smallest asymptotic variance for MALA.

Conclusion

We presented a class of unbiased algorithm that enables to benefit from the variance reduction of the nonreversible Langevin equations (1.3) with respect to the reversible dynamics (1.2). More precisely, we presented two variations of these algorithms. The first one (GMALA) can be viewed as a lifting method, and more specifically as a generalized Metropolis Hastings methods on a lifted state space. The second one (GHMALA), similar to the first one, can be viewed as a Generalized Hybrid Monte-Carlo method.

Numerical experimentations show that variance reductions (compared with classical MALA) of several orders of magnitude can be achieved for potentials concentrated on a lower dimensional submanifold. We also expect these algorithms to perform better in the case of entropic barriers. The main difficulty is, in the case of GMALA, to use a proposal that enables to achieve a sufficiently high average acceptance ratio (to compete with MALA). For example this can be done by using a mid-point discretization. Even though this scheme is implicit, the computation of the Metropolis-Hastings acceptance probability does not require the computation of the Hessian of log π. In the case of GHMALA, numerical experimentations show that the choice of a suitable Hamiltonian integrator may lead to large improvements and computational cost reduction compared with the mid-point method.

Algorithm 3 . 1 (( 1 )

 311 Generalized Hybrid MALA). Let x 0 be an initial point. Set ξ = ±1. Let h > 0. Iterate on n ≥ 0, Integration of the reversible part (3.1):

Figure 1 .

 1 Figure 1. Comparison of average rejection ratio for GMALA with proposal kernels Q ξ 1 and Q ξ 2 .

Figure 2 .

 2 Figure 2. Comparison of average rejection ratio for MALA, GMALA (with proposal kernel Q ξ 2 ) and GHMALA.

Figure 3 .

 3 Figure 3. Variance comparison of MALA, GMALA (Q ξ 2 ) and GHMALA on the anisotropic distribution

Figure 4 .

 4 Figure 4. Variance comparison of MALA, GMALA and GHMALA on the warped Gaussian distribution

  non Lipschitz nonlinearity of ∇V is not an obstacle to the well-posedness of this integrator.

Figure 5 .

 5 Figure 5. Variance comparison of MALA and GHMALA on the quartic Gaussian distribution

Acknowledgement : This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH.