OPTIMAL SCALING OF THE MALA ALGORITHM WITH
IRREVERSIBLE PROPOSALS FOR GAUSSIAN TARGETS

arXiv:1702.01777v2 [stat.ME] 27 Nov 2017

MICHELA OTTOBRE, NATESH S. PILLAI, AND KONSTANTINOS SPILIOPOULOS
Abstract. It is well known in many settings that reversible Langevin diffusions in confining potentials converge to equilibrium exponentially fast. Adding irreversible perturbations
to the drift of a Langevin diffusion that maintain the same invariant measure accelerates
its convergence to stationarity. Many existing works thus advocate the use of such nonreversible dynamics for sampling. When implementing Markov Chain Monte Carlo algorithms (MCMC) using time discretisations of such Stochastic Differential Equations (SDEs),
one can append the discretization with the usual Metropolis-Hastings accept-reject step and
this is often done in practice because the accept–reject step eliminates bias. On the other
hand, such a step makes the resulting chain reversible. It is not known whether adding
the accept-reject step preserves the faster mixing properties of the non-reversible dynamics.
In this paper, we address this gap between theory and practice by analyzing the optimal
scaling of MCMC algorithms constructed from proposal moves that are time-step Euler
discretisations of an irreversible SDE, for high dimensional Gaussian target measures. We
call the resulting algorithm the ipMALA , in comparison to the classical MALA algorithm
(here ip is for irreversible proposal). In order to quantify how the cost of the algorithm
scales with the dimension N , we prove invariance principles for the appropriately rescaled
chain. In contrast to the usual MALA algorithm, we show that there could be two regimes
asymptotically: (i) a diffusive regime, as in the MALA algorithm and (ii) a “fluid” regime
where the limit is an ordinary differential equation. We provide concrete examples where the
limit is a diffusion, as in the standard MALA, but with provably higher limiting acceptance
probabilities. Numerical results are also given corroborating the theory.

Keywords: MALA Algorithm, Langevin Diffusions, Optimal Scaling, Non-reversible.
1. Introduction
In this paper, we analyze the scaling properties of high dimensional Markov Chain Monte
Carlo (MCMC) algorithms constructed using non-reversible Langevin diffusions. Consider a
target measure
Z
1 −U (x)
dx, Z =
e−U (x) dx.
(1.1)
π(dx) = e
Z
d
R
It is known that, under mild assumptions on the potential U (x), the Langevin stochastic
differential equation (SDE)
√
(1.2)
dXt = −∇U (Xt )dt + 2dWt
has π as its unique invariant measure and is π-ergodic. For appropriate test functions
f : Rd 7→ R, by the ergodic theorem, we have
Z
Z
1 t
f (Xs )ds −→ f (x)π(dx).
(1.3)
t 0
Date: November 28, 2017.
1

Thus the Langevin diffusion (1.1) is a fundamental
R tool for sampling from the target measure
π or compute expectation of various functionals f dπ, in view of (1.3). The Langevin SDE
is time reversible – its generator is a self-adjoint operator in L2 (π).
The drift of (1.2) can be modified without altering the invariant measure. Indeed, if Γ is
a vector field such that div(Γe−U ) = 0, then diffusions of the form
√
(1.4)
dXt = [−∇U (Xt ) + Γ(Xt )] dt + 2dWt ,
also have π as their invariant measure. The divergence free condition can be written as
divΓ = Γ∇U.

(1.5)

Observe that, for Γ 6= 0, the diffusions in (1.4) are non-reversible.
1.1. The problem. It is known that adding a non-reversible component in a Langevin
SDE could accelerate its convergence to stationarity. Indications of this phenomenon are
the main results of [10, 24, 25]. The main result in [10] states that, among the family of
non-reversible diffusions (1.4), the one with the smallest spectral gap corresponds to Γ = 0.
Similar results were established from an asymptotic variance and large deviations point
of view in [24, 25]. Broadly speaking, it is a well documented principle that non-reversible
dynamics have better ergodic properties than their reversible counterparts. This observation
has sparked a significant amount of research work in recent years and several papers have
advocated the use of non-reversible diffusions for sampling.
Numerical discretisations of (1.2) or (1.4) do not necessarily inherit the ergodic properties
of the continuous-time dynamics [28]. In particular, discretised processes may not converge
at all, or will have an invariant measure different from π. To circumvent these issues, in
MCMC algorithms, practitioners often perform an additional Metropolis-Hastings acceptreject step for proposals constructed from time discretisations of Langevin diffusions. For
instance, the standard MALA algorithm is obtained when the discretised SDE is (1.2) (see
e.g., [26] and [29]). To the best of our knowledge, not much is known about whether making
the chain reversible by adding the accept reject step preserves the faster mixing enjoyed by
the non-reversible dynamics. In this article, we seek to address this important gap between
theory and practice.
1.2. Previous Work. Many different approaches were pursued in recent works in order
to exploit and analyse the beneficial effects of irreversibility in the algorithmic practice
of MCMC methodology. In particular: i) Irreversible algorithms have been proposed and
analyzed in [2, 3, 5, 17, 13, 23] and references therein (on the matter see also [19, 20] ); ii)
algorithms that are obtained by discretising irreversible Markov processes in a way that the
resulting Markov chain is still irreversible are studied in [9], [21]; iii) numerical algorithms
that take advantage of the splitting of reversible-irreversible part of the equation are analyzed
in [16, 7]. In addition, comparisons of MALA and of Langevin samplers without the acceptreject step have been performed in [6, 7]. In many cases of interest it has been observed
that irreversible Langevin samplers have smaller mean square error when compared to the
MALA algorithm, see [7]. The latter fact is related to the consideration that the variance
reduction achieved by the irreversible Langevin sampler can be more significant than the the
error due to the bias of the irreversible Langevin sampler.
2

1.3. Our Contribution. In this paper we take a different standpoint and analyse the exact
performance of a Metropolis-Hastings algorithm, on certain Gaussian target densities, where
the proposal is based on discretising (1.4) as opposed to (1.2). We call this algorithm
the irreversible proposal MALA or ipMALA for short. To quantify the efficiency of the
ipMALA algorithm in high dimensions and compare it to the usual MALA algorithm,
we study its optimal scaling properties and its limiting optimal acceptance probability [26,
27, 22]. Optimal scaling aims to find the “optimal size” of the local proposal variance as
a function of the dimension. The optimality criteria varies for different algorithms, but a
natural choice for algorithms that have a diffusion limit is the expected square jumping
distance [27].
The basic mechanism for Metropolis-Hastings algorithms consists of employing a proposal
transition density q(x, y) in order to produce a reversible chain {xk }∞
k=0 which has the target
measure π as invariant distribution. At step k of the chain, a proposal move yk+1 is generated by using q(x, y), i.e., yk+1 ∼ q(xk , ·). Then such a move is accepted with probability
α(xk , yk+1 ), where


π(yk+1 )q(yk+1 , xk )
α(xk , yk+1 ) = min 1,
.
(1.6)
π(xk )q(xk , yk+1 )
The MALA algorithm is a Metropolis-Hastings algorithm with proposal move generated
by a time-step discretisation of the Langevin equation (1.2). The ipMALA algorithm
we wish to analyze obtains proposals by discretising (1.4). As explained before, any nontrivial Γ satisfying the divergence free condition (1.5) will preserve the invariant measure. A
convenient choice is to pick Γ such that
divΓ = 0

and Γ∇U = 0 .

(1.7)

A standard choice of Γ(x) is
Γ(x) = S∇U (x),

(1.8)

where S is any antisymmetric matrix. A more elaborate discussion on other possible choices
of Γ(x) can be found in [24]. The meaning of the conditions (1.7) is straightforward: the
flow generated by Γ must preserve Lebesgue measure since it is divergence-free; moreover,
the micro-canonical measure on the surfaces {U = z} is preserved as well.
We make one further important assumption. For the rest of the paper, we focus exclusively
on Gaussian target measures. We believe most of our analysis should carry over to the
important case in which the target measure has a Radon-Nikodym derivative with respect
to a Gaussian measure using the methods in [18, 22].
The main result of the paper is as follows. We consider Gaussian target measures π N ∼
N (0, C N ) on RN , where C N := diag{λ21 , ..., λ2N } (see Section 2 for more details). Such a
measure is clearly of the form (1.1) for a quadratic potential U . Therefore, in this case, the
general form of a Euler-discretisation of (1.4) is given by
2
σN
α N N N
N
xN + σN
C S xk + σN (C N )1/2 zk+1
,
2 k
N
where zk+1
∼ N (0, IN ), σN = `/N γ and α > 0. The notation S N is used to stress that here
N
N
S is an N ×N antisymmetric matrix. The quantity yk+1
is the proposal, and then a Markov
chain is formed using the usual accept-reject mechanism (more details about the algorithm
can be found in Section 3). Let xN
k be the resulting ipMALA Markov chain. Interestingly,
N
yk+1
= xN
k −

3

depending on whether α is bigger, equal or smaller than two, we will have different limiting
behaviours. Broadly speaking, we show the following:
• if γ < 1/6 then the acceptance probability degenerates to zero exponentially quickly,
and thus this case is not of practical interest.
• if α ≥ 2 then the optimal value of γ is γ = 1/6. We prove that the continuous
interpolant of two subsequent steps of the MALA algorithm, see (4.6) for proper
definition, has a diffusion limit. The cost of the algorithm is still N 1/3 like in the
standard MALA case. However, if S N and α are chosen appropriately (see specific
examples in Sections 6 and 7), the limiting acceptance probability is higher. That
is, the ipMALA will accept moves more frequently than MALA.
• if 1 ≤ α < 2 then take γ ≥ 1/6. In this regime we show that the cost of the algorithm
is of the order N αγ , and one can choose αγ < 1/3. In addition, we show that the
continuous interpolant, see (4.6), converges weakly to the solution of a deterministic
ODE. It is also interesting that the ODE that we get in the limit in this regime, can
be related to Hamilton’s ODE in HMC (see Theorem 4.2) and we plan to investigate
this in the future.
We would like to stress that in this paper we only study different scaling regimes when the
algorithm is started in stationarity. In the case of the MALA algorithm it is a known fact
that the optimal scaling out of stationarity differs from the optimal scaling in stationarity
(the former being O(N 1/2 ) and the latter being O(N 1/3 ), see [4, 11, 12, 15]). It would be
relevant to address the same issue for the ipMALA algorithm (especially in view of the fact
that non-reversibility is known to speed up convergence, although we do realise that the
measure of efficiency we use here may well be unrelated to convergence rate) but we do not
do it in this paper, as this would involve substantial further analysis.
The goal of this paper is to explore in a rigorous mathematical way what happens when
a Metropolis-Hastings accept-reject step is applied to a proposal coming from discretisation
of an irreversible Langevin diffusion. We find that there are different possible regimes and
we characterize explicitly what is the limit of the continuous interplant of the chain as
the number of steps goes to ∞, see Theorem 4.2. We find that the irreversible perturbation
matrix has non-trivial effects on the limiting dynamics (see Theorem 4.2) and on the limiting
average optimal acceptance probability, see (5.29) and Section 6. In terms of applications,
the conclusion is that even though the introduction of the accept-reject step potentially
offsets some of the advantages of irreversible perturbations [10, 24, 25], some advantages
may appear in certain cases. In particular, by appropriately choosing S and α (e.g. S = S1
and α > 4 as in Section 6), we can obtain a diffusion limit but with limiting acceptance
probabilities higher than that of MALA. In addition, even in the fluid regime, where the
limit of the algorithm is an ODE, the acceptance probability can be the same as that of
MALA, but the algorithm takes, in stationarity, N αγ steps to explore the state space with
αγ < 1/3 as opposed to N 1/3 steps for standard MALA.
Lastly, it is possible that for multi-modal targets the use of irreversible proposals may be
beneficial as this is the case in the absence of the accept-reject step, see [25]; an analysis
of this fact is though beyond the scope of this paper. The present paper certainly invites
further research on this topic in different directions, including multi-modal targets, out of
stationarity analysis and computational considerations of implementation issues.
4

The paper is organized as follows. In Section 2 we introduce the notation used throughout
the paper. In Section 3 we describe and motivate the algorithm that we will examine. Section
4 contains the rigorous statement of our main result and its main implications. In Section 5
we give a detailed heuristic argument to explain how the main result is obtained and why it
should hold. In Section 6 we present some examples of potential choices for the antisymmetric
matrix and the computation of the corresponding limiting acceptance probabilities. Section
7 contains extensive simulation studies that demonstrate the theoretical results. Rigorous
proofs are relegated to the Appendix.
2. Preliminaries and Notation
Let (H, h·, ·i, k · k) denote an infinite dimensional separable Hilbert space with the canonical norm derived from the inner-product. Let C be a positive, trace class operator on H and
{φj , λ2j }j≥1 be the eigenfunctions and eigenvalues of C respectively, so that
X
Cϕj = λ2j ϕj ,
λ2j < ∞.
(2.1)
j

The eigenvalues λj are non-decreasing. We assume a normalization under which {φj }j≥1
forms a complete orthonormal basis in H. Let π be a Gaussian probability measure on H
with covariance operator C, that is,
π ∼ N (0, C).
If X N is the finite dimensional space
H ⊃ X N := span{ϕj }N
j=1
spanned by the first N eigenvectors of the covariance operator (notice that the space X N is
isomorphic to RN ), then for any fixed N ∈ N and x ∈ H we denote by P N (x) the projection
of x on X N . Moreover, the finite dimensional projection π N of the measure π on X N is given
by the Gaussian measure
π N ∼ N (0, C N )
where C N := diag{λ21 , . . ., λ2N } (or, more precisely, C N := P N ◦ C ◦ P N ). Thus we have,
i,N 2

N
|x |
Y
1
1 − 2λ2
N
N
1,N
N,N
i
.
π := π(x ) = π(x , . . ., x ) = √
e
( 2π)N i=1 λi

(2.2)

Let S̃ : H → H be a bounded linear operator. Thus there exists a constant κ > 0 such
that
kS̃xk ≤ κkxk,
x ∈ H.
(2.3)
For any N ∈ N we can consider the projected operator S̃ N , defined as follows:
S̃ N x := (P N ◦ S̃ ◦ P N )x.
The operator S̃ N is also bounded on H. Since X N is isomorphic to RN , S̃ N can be represented
by an N × N matrix. Throughout the paper we require the following: S̃ is such that for any
N ∈ N, the matrix S̃ N can be expressed as the product of a symmetric matrix, namely C N ,
and an antisymmetric matrix, S N :
S̃ N = C N S N ,

for every N ∈ N.
5

(2.4)

Throughout the paper we will use the following notation:
• x and y are elements of the Hilbert space H;
• the letter N is reserved to denote the dimensionality of the space X N where the
target measure π N is supported;
• xN is an element of X N ∼
= RN (similarly for y N and the noise ξ N ); the j-th compoj,N
nent of the N -dimensional vector xN (in the basis {ϕj }N
.
j=1 ) is denoted by x
We analyse Markov chains evolving in RN . Because the dimensionality of the space in
which the chain evolves will be a key fact, we want to keep track of both the dimension N
and the step k of the chain. Therefore,
N
N
• xN
k will denote the k-th step of the chain {xk }k∈N evolving in R ;
i,N
• compatibly with the notation set above, xk is the i-th component of the vector
i,N
N
= hxN
xN
k , ϕi i;
k ∈ R ; i.e., xk
N
N
N
• two (double) sequences of real numbers {AN
k } and {Bk } satisfy Ak . Bk if there
exists a constant K > 0 (independent of N and k) such that
N
AN
k ≤ KBk ,
N
for all N and k such that {AN
k } and {Bk } are defined. The notation . is used for
functions and random functions as well, with the constant K being independent of
the argument of the function and on the randomness.

As we have already mentioned, k · k is the norm on H, namely
2

kxk :=

∞
X

hx, ϕi i2 .

i=1

With abuse of notation, we will also write
N 2

kx k =

N
X

xi,N

2

.

i=1

We will also use the weighted norm k · kC , defined as follows:
kxk2C

:= hx, C

−1

xi =

∞
X
|hx, ϕi i|2

λ2i

i=1

,

for all x ∈ H such that the above series is convergent; analogously,
kxN k2C N

N
X
xi,N
:=
λ2i
i=1

2

.

If AN is an N × N matrix, AN = {AN
i,j }, we define
VAN

:=

N
X
i=1

λ2i

N
X

2

2
AN
i,j λj

j=1

= EπN k(C N )1/2 AN xN k2 = EπN h(C N )1/2 z N , AN xN i
6

2

(2.5)

In the above z N ∼ N (0, IN ), xN ∼ π N and EπN denotes expectation with respect to all the
sources of noise contained in the integrand. We will often also write Ek and Ex , to mean
Ek = E[·|xk ],

Ex = E[·|xk = x].

Finally, many of the objects of interest in this paper depend on the matrix S N and on the
parameters α and γ. When we want to stress such a dependence, we will add a subscript p ,
see for example the notation for the drift dp in Theorem 4.2.
3. The ipMALA algorithm
As we have already mentioned, the classical MALA algorithm is a Metropolis-Hastings algorithm which employs a proposal that results from a one-step discretisation of the Langevin
dynamics (1.2). This is motivated by the fact that the dynamics (1.2) is ergodic with unique
invariant measure given by π (1.1) and can therefore be used to sample from π. The Langevin
dynamics that samples from our target of interest, the Gaussian measure π N in (2.2), reads
as follows:
√
Xt ∈ X N .
dXt = −(C N )−1 Xt dt + 2dWt ,
If ∆ is any positive definite, N -dimensional symmetric matrix, then the equation
√
dXt = −∆(C N )−1 Xt dt + 2∆dWt
is still ergodic with invariant invariant measure π N in (2.2). Now notice that if C is a trace
class operator, then C −1 is unbounded, so in the limit as N → ∞ any of the above two
dynamics (and their discretisations) would lead to numerical instabilities. To avoid the
appearance of unbounded operators, we can choose ∆ = C N /2 and therefore consider the
SDE
1
dXt = − Xt dt + (C N )1/2 dWt .
2
Discretizing the above and using such a discretisation as a Metropolis-Hastings proposal
would result in a well-posed MALA algorithm to sample from (2.2). However, as explained
in the Introduction, here we want to analyze the MALA algorithm with irreversible proposal.
As in (1.8), we next consider the non-reversible SDE:


1
N
N −1
dXt = − Xt + S (C ) Xt dt + (C N )1/2 dWt ,
2
where S N is any N × N antisymmetric matrix. Again to avoid the appearance of unbounded
operators, we modify the irreversible part of drift term and, finally, obtain the dynamics


1
N N
(3.1)
dXt = − Xt + C S Xt dt + (C N )1/2 dWt .
2
Notice that, for any xN ∈ X N ,

∇ · (C N S N xN ) π N = Trace(C N S N )π N − hC N S N xN , (C N )−1 xN iπ N = 0,

(3.2)

having used the antisymmetry of S N and the symmetry of C N . Therefore, π N is invariant
for the dynamics (3.1). This justifies using the Metropolis-Hastings proposal (3.3) below,
which is a one-step Euler discretisation of the SDE (3.1).
7

We now describe the ipMALA algorithm. If the chain is in xN
k at step k, the ipMALA
algorithm has the proposal move:
N
yk+1
= xN
k −

2
σN
α N N N
N
xN + σN
C S xk + σN (C N )1/2 zk+1
2 k

(3.3)

where
`
,
`, γ > 0
(3.4)
Nγ
and α > 0. The proposal (3.3) is a generalised discretisation of (3.1). We may indeed choose
α so that, asymptotically, relative to the reversible drift and diffusion, the non-reversible
drift dominates (1 ≤ α < 2), vanishes (α > 2) or is balanced (α = 2). This will result in
different scaling limits. The choice of the parameters α, γ and ` will be further discussed in
1,N
N,N
N
= (zk+1
, . . ., zk+1
) is
Section 4 below. S N can be any antisymmetric N × N matrix and zk+1
a vector of i.i.d standard Gaussians. With this proposal, the acceptance probability is
σN =

β N (xN , y N ) = 1 ∧ eQ

N (xN ,y N )

(3.5)

where
2
 1 2α−2

σN
ky N k2C N − kxN k2C N + σN
kC N S N xN k2C N − kC N S N y N k2C N
8
2
α−2
N N
N
+ 2σN hS y , x i.

QN (xN , y N ) = −

The above expression for QN is obtained with straightforward calculations, after observing
that the proposal kernel q(xN , y N ) implied by (3.3) is such that


2
σN
1
N
N
N
α N N N 2
N
N
x − σN C S x kC N .
q(x , y ) ∝ exp − 2 ky − x +
2σN
2
If β̃ N ∼ Bernoulli(β N (xN , y N )), then the Metropolis-Hastings chain {xN
k }{k∈N} resulting from
using the proposal (3.3) can be written as follows
N N
N
N
N
N
N
N
xN
k+1 = β̃ yk+1 + (1 − β̃ )xk = xk + β̃ (yk+1 − xk ).

(3.6)

4. Main result and Implications
To understand the behaviour of the chain (3.6), we start by decomposing it into its drift
and martingale part; that is, let ζ > 0 and write
N
xN
k+1 = xk +

1
1 N N
N
dp (xk ) + ζγ/2 Mk,p
ζγ
N
N

(4.1)

where the approximate drift dN
p is defined as
N
ζγ
N
N
dN
p (xk ) := N Ek (xk+1 − xk )
N
and the approximate diffusion Mk,p
is given by


N
N
N
N
Mk,p
:= N ζγ/2 xN
k+1 − xk − Ek (xk+1 − xk ) .

8

(4.2)

Using (3.3)-(3.4) and (3.6), we rewrite the approximate drift (4.2) as follows:



`2 N
`α N N N
N
N
ζγ
QN
dp (xk ) = N Ek (1 ∧ e ) − 2γ xk + αγ C S xk
2N
N


`
N
N
.
+ N ζγ Ek (1 ∧ eQ ) γ (C N )1/2 zk+1
N

(4.3)
(4.4)

Looking at (4.3) it is clear that we need to examine two different cases, namely
i) Diffusive regime : 2γ ≤ αγ, i.e., α ≥ 2
ii) Fluid regime : 2γ > αγ, i.e., α < 2.
The names of the above regimes will be clear after the statement of our main results, Theorem
4.1 and Theorem 4.2, see also Remark 4.4. We will show that, in order for the algorithm to
have a well defined non-trivial limit, in the case i) we should choose ζ = 2, whereas in case
ii) we should choose ζ = α. For this reason, it is intended from now on that
ζ = α if α < 2

and

ζ = 2 if α ≥ 2.

(4.5)

With this observation in mind, we introduce the continuous interpolant of the chain {xN
k }:
ζγ
N
x(N ) (t) = (N ζγ t − k)xN
k+1 + (k + 1 − N t)xk ,

tk ≤ t < tk+1 ,

(4.6)

where tk = k/N ζγ .
We now come to state the main results of this paper, Theorem 4.1 and Theorem 4.2.
Theorem 4.2 is the most general result. The assumptions under which Theorem 4.2 holds
are a bit involved; we detail and motivate them in Section 5.1. In this section we first
state Theorem 4.1, which is just Theorem 4.2, adapted to the case in which the sequence of
matrices S N is chosen from a specific class. Namely, we restrict our attention to the case
N
= Ji and
in which each of the matrices S N is in Jordan block form: we assume that Si,j
N
N
Sj,i = −Ji for j = i + 1 and i = 1, 3, 5, 7, · · · , N − 1 and Si,j = 0 otherwise; here {Ji }i is
an aribitrary sequence of real numbers. To be more clear, when we refer to matrices S N in
Jordan block form, we are referring to matrices of the form


0 J1 0
0 0 0
 −J1 0
0
0 0 0 


 0
0
0 J2 0 0  .
(4.7)


0 −J2 0 0 0 
 0
.
0
0
0
0 ..
In this case, we can carry out explicit calculations, construct examples and demonstrate our
results.
Theorem 4.1. Assume that the family of anti-symmetric matrices S N is in Jordan block
form (i.e. of the form 4.7 as discussed above) and x0 ∼ π. Assume the constant
EπN kS̃ N xN k2C N
c1 = lim
is finite.
(4.8)
N →∞
N 2(α−1)γ
Then, as N → ∞, the continuous interpolant x(N ) (t) of the chain {xN
k } (defined in (4.6)-(4.5)
N
and (3.6), respectively) started from the state xN
=
P
(x
),
converges
weakly in C([0, T ]; H)
0
0
to the solution of the following equation
dx(t) = dp (x(t))dt + Dp dWtC ,
9

x(t) ∈ H,

(4.9)

where WtC is a H-valued Brownian motion with covariance C. The drift coefficient dp (x) :
H → H and the diffusion constant Dp ≥ 0 are defined as follows:

2
 − `2 hJp x if γ = 1/6 and α ≥ 2
(4.10)
dp (x) :=

0
if γ ≥ 1/6 and 1 ≤ α < 2,
and
( q
` hJp if γ = 1/6 and α ≥ 2
Dp :=
(4.11)
0
if γ ≥ 1/6 and 1 ≤ α < 2 .
The real constant hJp is defined as
hJp := 2Φ

−`6 /32 − a
p
(`6 /16) + 2a

!
,

(4.12)

where a = 2`2(α−1) c1 . Lastly, as measured by the number of steps, in stationarity, to explore
the state space the cost of the algorithm in the diffusion case, α ≥ 2, is N 1/3 whereas in the
fluid case, 1 ≤ α < 2, it is N αγ .
We now move on to stating Theorem 4.2. The precise statement of the assumptions of the
theorem is deferred to Section 5. After stating Theorem 4.2, we first compare it with Theorem
4.1 (Remark 4.3) and then make several observations regarding practical implications of the
statement (Remark 4.4). We stress that the operator S̃ appearing in the statement of the
theorem is as in Section 2 (see (2.3)-(2.4)). In particular, the projected operator S̃ N is, for
every N , the product of the matrix C N and of a generic antisymmetric matrix S N , appearing
in the proposal (3.3).
Theorem 4.2. Let Assumption 5.1, Assumption 5.2, Condition 5.3 and Assumption 5.4
hold and let x0 ∼ π. Then, as N → ∞, the continuous interpolant x(N ) (t) of the chain {xN
k }
N
N
(defined in (4.6)-(4.5) and (3.6), respectively) started from the state x0 = P (x0 ), converges
weakly in C([0, T ]; H) to the solution of the following equation
dx(t) = dp (x(t))dt + Dp dWtC ,

x(t) ∈ H,

(4.13)

where WtC is a H-valued Brownian motion with covariance C. In addition, the drift coefficient
dp (x) : H → H and the diffusion constant Dp ≥ 0 are defined as follows:
 `2
− 2 hp x
if γ = 1/6 and α > 2





2
dp (x) :=
(4.14)
− `2 hp x + τp `2 S̃x if γ = 1/6 and α = 2





τp `α S̃x
if γ ≥ 1/6 and 1 ≤ α < 2,
and
 p
` hp if γ = 1/6 and α ≥ 2
Dp :=
0
if γ ≥ 1/6 and 1 ≤ α < 2

(4.15)

The real constants hp and τp are defined in (5.28) and in (5.32) respectively.
Proof. See Appendix A and B.


10

Remark 4.3. If the sequence of matrices S N is in Jordan block form, then all the assumptions of Theorem 4.2 are satisfied, provided (4.8) holds. 1 It should be clearly said that,
while examples and calculations in the case of Theorem 4.1 are very explicit, we have not
been able to find examples of matrices S N such that the constant τp in Theorem 4.2 is nonvanishing. However it is important to stress that even when τp = 0 it is still possible that
the average acceptance probability (given, in the limit, by hp ) is not tending to one and it
is of order 1 instead, see the examples in Sections 6 and 7. Notice that τp = 0 when b = 2a,
with b and a are parameters depending on the choice of S N , α and γ, defined in (5.26) and
(5.27). Cfr Table 2 and Table 1 to have a summary of how hp and τp depend on a and b.
We now comment on the practical implications of Theorem 4.2.
Remark 4.4. The limiting drifts (and corresponding diffusion coefficients) appearing in
(4.14) correspond to the two different regimes i) and ii), which we identified before the
statement of the theorem. The choice γ = 1/6 in case (i) (i.e., in the diffusive regime)
and γ ≥ 1/6 in the fluid regime will be fully motivated in Section 5.1. The various results
obtained in the different regimes are summarised in Tables 1 and 2 below.
• In the regime (i), the effective time-step implied by the interpolation (4.5)–(4.6), is
N −2γ = N −1/3 . Therefore, if γ = 1/6 and α ≥ 2, Theorem 4.2 implies that the
optimal scaling for the proposal variance when the chain is in its stationary regime
is of the order N −1/3 . That is, the cost of the algorithm (in terms of number of steps
needed, in stationarity, to explore the state space) is O(N 1/3 ). This is the same
scaling as obtained in [26, 22] for the MALA algorithm. Therefore, in the regime
(i), the ipMALA algorithm has the same scaling properties of the usual MALA.
More discussions on this can be found in Section 5.2.
• In the case γ = 1/6 and α > 2, one can construct specific matrices S, which result
in considerably higher limiting acceptance probabilities than in the classical MALA
case, see Sections 6 and 7
• When α < 2 and γ ≥ 1/6, the (rescaled) chain converges to a fluid limit; i.e., the
limiting behaviour of the chain is described by an ODE. Such an ODE still admits
our target as invariant measure. With the same reasoning as above, the main result
implies that, in this regime, the cost of the algorithm is of the order N αγ and one
can choose αγ < 1/3. In the fluid regime, one may expect that the chain gets to
get stuck in the limit as N → ∞. However, we also comment here that in the
simulations that we performed we did not see the chain getting stuck in this regime,
which may imply that one may need to go to very high dimensions to observe such
issues.
5. Heuristic derivation of the diffusion limit
In this section we give heuristic arguments to explain how one can formally obtain the
diffusion limit of Theorem 4.2 for the chain {xN
k }k . We stress that the arguments of this
section are only formal; therefore, we often use the notation “ ' ”, to mean “approximately
equal”. We write A ' B when A = B+ “terms that are negligible” as N tends to infinity;
1Indeed,

in this case Assumption 5.1 and Assumption 5.2 are satisfied with c1 as in (4.8) and c2 = c3 = 0.
This makes Assumption 5.4 easy to verify. Moreover, Condition 5.3 is trivially satisfied for matrices in
Jordan block form. Detailed comments on this can be found in Section 5, see comments after (5.29).
11

Setting of Theorem 4.1: cj = 0 for j = 2, 3 and b = 2a
γ = 1/6 ⇒ acc. prob
q is O(1)
α≥2

2

dxt = − `2 hJp xt + `

hJp dWtC

N 1/3

dxt = 0
N αγ
γ > 1/6 ⇒ acc. prob tends to one
α≥2
too costly N αγ (> N 1/3 )
1≤α<2
dxt = 0
N αγ
Table 1. The possible different regimes when the antisymmetric matrix S
takes the Jordan block form. The constants cj are those appearing in Assumption 5.2; such constants also determine the value of the parameters a and
b (defined after (5.25))
1<α<2

Setting of Theorem 4.2: case cj 6= 0 for j = 1, 2, 3 and b 6= 2a
γ = 1/6 ⇒ acc. prob is O(1)
p
2
α>2
dxt = − `2 hp xt + ` hp dWtC
N 1/3
p
`2
α=2
dxt = − 2 hp xt + τp `2 S̃xt + ` hp dWtC
N 1/3
α
1≤α<2
dxt = τp ` S̃xt
N αγ
γ > 1/6 ⇒ acc. prob is O(1)
α≥2
too costly N αγ (> N 1/3 )
1≤α<2
dxt = τp `α S̃xt
N αγ
Table 2. The above table summarizes the general scenario described by Theorem 4.2 in the case in which cj 6= 0 for j = 1, 2, 3 and b 6= 2a. The constants
cj are those appearing in Assumption 5.2; such constants also determine the
value of the parameters a and b (defined after (5.25)). Table 1 is a subcase of
this table. Indeed, Theorem 4.2 also covers the following: if b = 2a and cj 6= 0
for at least one j then all of the above holds with τp = 0. If b = 2a and all the
cj ’s vanish, then all of the above holds with τp = 0.

we then rigorously justify these approximations, and the resulting limit theorems, in the
Appendix.
5.1. Study of the acceptance probability and statement of main assumptions. In
order to understand the behaviour of the chain, it is crucial to gain intuition about the
acceptance probability β N . While attempting to improve such an intuition in this section,
we also present the motivation behind the assumptions needed for Theorem 4.2.
Using (3.3), a more useful and detailed expression for QN (which was introduced in (3.5))
is
N
N
QN (xN , y N ) = Q̄N (xN , y N ) + QN
α (x , y )

12

(5.1)

where Q̄N contains all the terms that come from the reversible part of the proposal while
QN
α contains all the terms that come from the irreversible part of the proposal; namely,

 3
5

σN
σ4
σN
N
+
hxN , (C N )1/2 z N iC N + N kxN k2C N − kz N k2
Q̄ := −
4
8
8
6
σ
− N kxN k2C N
32
and
QN
α


 
1 3+α
1 2α 
α−1
N
N 1/2 N 2
N 2
:= σN kS̃x kC N − kS̃ (C ) z kC N − 2σN + σN
h(C N )1/2 z N , S N xN i
2
4
(5.2)




1 2α+1
1 2α+2
2α−1
2α−2
− σN
− 2σN
+ σN
kS̃ N xN k2C N − σN
hS̃ N (C N )1/2 z N , S̃ N xN iC N
4
2
(5.3)
1 4α−2
3α−1
− σN
hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC N − σN
k(S̃ N )2 xN k2C N
(5.4)
2


2
σN
3α−2
hS̃ N xN , (S̃ N )2 xN iC N .
(5.5)
− σN
1−
2

α−1
3+α
For N large, from (3.4) we have σN
>> σN
; so we expect that, asymptotically, the last
term in (5.2) will be negligible. The same reasoning can be applied to the terms in (5.3).
Therefore, irrespective of the choice of γ and α, we have the approximation


1 2α  N N 2
α−1
N
N 1/2 N 2
QN
'
σ
k
S̃
x
k
−
k
S̃
(C
)
z
k
− 2σN
h(C N )1/2 z N , S N xN i
N
N
α
C
C
2 N
2α−2
2α−1
− 2σN
kS̃ N xN k2C N − σN
hS̃ N (C N )1/2 z N , S̃ N xN iC N
1 4α−2
3α−1
k(S̃ N )2 xN k2C
− σN
hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC − σN
2


2
σN
3α−2
− σN
1−
hS̃ N xN , (S̃ N )2 xN iC N .
2

(5.6)

We further observe that, in stationarity (and again irrespective of the choice of γ, α and S N )
the term

1 2α  N N 2
σN kS̃ x kC N − kS̃ N (C N )1/2 z N k2C N
2
is smaller than the term
2α−2
−2σN
kS̃ N xN k2C N

in the sense that it is always of lower order in N . Moreover, due to the skew-symmetry of
S N , the term (5.6) is identically zero:
hS̃ N xN , (S̃ N )2 xN iC N = hC N S N xN , C N S N C N S N xN iC N = hC N S N xN , S N C N S N xN i = 0 .
13

Therefore we can make the further approximation
α−1
2α−2
N 1/2 N
QN
z , S N xN i − 2σN
kS̃ N xN k2C N
α ' − 2σN h(C )
1 4α−2
2α−1
k(S̃ N )2 xN k2C N
− σN
hS̃ N (C N )1/2 z N , S̃ N xN iC N − σN
2
3α−1
− σN
hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC N =: RαN .

(5.7)
(5.8)
(5.9)

As a result of the above reasoning we then have the heuristic approximation
N
QN
α ' Rα .

(5.10)

We now heuristically argue that the only sensible choice for γ is γ = 1/6 if α ≥ 2 and
γ ≥ 1/6 if α < 2. In order to do so we look at the decomposition (5.1) of QN and recall the
definition (3.5) of the acceptance probability; we then write
N

N

β N = 1 ∧ (eQ̄ eQα ).
N
Let us start by looking at Q̄N : if we start the chain in stationarity, i.e. xN
then
0 ∼ π
N
N
N
N
N
xk ∼ π for every k ≥ 0. In particular, if x ∼ π then x can be represented as
P
xN = N
i=1 λi ρi ϕi , where ρi are i.i.d. N (0, 1). We can therefore use the law of large numbers
and observe that
N
X
kxN k2C N =
|ρi |2 ' N.
i=1

Similarly, by the Central Limit Theorem (CLT) the term hxN , (C N )1/2 z N iC N is O(N 1/2 ) and
converges to a standard Gaussian as N → ∞. Again by the CLT we can see that also the
term kxN k2C N − kz N k2 is O(N 1/2 ) in stationarity. Therefore, recalling (3.4), one has
`3
`6
N
N 1/2 N
N
hx
,
(C
)
z
i
−
kxN k2C N .
C
4N 3γ
32N 6γ
With these observations in place we can then argue the following:
• If γ = 1/6 then Q̄N is O(1) (in particular, for large N it converges to a Gaussian
N
with finite mean and variance, see (5.19)); therefore eQ̄ is O(1) as well
N
• If γ < 1/6 then Q̄N → −∞ (more precisely EQ̄N → −∞), hence eQ̄ → 0
N
• If γ > 1/6 then Q̄N → 0, hence eQ̄ → 1
Let us now informally discuss the acceptance probability β N in each of the above three cases,
taking into account the behaviour of QN
α as well. To this end it is worth noting the following:
Q̄N ' −

EπN kS̃ N xN k2C N = EπN k(C N )1/2 S N xN k2 = VSN = EπN hC 1/2 z N , S N xN i
EπN k(S̃ N )2 xN k2C N = VS S̃ = EπN hS̃ N (C N )1/2 z N , S̃ N xN iC N

2

(5.11)

2

,

(5.12)

where the first equality in (5.11) holds because
k(C N )1/2 S N xN k2 = kS̃ N xN k2C N

for all x,

and the others follow from (2.5). By (5.11)- (5.12) and (5.7) - (5.9), as N → ∞, the (average,
given x, of the) quantity QN
α can only tend to zero, go to −∞ or be O(1), but it cannot
diverge to +∞. With this in mind,
a): If γ = 1/6 (that is, Q̄N is O(1)) then one has the following two sub-cases
14

N
a1): γ = 1/6 and either QN
α → 0 or Qα is O(1). In this case the limiting
acceptance probability does not degenerate to zero (and does not tend to one)
N
a2): γ = 1/6 and QN
→ 0.
α → −∞. In this case β
N
Q
N
b): If γ < 1/6 then e can only tend to zero (if Qα → 0, to −∞ or is O(1)). This is
not a good regime as in this case β N would overall tend to zero.
c): If γ > 1/6 we have three sub-cases
N
c1): γ > 1/6 and QN
is O(1).
α is O(1). In this case β
N
N
c2): γ > 1/6 and Qα → 0, in which case β → 1.
N
degenerates to zero.
c3): γ > 1/6 and QN
α → −∞ , in which case β
N
We clearly want to rule out all the cases in which β tends to zero, (because this implies
that the chain is getting stuck). Therefore, if we want to “optimize” the behaviour of the
chain, the only good choices are a1), c1) and, potentially c2). However, considering that the
cost of running the chain is either 2γ (when α ≥ 2) or αγ (when α < 2), the case γ > 1/6
is always going to be sub-optimal if α ≥ 2; which is why in Theorem 4.2 we only consider
γ > 1/6 if α < 2.
We therefore want to make assumptions on S N , γ and α to guarantee that our analysis
falls only into one of the good regimes, i.e., to guarantee that QN
α is of order one (or tends
to zero). Assumption 5.1 and Assumption 5.2 below are enforced in order to guarantee that
this is indeed the case.
To explain the nature of Assumption 5.1 below, we recall the approximation (5.10) and
the expression for RαN , namely:

`α−1
`2(α−1)
N 1/2 N
N N
h(C
)
z
,
S
x
i
−
2
kS̃ N xN k2C N
(α−1)γ
2(α−1)γ
N
N
1 `2(2α−1)
`2α−1
− (2α−1)γ hS̃ N (C N )1/2 z N , S̃ N xN iC N −
k(S̃ N )2 xN k2C N
2(2α−1)γ
N
2N
`3α−1
− (3α−1)γ hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC N
(5.13)
N
Roughly speaking, Assumption 5.1 requires that a CLT should hold for all the scalar
products in the above expressions, Assumption 5.2 requires a Law of Large Numbers (LLN)
to hold for the remaining terms. We will first state such assumptions and then comment on
them in Remark 5.5.
We recall that S̃ N = C N S N so, since C N is fixed by the problem and invertible, the
following assumptions can be equivalently rephrased in terms of either S̃ N or S N .
RαN = − 2

Assumption 5.1. The entries of the matrix S N do not depend on N . There exist positive
constants 0 ≤ d1 , d2 , d3 < ∞ such that the sequence of matrices {S N }N (and the related
{S̃ N }N ) satisfies the following conditions, for some α ≥ 1:
1
D
i)
h(C N )1/2 z N , S N xN i −→ N (0, d1 )
(α−1)γ
N
1
D
ii)
hS̃ N (C N )1/2 z N , S̃ N xN iC N −→ N (0, d2 )
N (2α−1)γ
1
D
iii)
hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC N −→ N (0, d3 ),
(3α−1)γ
N
D

where −→ denotes convergence in distribution.
15

Assumption 5.2. The parameter α ≥ 1 and the sequence of matrices {S N }N are such that
2

i) lim EπN

kS̃ N xN k2C N
− c1
N 2(α−1)γ

ii) lim EπN

k(S̃ N )2 xN k2C N
− c2
N 2(2α−1)γ

N →∞

N →∞

iii) lim

N →∞

EπN

=0

(5.14)

2

=0

k(S̃ N )T (S̃ N )2 xkC N
− c3
N 2(3α−1)γ

2

= 0,

for some constants 0 ≤ c1 , c2 , c3 < ∞. We do not consider sequences {S N }, such that (5.14)
is satisfied with c1 = 0 and α = 1. That is, we consider sequences {S N }N such that (5.14)
is satisfied with c1 = 0, only in the cases in which this happens for α > 1, see Remark 5.5.
In view of (5.11)-(5.12), Assumption 5.2 implies the following three facts
EπN kS̃ N xN k2C N
= c1
lim
N →∞
N 2(α−1)γ
EπN k(S̃ N )2 xN k2C N
= c2
lim
N →∞
N 2(2α−1)γ
k(S̃ N )T (S̃ N )2 xN kC N
lim EπN
= c3 .
N →∞
N 2(3α−1)γ

(5.15)
(5.16)
(5.17)

Let us state a sufficient condition under which Assumption 5.1 holds.
Condition 5.3. Let z N ∼ N (0, IdN ) and xN ∼ π N . The sequence of matrices {S N }N is
such that for every h = 1, . . ., N the σ-algebra generated by the random variable z h,N (S N xN )h
is independent of the σ-algebra generated by the sum
h−1
X

z j,N (S N xN )j .

j=1

In the above we stress that (S N xN ) is a vector in RN and (S N xN )h is the h-th component of
such a vector. In order for Assumption 5.1 to be satisfied, the same should hold for (S̃ N )T S̃ N
and for (S̃ N )T (S̃ N )2 as well.
Notice that, because C N is diagonal, if (S N )2 and (S N )3 satisfy Condition 5.3, then also
(S̃ N )T S̃ N and (S̃ N )T (S̃ N )2 do. We give an example of a matrix satisfying all such conditions
in Remark 5.5. Finally, for technical reasons (see proof of statement (ii) of Lemma B.2, we
also assume the following.
Assumption 5.4.
(i) Under π N , the variables appearing in (5.15)- (5.17) have exponential moments, i.e., there exists a constant m, independent of N , such that
EπN e

kS̃ N xN k2 N
C
N 2(α−1)γ

, EπN e

k(S̃ N )2 xN k2 N
C
N 2(2α−1)γ

16

, EπN e

k(S̃ N )T (S̃ N )2 xN k N
C
N 2(3α−1)γ

< m.

(ii) There exists some r > 1 such that
N
X
λ6j (S N xN )j
N 2γ(α−1)
j=1

EπN

4

!r
< ∞.

Remark 5.5. Let us clarify the meaning of Assumption 5.1, Assumption 5.2 and Condition
5.3.
• In Assumption 5.1, we do not allow the entries of S N to depend explicitly on N
because this can be achieved by simply changing scaling (i.e., values of α).
• Recall that in stationarity xi,N ∼ N (0, λ2i ) and observe that without any assumption
on S N the addends in the sum
N 1/2 N

h(C )

N

N

z ,S x i =

N
X

λi z

i,N

N

N i

(S x ) =

i=1

N
X

λi z

i,N

N
X

N j,N
Si,j
x

(5.18)

j=1

i=1

are non-independent and non-identically distributed random variables, so the Central Limit theorem does not necessarily hold and in general the behaviour of the
above sum will depend on S N . However, if we want QN
α to be order one, the validity
of the CLT is a reasonable requirement to impose on S N .
• Condition 5.3 allows simplifications in the calculation of the asymptotic variance for
terms of the form (5.18). Indeed a generalized form of the CLT (see [8, Theorem
2.1]) implies that (under suitable assumptions) the asymptotic variance of the sum
(5.18) is related to the following sum of conditional expectations:
!
k−1
X
X

2
a2k,N , where a2k,N := E λ2k z k,N (S N xN )k
z i,N (S N xN )i .
i=1

k

This may not be easy to calculate in general. If Condition 5.3 holds then clearly
the above simplifies and one can apply the Central Limit Theorem in its simplest
form to the summation (5.18). In particular, if Condition 5.3 holds then each term
is independent of the sum of the previous terms so we obtain that for large N (and,
again, assuming that we are in stationarity)
1
N (α−1)γ

h(C N )1/2 z N , S N xN iC N ∼ N (0, c1 ).

Notice that we are requiring that also the structure of (S N )2 and (S N )3 allows for
a similar simplification. So the main effect of Condition 5.3 is to obtain ci = di
for i ∈ {1, 2, 3}. A matrix that satisfies all the requirements of Condition 5.3 is the
N
matrix S N = [SijN ]N
is defined by (4.7) with Ji = 1
i,j=1 of (4.7). If, for example, S
N 2
N 3
N
N 4
then (S ) = −Id so (S ) = −S and (S ) = Id.
• Roughly speaking, while Assumption 5.1 imposes the validity of a CLT, Assumption
5.2 requires that a Law of Large numbers should hold. This assumption is again
made to guarantee that QN
α is of order one (or zero), which is the only case of
N
interest, as we have argued (Because QN
α ' Rα , this is the same as assuming that
N
Rα is O(1). )
• We only consider values of α in the range α ≥ 1 for consistency: the sequence
EπN kS̃ N xN k2C N is positive (and non identically zero), and the entries of S N do not
17

depend on N . Therefore there is no 0 ≤ α < 1 such that Assumption 5.2 is satisfied
for some finite constant c1 .
• Finally, we make the obvious observation that, under Assumption 5.1, Assumption
5.2 and Condition 5.3 if the constant c1 is assumed to be vanishing, c1 = 0, then
2 2
also c2 = c3 = 0. Therefore in this case QN
If c1 = 0 then we
α tends to zero in L .
N N
rule out the case α = 1 because this would imply EπN kS̃ x kC N → 0, which could
be achieved with a different rescaling. If at least one of these constants is non-zero
then QN
α is O(1).
To recap, we have that for N large, then
 6 6
` `
N
Q̄ ≈ N −
if γ = 1/6,
32 16

Q̄N ≈ 0 if γ > 1/6 .

(5.19)

On the other hand, if Assumption 5.1 and Assumption 5.2 hold, then, from (5.7)-(5.12), we
deduce after some calculus 3 that asymptotically
QN
α ≈ N (−ã, b̃),

γ ≥ 1/6,

(5.20)

where
1
ã := 2`2(α−1) c1 + `2(2α−1) c2
2
2(α−1)
b̃ := 4`
d1 + 5`2(2α−1) d2 + `2(3α−1) d3 .

(5.21)
(5.22)

From this we can heuristically deduce that if Z̃ is a normally distributed random variable
with
 6

`
`6
+ b̃ ,
Z̃ ∼ N − − ã,
32
16
then, for N large one has
N
Q̄N + QN
α = Q ≈ Z̃

if γ = 1/6,

QN ≈ N (−ã, b̃) if γ > 1/6 .

(5.23)

As explained in Remark 5.5 (see third bullet point), if also Condition 5.3 is enforced, then
ci = di for i ∈ {1, 2, 3}; therefore, for N large, one finds that QN is approximately distributed
like Q,
QN ≈ Q,
(5.24)
where
 6

`
`6
Q ∼ N − − a,
+b
if γ = 1/6,
32
16

Q ∼ N (−a, b) if γ > 1/6,

(5.25)

with
1
a := 2`2(α−1) c1 + `2(2α−1) c2
2
2(α−1)
b := 4`
c1 + 5`2(2α−1) c2 + `2(3α−1) c3 .
2This

(5.26)
(5.27)

does not mean that the effect of the irreversible term is destroyed. It simply means that the
acceptance probability will not feel it.
3These calculations are a bit long but straightforward and follow the lines of the calculation done in B.1,
proof of vi
18

Before moving to the heuristic analysis of the drift coefficient, we set


hp := E (1 ∧ eQ ) ,

(5.28)

with Q as in (5.25). The constant hp represents the limiting average acceptance probability
and it can be calculated by using the following lemma.
Lemma 5.6. If G is a normally distributed random variable with G ∼ N (µ, δ 2 ) (µ, δ 6= 0)
then

µ
 µ
2
,
E(1 ∧ eG ) = eµ+δ /2 Φ − − δ + Φ
δ
δ
and
 µ

2
E(G1{G<0} ) = eµ+δ /2 Φ − − δ ,
δ
where Φ is the CDF of a standard Gaussian.
Proof. By direct calculation.



Therefore,
b

hp = E(1 ∧ eQ )=e 2 −a Φ

6

−` /32 − b + a
p
(`6 /16) + b

!

6

+Φ

−` /32 − a
p
`6 /16 + b

!
,

if γ = 1/6, a, b, ∈ R
(5.29)

and




√
a
a
hp = e Φ √ − b + Φ − √ , if γ > 1/6, a, b 6= 0 .
b
b
The above discussion makes it clear that the acceptance probability hp depends on the
constants ci and therefore, ultimately, on the choice of the antisymmetric matrix S N . Notice
that if S N takes the Jordan block diagonal form then c2 = c3 = 0 and b = 2a resulting in
limiting acceptance probability hp = hJp , as in (4.12).
b
−a
2



5.2. Heuristic analysis of the drift coefficient. In view of the discussion in Section 5.1,
we comment separately on the case α > 2, α = 2 and α < 2.
• Let us first assume that γ = 1/6 and suppose α > 2. Then the approximate drift of
the chain is given by (4.3)-(4.4), where we take ζ = 2. We will prove that when α > 2 the
addend (4.4) is asymptotically small. This happens because QN and z N are asymptotically
independent; in particular we will show that the correlations between QN and z N decay
faster than N −1/6 , that is
h

i
1/6
QN
1/2 N
N Ek 1 ∧ e
C zk+1 → 0
(the formal statement of the above heuristic is contained in B.17 and B.18. Therefore



`α N N
`2
N
N
N
1/3
QN
d (xk ) ' N Ek (1 ∧ e ) − 1/3 xk + α/6 S̃ xk
2N
N

 2 


` N
`α
QN
QN
N N
= Ek (1 ∧ e ) − xk
+ Ek (1 ∧ e ) (α−2)/6 S̃ xk .
(5.30)
2
N
19

We now use the definition (5.28) and the heuristic approximation (5.24) and observe that if
α > 2 the second addend in (5.30) disappears in the limit; therefore the drift coefficient of
the limiting diffusion is
`2
dp (x) = − hp x .
2
• If instead γ = 1/6 and α = 2 then the second addend in (5.30) does not disappear in the
limit; on the contrary, with the same reasoning as in the previous case, it will converge to
hp `2 S̃x. Moreover, the term in (4.4) is no longer asymptotically small and contributes to the
limiting drift. The reason why this happens can be seen by recalling that QN = Q̄N + QN
α:
clearly, the correlations between Q̄N and z N are not affected by the choice of α; however the
N
value of α does affect the decay of the correlations between QN
α and z .
In B.2 we present a calculation which shows that (under appropriate conditions on S N ),
if α = 2 then
h

i
N 1/6 `Ek

1 ∧ eQ

N

N
C 1/2 zk+1
' −2`2 νp S̃xk

where νp is a real constant, namely
νp := EeQ 1{Q<0}

if

γ=1/6

=

e

b
−a
2

Φ

−`6 /32 − b + a
p
(`6 /16) + b

!
.4

(5.31)

Thus we find that the drift of the limiting diffusion is
d(x) = −

`2
hp x + τp `2 S̃x,
2

where
τp := hp − 2νp .
By Lemma 5.6 one can also see that, again if γ = 1/6,


`6
− −a
.
τp = −hp + 2Φ  q32
`6
+b
16

(5.32)

(5.33)

If the constant c1 = 0 or if more generally b = 2a (e.g., in the Jordan block diagonal form
case), then τp = 0 and the limiting drift reduces to
`2
hp x.
2
Notice that the constant hp can never vanish in this regime.
• As already observed, if γ ≥ 1/6 and 1 ≤ α < 2 we need to scale the algorithm differently;
in particular, the approximate drift is given this time by





`2 N
`α N N
`
N
αγ
αγ
QN
N
N 1/2 N
QN
dp (xk ) ' N Ek (1 ∧ e ) − 2γ xk + αγ S̃ xk
+ N Ek (1 ∧ e ) γ (C ) zk+1
2N
N
N
(5.34)
d(x) = −

4Having

used Lemma 5.6 for the last equality.
20

Therefore in this case (α < 2) the first term in (5.34) asymptotically disappears. As for the
third addend, this is again asymptotically not small (see Lemma B.2, part (ii) of B) and the
second addend contributes with a term `α hp S̃x. Therefore, the limiting drift is
d(x) = τp `α S̃x,
with τp defined as in (5.32). Again, if b = 2a (which happens e.g. for c1 = c2 = c3 = 0),
then τp = 0. Therefore, in this case the limit of the chain would simply be ẋ = 0.
5.3. Heuristic analysis of the diffusion coefficient. The heuristic argument for approximate diffusion coefficient is relatively straightforward.


N
N
N
−
x
−
E
(x
−
x
)
MkN = N ζγ/2 xN
k
k+1
k
k+1
k
h
i
ζγ/2
N
N
N
N
N
N
=N
β̃ (yk+1 − xk ) − Ek [β̃ (yk+1 − xk )]
 

`2 N
`α N N
` 1/2 N
ζγ/2
N
=N
β̃
− 2γ xk +
S̃ xk + γ C zk+1
2N
2N αγ
N


α
2
`
` 1/2 N
`
N
N N
N
S̃ xk + γ C zk+1 ]
−Ek [β̃
− 2γ xk +
2N
2N αγ
N

 2




` N N
(ζ/2−2)γ
(ζ/2−α)γ
α N N N
(ζ/2−1)γ
N 1/2 N
=N
` β̃ S̃ xk + N
`β̃ C zk+1
− β̃ xk + N
2
− N −ζγ/2 dN (xN
k )

(5.35)

By Lemma A.2, we get that in the case γ = 1/6, α ≥ 2 and ζ = 2, MkN behaves like
N (0, `2 hp C); in the case γ ≥ 1/6, α < 2 and ζ = α and MkN converges to zero.
6. On numerical computation of acceptance probability
In this section we aim to show via concrete examples how to construct sequences of matrices
S that satisfy the assumptions of the paper. In addition we calculate in these examples
the constants di , ci for i = 1, 2, 3 that are necessary in order to numerically approximate the
acceptance probabilities in each case.
Let us first discuss how one can construct S N so that it satisfies the assumptions of the
paper, but at the same time leads to non-zero constant c1 . Let us again consider a Jordan
block-diagonal form for S N of (4.7) but now set S1N = [Sij ]N
i,j=1 with Sij = Ji and Sji = −Ji
for j = i + 1 and i = 1, 3, 5, 7, · · · , N − 1 and Sij = 0 otherwise with Ji to be chosen.
Let us set λj = j −k for k > 1/2. Then, we can compute
N

EπN kS̃1 xN k2C =

N
X

λ2i

i=1

=2

N
X

N
X
j=1

(Sij )2 λ2j =

X

i−2k |Si,i+1 |2 (i + 1)−2k +

X

i−2k |Si,i−1 |2 (i − 1)−2k

i even

i odd

(2i − 1)−2k Ji2 (2i)−2k

i=1

In order to make sure that c1 6= 0, let us choose Ji = (2i − 1)k (2i)k i[2(α−1)γ−1]/2 with α > 1.
P
2(α−1)γ−1
Then, we have that EπN kS̃1 xN k2C = 2 N
.
i=1 i
21

We obtain that
kS̃1 xN k2C N
12
c1 = lim EπN 2(α−1)γ
= lim
N →∞
N →∞ N
N

PN

i 2(α−1)γ−1
N
N 2(α−1)γ



i=1

Z

1

=2

x2(α−1)γ−1 dx =

0

1
.
(α − 1)γ

Next we compute c2 for this particular choice of the S matrix. We have
EπN kS̃12 xN k2C N =

N
X

λ2i

i=1

=2

N
X

((S S̃)ij )2 λ2j =

j=1

N
X

i

−4k

i

−4k

N
X

λ2i ((S S̃)ii )2 λ2i = 2

i=1

(i + 1)

−4k

(i + 1)

−4k

N
X

4
λ4i λ4i+1 J(i+1)/2

i odd
4k

(2(i + 1)/2 − 1) (2(i + 1)/2)

4k



i odd

=2

N
X

4k

(i) (i + 1)

i odd

4k



i+1
2

i+1
2

2[2(α−1)γ−1]

2[2(α−1)γ−1]

2[2(α−1)γ−1]
N 
X
i+1
=2
2
i odd
=2

N
X

i2[2(α−1)γ−1]

i=1

Similarly to the computation for c1 , we obtain that c2 = 0. As far as c3 is concerned, we
obtain by a computation similar to the computation for c2 that
kS̃1T S̃32 xkC
=0
N →∞
N 2(3α−1)γ
For the same reasons as in Remark 5.5 we also have di = ci for i = 1, 2, 3. Hence,
1
we have obtained that in this example d1 = c1 = (α−1)γ
and d2 = c2 = d3 = c3 = 0.
 6

6
`
`
1
Then, this implies that Q ∼ N − 32
− a, 16
+ b in (5.25), with a = 2`2(α−1) (α−1)γ
and
c3 = lim EπN

1
. Hence, in this case b = 2a and thus τp = 0.
b = 4`2(α−1) (α−1)γ
Next we discuss optimal acceptance probabilities. Since in this case b = 2a which then
implies τp = 0, we can answer this question in the case of regimes i), i.e., α > 2 and ii)
α = 2. In these cases we can find the ` that maximizes the speed of the limiting diffusion,
i.e., `2 hp and then for the maximizing ` compute the limiting acceptance probability hp as
given by (5.29). In Table 3 we record optimal acceptance probabilities hp derived in (5.29)
for different choices for α when γ = 1/6 in the case the matrix S being S1 as defined above.

α
optimal hp

2
4
6
8
10
15
30
0.234 0.574 0.702 0.767 0.803 0.848 0.884

Table 3. Optimal acceptance probability when S N = S1N for different α0 s
when γ = 1/6.

In particular, Table 3 shows that by increasing α, we can increase the optimal acceptance
probability as desired.
22

Next, we discuss two cases where c1 = c2 = c3 = 0. In particular, one can define the
matrix S2N = [Sij ]N
i,j=1 such that Sij = 1 and Sji = −1 for j = i + 1 and i = 1, 2, 3, 4, · · · , N
and Sij = 0 otherwise. A similar alternative construction is to set S3N = [Sij ]N
i,j=1 with
Sij = 1 and Sji = −1 for j = i + 1 and i = 1, 3, 5, · · · , N − 1 and Sij = 0 otherwise. In both
cases S2N and S3N one can check that ci = di = 0 for all i = 1, 2, 3. This then implies that
`6 `6
, 16 in (5.25), as a = b = 0 in this case. Hence, one recovers that in this case
Q ∼ N − 32
the optimal acceptance probability for both choices S1N and S2N is achieved at 0.574 as in
the classical work [26].
In Section 7 we present simulation studies based on using S1N , S2N and S3N confirming and
illustrating the theoretical results.
7. Simulation studies
The goal of this section is to illustrate the theoretical findings of this paper via simulation
studies. Upon close inspection of Theorem 4.2 it becomes apparent that closed form expressions of quantities such as optimal acceptance probabilities and optimal scalings are hard to
obtain and dependent on the choice of the sequence S N . For this reason we will resort to
simulations.
As a measure of performance we choose similarly to [26, 27] the expected squared jumping
distance (ESJD) of the algorithm defined as
1
ESJD = E Xk+1
− Xk1

2

It is easy to see that maximizing ESJD is equivalent to minimizing the first order lag
autocorrelation function, ρ1 = Corrπ (Xk+1 , Xk ). Hence, motivated by [27] and based on
Theorem 4.2 we define the computational time as
1
CTp = − ζγ
(7.1)
N log ρ1
where ζ = 2 if α ≥ 2 and ζ = α if α = 2. We seek to minimize CTp . A word of caution
is necessary here. Even though CTp as a measure of performance is easily justifiable in
the case that one has a diffusion limit, i.e., when α ≥ 2, the situation is more complex for
α < 2. By Theorem 4.2 if α < 2, the algorithm converges to a deterministic ODE limit
1
and not to a diffusion limit. Hence, in this case the absolute jumping distance Xk+1
− Xk1
may make more sense as a measure of performance in the limit as N → ∞. However, in
the deterministic case, maximizing the absolute jumping distance or the squared jumping
distance are equivalent tasks. Hence, working with the acceptance probability and with
ESJD and as a consequence with CTp as defined in (7.1) may be reasonable criterion on
quantifying convergence.
The simulations presented below were performed using a MPI C parallel code. The numerical results that follow report sample averages and sample standard deviations of quantities
such as acceptance probability and different observables. The reported values were obtained
based on 1600 independent repetitions of MALA and ipMALA algorithm with 104 steps
for each independent realization. For each independent realization of the algorithm the numerical values selected are the ones corresponding to minimizing CTp . We also remark here
that in all of the simulation studies CTp was a convex function of the acceptance probability
(this resembles the behavior in the standard MALA case, see [27]).
23

Below we report a statistical estimator and its corresponding empirical standard deviation
for the optimal acceptance probability. In addition, to demonstrate convergence to the
equilibrium we also computed a few observables. We report below statistical estimators
along with their corresponding statistical standard deviations for
θ2 = EπN kXk2 , and θ3 = EπN

N
X

Xi3 ,

i=1

where with some abuse of notation we have denoted by Xi the ith component of the X vector.
In regards to the measure π N ∼ N (0, C N ), where C N = diag(λ21 , · · · , λ2N ) with λ2i = i−2 .
Let us first study the behavior of the algorithm when S N = S1N , as defined in Section 6.
1
In this case c1 = (α−1)γ
6= 0 and c2 = c3 = 0. However, as we saw in Table 3 the effect of the
value of α on the optimal acceptance probabilities is significant.
Let us demonstrate the behavior via a number of simulation studies reported in Tables
4-6. As we described in Section 6, such a computation makes sense mainly in the case α ≥ 2
and γ = 1/6. Of course one can compute via simulation empirical values for any value of
α > 1 and γ > 0. The acceptance probabilities in the case α ∈ (1, 2) were very small,
leading to estimators with large variances. Hence, we only report values for larger values of
α where one may expect to have results that are comparable to what standard MALA gives.
For completeness, we also present data for γ = 1/2, even though for standard MALA the
optimal choice is γ = 1/6.
Notice that the estimates for acceptance probabilities b
hp align very well with the theoretical predictions that appear in Table 3. In particular, as α increases the optimal acceptance
probability also increases, in accordance to what the theory predicts.
γ = 1/6
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30
γ = 1/2
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30

\
(b
hp , sd(h
p ))
(0.541, 0.043)
(0.229, 0.056)
(0.566, 0.065)
(0.686, 0.063)
(0.751, 0.051)
(0.784, 0.052)
(0.821, 0.045)
(0.855, 0.026)

\
(θb2 , sd(θ
2 ))
(1.073, 0.022)
(0.935, 0.081)
(1.026, 0.034)
(1.045, 0.072)
(1.051, 0.026)
(1.056, 0.061)
(1.059, 0.031)
(1.069, 0.025)

\
(θb3 , sd(θ
3 ))
(−0.001, 0.061)
(−0.001, 0.211)
(0.006, 0.141)
(0.003, 0.101)
(0.0001, 0.101)
(0.004, 0.172)
(−0.002, 0.096)
(−0.003, 0.087)

dp
CT
0.725
8.139
3.367
1.946
1.418
1.151
0.951
0.798

(0.753, 0.008)
(0.205, 0.066)
(0.555, 0.092)
(0.691, 0.083)
(0.761, 0.068)
(0.808, 0.064)
(0.866, 0.054)
(0.924, 0.046)

(1.072, 0.020)
(0.869, 0.118)
(0.979, 0.087)
(0.997, 0.056)
(1.001, 0.038)
(1.015, 0.036)
(1.023, 0.032)
(1.038, 0.067)

(−0.002, 0.068)
(0.0002, 0.339)
(−0.012, 0.242)
(0.002, 0.174)
(−0.008, 0.166)
(0.002, 0.168)
(−0.0003, 0.156)
(−0.002, 0.206)

0.198
9.419
2.022
1.411
1.232
1.051
0.551
0.213

Table 4. Dimension N = 10 and S N = S1N .

To visualize the situation, in Figure 1 we have plotted α versus the empirical acceptance
probabilities b
hp for γ = 1/6. As expected by the theory, as α increases, b
hp increases.
24

γ = 1/6
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30
γ = 1/2
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30

\
(b
hp , sd(h
p ))
(0.538, 0.052)
(0.216, 0.077)
(0.563, 0.074)
(0.689, 0.068)
(0.748, 0.064)
(0.784, 0.063)
(0.824, 0.045)
(0.849, 0.035)

\
(θb2 , sd(θ
2 ))
(1.054, 0.025)
(0.882, 0.177)
(1.004, 0.041)
(1.026, 0.035)
(1.036, 0.042)
(1.048, 0.098)
(1.048, 0.031)
(1.054, 0.025)

\
(θb3 , sd(θ
3 ))
(0.004, 0.089)
(−0.035, 0.451)
(0.004, 0.180)
(−0.009, 0.141)
(0.003, 0.157)
(0.006, 0.181)
(−0.005, 0.123)
(−0.0002, 0.114)

dp
CT
0.912
8.453
3.143
2.024
1.512
1.269
0.991
0.843

(0.849, 0.002)
(0.222, 0.102)
(0.564, 0.013)
(0.699, 0.023)
(0.774, 0.071)
(0.816, 0.034)
(0.873, 0.088)
(0.932, 0.073)

(1.038, 0.037)
(0.681, 0.111)
(0.837, 0.107)
(0.887, 0.139)
(1.001, 0.019)
(0.967, 0.056)
(0.948, 0.101)
(0.958, 0.071)

(0.0004, 0.147)
(0.025, 0.437)
(0.008, 0.378)
(0.003, 0.272)
(0.005, 0.267)
(−0.002, 0.263)
(0.008, 0.264)
(−0.001, 0.114)

0.211
3.396
1.696
1.405
1.269
1.101
0.781
0.583

Table 5. Dimension N = 50 and S N = S1N .

γ = 1/6
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30
γ = 1/2
standard MALA
α=2
α=4
α=6
α=8
α = 10
α = 15
α = 30

\
(b
hp , sd(h
p ))
(0.538, 0.057)
(0.225, 0.071)
(0.559, 0.083)
(0.688, 0.073)
(0.749, 0.068)
(0.784, 0.059)
(0.823, 0.051)
(0.851, 0.034)

\
(θb2 , sd(θ
2 ))
(1.047, 0.036)
(0.891, 0.178)
(1.003, 0.039)
(1.019, 0.047)
(1.032, 0.036)
(1.034, 0.038)
(1.043, 0.038)
(1.053, 0.033)

\
(θb3 , sd(θ
3 ))
(0.001, 0.093)
(0.003, 0.352)
(−0.008, 0.180)
(−0.001, 0.168)
(−0.006, 0.155)
(−0.001, 0.144)
(−0.008, 0.131)
(−0.0016, 0.124)

dp
CT
0.938
7.654
3.467
2.012
1.488
1.245
0.981
0.942

(0.948, 0.002)
(0.221, 0.101)
(0.562, 0.015)
(0.696, 0.034)
(0.771, 0.082)
(0.817, 0.067)
(0.873, 0.093)
(0.937, 0.075)

(1.019, 0.045)
(0.782, 0.113)
(0.843, 0.139)
(0.841, 0.121)
(0.951, 0.082)
(0.989, 0.063)
(1.003, 0.102)
(0.908, 0.101)

(0.0001, 0.183)
(0.015, 0.357)
(0.012, 0.313)
(0.001, 0.275)
(−0.001, 0.217)
(−0.003, 0.311)
(−0.002, 0.235)
(−0.013, 0.111)

0.211
3.145
1.551
1.299
1.187
1.127
0.967
0.486

Table 6. Dimension N = 100 and S N = S1N .

dp . In each case, standard MALA is
In Figure 2 we have plotted α versus the statistic CT
represented by the rightmost value for α. As expected by the theory, in each of the cases,
dp converges to the corresponding value of standard MALA.
as α increases CT
Next we record simulations for the antisymmetric matrix S N = S2N , defined in Section 6.
We remark here that we also did the same simulations with the S = S3N and the numerical
results were statistically the same. So, we only report the ones for S N = S2N , see Tables 7-9.
In the case of S N = S2N , we have that ci = di = 0 and thus the optimal acceptance
probability in the limit as N → ∞ is the same as with standard MALA, i.e. around 0.574.
However, as the Tables 7-9 demonstrate when γ = 1/6 (the optimal choice for standard
MALA) this is being realized for much larger values of N for ipMALA compared to MALA.
25

0.0

0.2

0.4

h_p

0.6

0.8

1.0

Acceptance probabilities for different values of alpha
gamma=1/6, S=S1 and N=10
gamma=1/6, S=S1 and N=50
gamma=1/6, S=S1 and N=100

5

10

15

20

25

30

alpha

Figure 1. Empirical acceptance probability b
hp versus values of α for S N =
N
S1 and γ = 1/6. For reference the limiting optimal acceptance probability for
standard MALA is 0.574.

10

Computational time as defined via (7.1) for different values of alpha

2

4

CT_p

6

8

gamma=1/2, S=S1 and N=10
gamma=1/2, S=S1 and N=50
gamma=1/2, S=S1 and N=100

0

0

2

4

CT_p

6

8

10

Computational time as defined via (7.1) for different values of alpha
gamma=1/6, S=S1 and N=10
gamma=1/6, S=S1 and N=50
gamma=1/6, S=S1 and N=100

10

20

30

40

50

10

alpha

20

30

40

50

alpha

(b) S N = S1N and γ = 1/2

(a) S N = S1N and γ = 1/6

dp versus values of α. In each of the figures the rightFigure 2. Measure CT
most value for α corresponds to standard MALA

On the other hand if we take γ > 1/6 (which is not optimal for standard MALA), e.g.,
γ = 1/2 in Table 9, and focus on the fluid regime, regime iii), we notice that CTp is
estimated to be smaller with ipMALA compared to MALA, as predicted by our theory.
dp . In
To visualize the situation, in Figure 3 we have plotted α versus the statistic CT
each case, standard MALA is represented by the rightmost value for α. As we observe from
Tables 7-9 and from Figure 3 and also indicated by the theory, when γ > 1/6 (in particular
dp is smaller when 1 < α < 2. In particular, notice that for γ = 1/2 > 1/6,
γ = 1/2 here), CT
dp seems to obtain a minimum value around α ≈ 1.8 < 2, which is in the fluid regime.
CT
26

γ = 1/6
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0
γ = 1/2
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0

\
(b
hp , sd(h
p ))
(0.537, 0.043)
(0.464, 0.039)
(0.472, 0.041)
(0.481, 0.041)
(0.488, 0.040)
(0.507, 0.043)
(0.535, 0.044)

\
(θb2 , sd(θ
2 ))
(1.071, 0.022)
(1.064, 0.025)
(1.063, 0.025)
(1.061, 0.025)
(1.061, 0.025)
(1.059, 0.024)
(1.061, 0.023)

\
(θb3 , sd(θ
3 ))
(−0.0002, 0.062)
(0.004, 0.071)
(−0.0016, 0.073)
(−0.0015, 0.074)
(0.0005, 0.077)
(−0.0041, 0.076)
(0.00068, 0.077)

dp
CT
0.725
1.289
1.175
1.070
1.005
1.038
1.063

(0.753, 0.008)
(0.646, 0.008)
(0.649, 0.008)
(0.653, 0.009)
(0.655, 0.009)
(0.662, 0.011)
(0.668, 0.012)

(1.071, 0.020)
(1.065, 0.024)
(1.066, 0.024)
(1.063, 0.024)
(1.064, 0.024)
(1.063, 0.023)
(1.064, 0.023)

(−0.0006, 0.066)
(−0.0025, 0.082)
(0.0025, 0.082)
(−0.0011, 0.081)
(0.0010, 0.083)
(0.0011, 0.084)
(0.00317, 0.078)

0.427
0.663
0.466
0.327
0.556
0.549
0.541

Table 7. Dimension N = 10 and S N = S2N .

γ = 1/6
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0
γ = 1/2
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0

\
(b
hp , sd(h
p ))
(0.538, 0.052)
(0.463, 0.047)
(0.476, 0.048)
(0.489, 0.049)
(0.496, 0.052)
(0.518, 0.051)
(0.528, 0.053)

\
(θb2 , sd(θ
2 ))
(1.054, 0.025)
(1.046, 0.028)
(1.046, 0.027)
(1.046, 0.027)
(1.046, 0.031)
(1.046, 0.026)
(1.047, 0.028)

\
(θb3 , sd(θ
3 ))
(−0.0035, 0.086)
(−0.0023, 0.100)
(−0.0021, 0.097)
(0.0016, 0.096)
(0.0017, 0.094)
(−0.0028, 0.096)
(0.00031, 0.093)

dp
CT
0.900
1.826
1.485
1.209
1.053
1.037
1.021

(0.949, 0.003)
(0.788, 0.006)
(0.832, 0.005)
(0.865, 0.005)
(0.883, 0.005)
(0.914, 0.004)
(0.931, 0.003)

(1.036, 0.038)
(1.021, 0.043)
(1.023, 0.041)
(1.028, 0.041)
(1.030, 0.041)
(1.031, 0.038)
(1.034, 0.037)

(−0.0015, 0.144)
(−0.0007, 0.169)
(0.0008, 0.166)
(−0.0043, 0.158)
(−0.0021, 0.158)
(−0.0001, 0.153)
(0.0002, 0.151)

0.776
1.326
0.684
0.361
0.873
0.828
0.803

Table 8. Dimension N = 50 and S N = S2N .
The numerical studies presented above illustrate the theoretical findings of the paper and
also show that the choice of the irreversible perturbations S can have serious consequences
on the performance of the algorithm.
In particular, in the case of S = S1N
• the theoretical acceptance probability for ipMALA is increasing as a function of α.
This means that by increasing α one can increase the optimal acceptance probability
considerably and combined with the fact that (a) the limit of the interpolation of
the chain is a diffusion as in the standard MALA, and that (b) for large values of
dp are the same for ipMALA and standard MALA, we have an
α the values of CT
algorithm with the same cost but much higher acceptance probabilities.
dp for ipMALA converges to that for the standard MALA case as
• the statistic CT
α increases. However, when 1 < α < 2 the performance seems to be worse than
standard MALA due to very low optimal acceptance probabilities in this case.
27

γ = 1/6
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0
γ = 1/2
standard MALA
α = 1.2
α = 1.5
α = 1.8
α = 2.0
α = 2.5
α = 3.0

\
(b
hp , sd(h
p ))
(0.549, 0.054)
(0.473, 0.049)
(0.484, 0.054)
(0.498, 0.054)
(0.507, 0.056)
(0.523, 0.054)
(0.533, 0.057)

\
(θb2 , sd(θ
2 ))
(1.045, 0.027)
(1.037, 0.029)
(1.036, 0.034)
(1.037, 0.036)
(1.041, 0.084)
(1.041, 0.027)
(1.039, 0.034)

\
(θb3 , sd(θ
3 ))
(0.0015, 0.102)
(−0.0012, 0.115)
(−0.0028, 0.111)
(−0.0021, 0.109)
(0.0074, 0.235)
(−0.0018, 0.106)
(0.00019, 0.104)

dp
CT
0.955
2.080
1.623
1.268
1.076
1.052
1.031

(0.975, 0.002)
(0.809, 0.006)
(0.865, 0.005)
(0.904, 0.004)
(0.922, 0.004)
(0.952, 0.003)
(0.966, 0.002)

(1.006, 0.050)
(0.981, 0.055)
(0.989, 0.053)
(0.993, 0.054)
(0.997, 0.049)
(1.002, 0.049)
(1.002, 0.049)

(0.0103, 0.206)
(0.0039, 0.228)
(0.0138, 0.216)
(−0.0066, 0.213)
(0.0004, 0.215)
(0.0018, 0.211)
(−0.0016, 0.204)

0.970
1.705
0.785
0.371
1.055
1.009
0.985

Table 9. Dimension N = 100 and S N = S2N .

2.0

Computational time as defined via (7.1) for different values of alpha
gamma=1/2, S=S2 and N=10
gamma=1/2, S=S2 and N=50
gamma=1/2, S=S2 and N=100

1.0

CT_p

1.5

0.0

0.0

0.5

0.5

1.0

CT_p

2.0

1.5

2.5

3.0

Computational time as defined via (7.1) for different values of alpha
gamma=1/6, S=S2 and N=10
gamma=1/6, S=S2 and N=50
gamma=1/6, S=S2 and N=100

2

3

4

5

2

alpha

3

4

5

alpha

(a) S N = S2N and γ = 1/6

(b) S N = S2N and γ = 1/2

dp versus values of α. In each of the figures the rightFigure 3. Measure CT
most value for α corresponds to standard MALA
• for large N , the estimated standard deviation of the statistics that were computed,
θ2 and θ3 , appear to be slightly larger than the corresponding estimated standard
deviation of the statistics for the standard MALA case, but certainly of the same
order.
On the other hand, in the case of S = S2N
• the theoretical acceptance probability for ipMALA is the same as that of standard
MALA, i.e., around 0.574, even though this seems to be realizable in practice in very
high dimensions as even for N = 100, the estimated optimal acceptance probabilities
were lower.
dp has lower values for γ = 1/2 > 1/6 and 1 < α < 2 when compared
• the statistic CT
to the corresponding values for the standard MALA case. This can be explained by
28

the cost of the algorithm which in that case is N αγ with αγ < 1/3 as opposed to
the N 1/3 which is the cost of the standard MALA case. As α increases, the value of
dp seems to approach that of standard MALA.
CT
• the estimated standard deviation of the statistics that were computed, θ2 and θ3 ,
are either the same or slightly larger than the corresponding estimated standard
deviation of the statistics for the standard MALA case, but again of the same order.
8. Acknowledgments
K.S. was partially supported by NSF CAREER AWARD nsf-dms 1550918.
Appendix A. Proof of Theorem 4.2
In this section we present the proof of our main results. The proof is based on diffusion
approximation techniques analogous to those used in [18]. In [18] the authors consider the
MALA algorithm with reversible proposal. That is, if we fix S = 0 in our paper and Ψ = 0 in
their paper, the algorithms we consider coincide. For this reason we try to adopt a notation
as similar as possible to the one used in [18] and, for the sake of brevity, we detail only the
parts of the proof that differ from the work [18] and just sketch the rest.
We start by recalling that, by Fernique’s theorem,
EπN kxN kp = Ek(C N )1/2 z N kp . 1,

for all p ≥ 1.

(A.1)

This fact will be often implicitly used without mention in the remainder of the paper.
We also recall that the chain {xN
k }k that we consider is defined in (3.6); the drift-martingale
decomposition of such a chain is given in equation (4.1) . Let us start by recalling the
definition of the continuous interpolant of the chain, equations (4.5)-(4.6), and by introducing
the piecewise constant interpolant of the chain xN
k , that is,
x̄(N ) (t) = xN
k

tk ≤ t < tk+1 ,

where tk = k/N ζγ . It is easy to see (see e.g. [21, Appendix A]) that
Z t
(N )
N
(N )
x (t) = x0 +
dN
(v))dv + wpN (t),
p (x̄

(A.2)

0

where
wpN (t)

1

:=

N ζγ/2

k−1
X

MjN + N ζγ/2 (t − tk )MkN .

j=1

For any t ∈ [0, T ], we set
ŵpN (t)

Z
:=
Z




(N )
(v)) − dp (x(N ) (v)) dv + wpN (t)
dN
p (x̄

0
t

=
Z

t

 N (N )

dp (x̄ (v)) − dp (x̄(N ) (v)) dv

0
t

+




dp (x̄(N ) (v)) − dp (x(N ) (v)) dv + wpN (t).

0

29

(A.3)

With the above notation, we can then rewrite (A.2) as
Z t
(N )
N
dp (x(N ) (v))dv + ŵpN (t).
x (t) = x0 +

(A.4)

0

Let now C([0, T ]; H) denote the space of H-valued continuous functions, endowed with the
uniform topology and consider the map
I : H × C([0, T ]; H) −→ C([0, T ]; H)
(x0 , η(t)) −→ x(t).
That is, I is the map that to every (x0 , η(t)) ∈ H × C([0, T ]; H) associates the (unique
solution) of the equation
Z t
x(t) = x0 +
dp (x(s))ds + η(t).
(A.5)
0
N
From (A.4) it is clear that x(N ) = I(xN
0 , ŵp ). Notice that, under our continuity assumption
on S̃, I is a continuous map. Therefore, in order to prove that x(N ) (t) converges weakly to
x(t) (where x(t) is the solution of equation (A.5) with η(t) = Dp W C (t)), by the continuous
mapping theorem we only need to prove that ŵpN converges weakly to Dp W C (t), where W C (t)
is a H-valued C-Brownian motion. The weak convergence of ŵpN to Dp W C (t) is a consequence
of (A.3) and of Lemmata A.1 and A.2. Then, we get the statement of Theorem 4.2. The
proof of Lemma A.1 and Lemma A.2 is contained in the remained of this Appendix.

Lemma A.1. Under Assumption 5.1, Assumption 5.2 and Condition 5.3, the following holds
Z T
(N )
EπN
kdN
(t)) − dp (x̄(N ) (t))kdt −→ 0, as N → ∞
(A.6)
p (x̄
0

and
Z
EπN

T

kdp (x̄N (t)) − dp (x(N ) (t))kdt −→ 0,

as N → ∞,

(A.7)

0

where the function dp (x) has been defined in the statement of Theorem 4.2.
Set now
N
p (x) := Ex

h

1 ∧ eQ

N



1/2

CN z N

i

(A.8)

and


QN
.
hN
(x)
:=
E
1
∧
e
x
p
While hp (see (5.28)) is the limiting average acceptance probability, hN
p (x) is the local average
acceptance probability. The above notation will be used in the proof of the next lemma.
Lemma A.2. If Assumption 5.1, Assumption 5.2 and Condition 5.3 hold, then wpN (t) converges weakly in C([0, T ]; H) to Dp W C (t), where W C (t) is a H-valued C-Brownian motion
and the constant Dp has been defined in the statement of Theorem 4.2.
30

Proof of Lemma A.1. We start by proving (A.7), which is simpler. The drift coefficient dp
is globally Lipshitz; therefore, using (3.6), (3.3) and (3.4), if tk ≤ t < tk+1 , we have
EπN kdp (x̄(N ) (t)) − dp (x(N ) (t))k . EπN kx̄(N ) (t) − x(N ) (t)k
N
N
N
. (N ζγ t − k) EπN kxN
k+1 − xk k . Eπ N kyk+1 − xk k


1
1
1
N
.
+ αγ EπN kxN
Ek(C N )1/2 zk+1
k → 0.
k k+
2γ
N
N
Nγ

Let us now come to the proof of (A.6). From (4.3)-(4.4), we have
N
N
N
dN
p (x) − dp (x) = A1 + A2 + A3 − dp (x),

where


`2 N
:= N Ex (1 ∧ e ) − 2γ x
2N
h
i
QN
N N
N
(ζ−α)γ α
A2 := N
` Ex (1 ∧ e )S̃ x
h
i
1/2 N (B.17)
(ζ−1)γ
QN
AN
:=
N
`E
(1
∧
e
)C
z
= N (ζ−1)γ `N
x
3
p (x).
N
AN
1



ζγ

QN

(A.9)

(A.10)

We split the function dp (x) in three (corresponding) parts:
dp (x) = A1 (x) + A2 (x) + A3 (x),
with

A1 :=

2

− `2 hp x if γ = 1/6 and α ≥ 2
0
otherwise



hp `α S̃x if γ ≥ 1/6 and α ≤ 2
A2 :=
0
otherwise

−2νp `α S̃x if γ ≥ 1/6 and 1 ≤ α ≤ 2
A3 :=
0
otherwise
We therefore need to consecutively estimate the above three terms.
• AN
1 − A1 : if α ≥ 2 (and γ = 1/6) we fix ζ = 2 and we have
h


i
QN
N
E
1
∧
e
−
h
kxk
+ EπN kxN − xk
EπN kAN
−
A
k
.
E
1
x
p
π
1
(A.1)

.


EπN

hN
p (x)

− hp


2 1/2

+ EπN kxN − xk −→ 0,

(A.11)

as the first addend tends to zero by Lemma B.2 of B and the second addend tends to zero
by definition (see also [22, equation (4.3)]). If α < 2 then ζ = α and we have
N

(α−2)γ
EπN kAN
EπN kEx (1 ∧ eQ )xN k . N (α−2)γ → 0 .
1 k . N

(A.12)

• AN
2 − A2 : if α ≤ 2 then, recalling (4.5), a calculation analogous to the one in (A.11)
gives the statement. If α > 2 then we can act as in (A.12).
• AN
3 − A3 : by Lemma B.2 of B (and (4.5)), we have
EπN kAN
3 − A3 k → 0 as N → ∞.
This concludes the proof.


31

Proof of Lemma A.2. The calculations here are standard so we only prove it for the case
γ = 1/6 and α > 2. Let us recall the martingale difference given by (5.35). In the case
2
γ = 1/6 and α > 2, we have ζ = 2 and dp (x) = − `2 hp x. Hence, the expression (5.35)
becomes



 

N
(1−α)/6
α N
N
N 1/2 N
N
N
−1/6
N 1
dp (xk ) + N
` β̃ S̃xk + `β̃ C zk+1 − N −1/6 dN
Mk = N
β̃
p (xk )
hp




h
i
N 1/2 N
−(α−1)/6
α N
N
−1/6 1
N
N
N
N
= `β̃ C zk+1 + N
` β̃ S̃xk − N
hp dp (xk ) − β̃ dp (xk )
hp
2
By Lemma A.1, we have that EπN kdN
p (x) − dp (x)k −→ 0 as N → ∞. This implies that
N
2
N −1/3 EπN khp dN
p (x) − β̃ dp (x)k → 0.

At the same time, we also notice that
N −(α−1)/3 `2α EπN kβ̃ N S̃xk2 . N −(α−1)/3 `2α EπN kxk2 −→ 0 if α > 1.


Hence, if we define MN (x) = Ex MkN ⊗ MkN |xN
k = x , we obtain that up to a constant
h
i
EπN Tr(MN (x)) − Ex k`β̃ N C 1/2 zk2 ≤ N −1/3
Then, as in Lemma 4.8 of [22] we obtain that
h
i
2
N 1/2
2
EπN ` hp Tr(C) − Ex k`β̃ C zk → 0, as N → ∞
which then immediately implies that
EπN `2 hp Tr(C) − Tr(MN (x)) → 0, as N → ∞.
The latter result implies that the invariance principle of Proposition 5.1 of [1] holds, which
then imply the statement of the lemma.

Appendix B. Auxiliary estimates
We first decompose QN as follows: let
`6
`3
−a−
hx, C 1/2 z N iC
32
4N 3γ
`α−1
`(2α−1)
N 1/2 N
N N
− 2 (α−1)γ h(C ) z , S x i − (2α−1)γ hS̃ N (C N )1/2 z N , S̃ N xN iC N
N
N
`(3α−1)
− (3α−1)γ hS̃ N (C N )1/2 z N , (S̃ N )2 xN iC N .
N

Z N := −

(B.1)

Then
QN = Z N + eN
? ,

(B.2)

where
`6
`6
−
kxN k2C N
32 32 N 6γ
`2(α−1)
`2(2α−1)
N N 2
+ a − 2 (α−1)2γ kS̃ x kC N −
k(S̃ N )2 xN k2C N
2γ(2α−1)
N
2N
+ iN (x, z) + eN (x, z),

eN
? :=

32

(B.3)

with
N
iN (x, z) := iN
1 (x, z) + i2 (x, z),

having defined

`4
kxN k2C N − kz N k2
4γ
8N

`2α  N N 2
N
N
N 1/2 N 2
i2 (x, z) :=
kS̃ x kC N − kS̃ (C ) z kC N ,
2N 2αγ
iN
1 (x, z) :=

(B.4)
(B.5)

and
`5
`2(α+1)
N
N 1/2 N
N
hx
,
(C
)
z
i
−
kS̃ N xN k2C
C
8N 5γ
4N 2(1+α)γ
`(2α+1)
`3+α
N 1/2 N
N N
h(C
)
z
,
S
x
i
+
hS̃(C N )1/2 z N , S̃ N xN iC N .
−
4N (3+α)γ
2N (2α+1)γ

eN :=

(B.6)

Finally, we set
`2(α+1)
kS̃ N xN k2C N .
(B.7)
4N 2(1+α)γ
That is, ẽN contains only the addends of eN that depend on the noise z.
Furthermore, we split QN (x, z) into the terms that contain z j,N and the terms that don’t,
N
Qj and QN
j,⊥ , respectively; that is
ẽN (x, z) := eN +

N
QN = QN
j + Qj,⊥ ,

where
N
N
N
N
QN
j := ẽ + (i1 )j + (i2 )j + Zj ,

(B.8)

N
N N
N
N
having denoted by (iN
1 )j , (i2 )j and Zj , the part of i1 , i2 and Z , respectively, that depend
on z j,N .

Lemma B.1. Let Assumption 5.1, Assumption 5.2 and Condition 5.3 hold; then,
2

i)

EπN ẽN

ii)

N 1/3 EπN

.

N
X

1
1
+ 4γ ,
2/3
N
N
λ2j (iN
1 )j

2

for all α ≥ 1, γ ≥ 1/6

−→ 0,

(B.9)

as N → ∞

j=1

iii)

EπN iN
2

iv)

1/3

N

2

EπN

.

N
X

1
, for all α ≥ 1, γ ≥ 1/6
N 4γ
λ2j ZjN

2

−→ 0,

as N → ∞, for all α > 2, γ = 1/6

(B.10)

j=1

1
1
+
,
for all α ≥ 1, γ ≥ 1/6
(B.11)
N 2/3 N 4γ
2
(Varx Z N )1/2 − (VarQ)1/2 −→ 0, as N → ∞, for all α ≥ 1, γ ≥ 1/6 .
(B.12)

v)

EπN eN

vi)

EπN

2

.

33

Proof of Lemma B.1. Recall that, under π N , xi,N ∼ λi ρi , where {ρi }i∈N are i.i.d standard
Gaussians. We now consecutively prove all the statements of the lemma.
Proof of i). Notice that
N

N 1/2 N

EπN hx , (C )

z iC N

2

=E

N
X

2
i i,N

ρz

=N.

i=1

Therefore, since γ ≥ 1/6,
EπN

`5
hxN , (C N )1/2 z N iC N
8N 5γ

2

. N −2/3 .

(B.13)

Furthermore, using (5.15) (which follows from point i) of Assumption 5.2) and (5.11), we
have
1
2
EπN h(C N )1/2 z N , S N xN i . N −8γ ;
(B.14)
(3+α)2γ
N
similarly, using (5.16) (which follows from point ii) of Assumption 5.2) and (5.12),
2
1
N 1/2 N
N N
N hS̃(C )
E
z
,
S̃
x
i
. N −4γ .
(B.15)
C
π
(2α+1)2γ
N
Now the first statement of the lemma is a consequence of (B.6)-(B.7) and the above (B.13),
(B.14) and (B.15).
Proof of ii). This is proved in [22], see calculations after [22, equation (4.18)], so we omit
it.
Proof of iii). This estimate follows again from Assumption 5.2, once we observe that if
x ∼ π N then kS̃ N xN k2C N and kS̃ N (C N )1/2 z N k2C N are two independent random variables with
the same distribution. With this observation in place, we have
1/2

EπN

kS̃xN k2C − kS̃CN z N k2C
N 2αγ

2

1/2

= Eπ N

kS̃xN k2C
c1
c1
kS̃CN z N k2C
−
+
−
N 2αγ
N 2γ N 2γ
N 2αγ

. EπN

kS̃xN k2C
c1
−
N 2αγ
N 2γ

2

2

2

=N

−4γ

EπN

kS̃xN k2C
− c1 ,
N 2(α−1)γ

which gives the claim.
Proof of iv). We recall that ZjN has been introduced in (B.8). Using the antisymmetry
of S N and the definition of S̃ N , we have
−hS̃ N (C N )1/2 z N , S̃ N xN iC = h(C N )1/2 z N , S N S̃ N xN i
and
−hS̃ N (C N )1/2 z N , (S̃ N )2 xiC N = h(C N )1/2 z N , S N (S̃ N )2 xN i.
We can therefore write an explicit expression for ZjN :
`3 xj,N z j,N
`α−1
−
2
λj z j,N (S N xN )j
4N 3γ λj
N (α−1)γ
`2α−1
`3α−1
+ (2α−1)γ λj z j,N (S N S̃ N xN )j + (3α−1)γ λj z j,N (S N (S̃ N )2 xN )j .
N
N

ZjN = −

34

(B.16)

For the sake of clarity we stress again that in the above (S N xN ) is an N -dimensional vector
and (S N xN )j is the j-th component of such a vector. Therefore, recalling (2.3), (2.1), (A.1)
and setting γ = 1/6, we have
N
X

λ2j EπN

2
ZjN

j=1

N
1 X 2
.
λ E ρj z j,N
N j=1 j

+

N
X

EπN
N (2α−1)/3

2

+

N
X

EπN
N (α−1)/3

λ4j (S N S̃ N xN )j

λ4j (S N xN )j

2

j=1

2

+

j=1

N
X

Eπ N
N (3α−1)/3

λ4j (S N (S̃ N )2 xN )j

2

j=1

1
1
1
1
+ (α−1)/3 EπN kS̃ N xN k2 + (2α−1)/3 EπN k(S̃ N )2 xN k2 + (3α−1)/3 EπN k(S̃ N )3 xN k2
N
N
N
N
1
1
.
+ (α−1)/3 EπN kxN k2 .
N
N

.

Therefore,
N

1/3

N
X

λ2j EπN ZjN

2

j=1

.

1
1
+ (α−2)/3 −→ 0,
2/3
N
N

when α > 2.

Proof of v). Follows from Assumption 5.2, from statement i) of this lemma and from (B.7).
Proof of vi). From (B.1),
Varx (Z N ) =Ex −
−

`3
`α−1
1/2 N
hx,
C
z
i
−
2
hC 1/2 z N , Sxi
C
4N 3γ
N (α−1)γ

`(2α−1)
`(3α−1)
1/2 N
1/2
h
S̃C
z
,
S̃xi
−
hS̃CN z N , S̃ 2 xiC
C
N
(2α−1)γ
(3α−1)γ
N
N

2

.

Therefore,
`2(α−1)
`6
1/2
2
E
kxk
+
4
Ex kCN Sxk2
x
C
6γ
2(α−1)γ
16N
N
2(2α−1)
`
`2(3α−1)
1/2
1/2
+ 2(2α−1)γ Ex kCN S S̃xk2 + 2(3α−1)γ Ex kCN S S̃ 2 xk2 + Ex rN
N
N
N
where r contains all the cross-products in the expansion of the variance. By direct calculation and using the antisymmetry of S, one finds that most of such cross products vanish
and we have
`(2α+2)
`(4α−2)
Ex rN :=
hSx,
S̃xi
+
4
Ex kS̃ 2 xk2C .
2N 3γ N (2α−1)γ
N 2(2α−1)γ
Observe that
hSx, S̃xi = kS̃xk2C ;
using this fact, Assumption 5.2 implies that the first addend in the above expression for
Ex rN vanishes as N → ∞. The second addend contributes instead to the limiting variance.
Now straightforward calculations give the result.

Varx (Z N ) =

We recall the definitions
N
p (x)

:= Ex

h

QN

1∧e
35



1/2
CN z N

i

(B.17)

and


QN
hN
(x)
:=
E
1
∧
e
.
x
p
Lemma B.2. Suppose that Assumption 5.1, Assumption 5.2 and Condition 5.3 hold. Then
(i) If α > 2 and γ = 1/6,
N →∞

2
N 1/3 EπN kN
p (x)k −→ 0 ;

(B.18)

(ii) if 1 ≤ α ≤ 2 and γ ≥ 1/6 then
N →∞

α−1
EπN kN γ(α−1) N
νp S̃xk2 −→ 0
p (x) + 2`

(B.19)

where the constant νp has been defined in (5.31).
(iii) if α ≥ 1, γ ≥ 1/6 and S N is such that (c1 , c2 , c3 ) 6= (0, 0, 0), then
EπN hN
p (x) − hp

2 N →∞

−→ 0 ;

(B.20)

(iv) finally, if 1 < α < 2, γ > 1/6 and S N is such that c1 = c2 = c3 = 0, then
EπN hN
p (x) − 1

2 N →∞

−→ 0 ,

i.e. the constant hp in (B.20) is equal to one. This means that, as N → ∞, the
acceptance probability tends to one.
Proof of Lemma B.2. • Proof of (i). Acting as in [22, page 2349], we obtain
hN
p (x), ϕj i

2

. λ2j Ex QN
j

2

.

Taking the sum over j on both sides of the above then gives
2
kN
p (x)k .

N
X

λ2j Ex QN
j

2

.

j=1

Therefore, if we show
N

1/3

N
X

λ2j Ex QN
j

2 N →∞

−→ 0,

j=1

(B.18) follows. From (B.8), it is clear that the above is a consequence of Lemma B.1 (in
particular, it follows from (B.9)-(B.10)).
• Proof of (ii). Let us split QN as follows:
QN = RN + eN + iN
2 ,
where eN and iN
2 are defined in (B.6) and (B.5), respectively, while
RN := I N + i1 + B N + H N ,
36

(B.21)

having set
`2(α−1)
`2(2α−1)
`6
N 2
N 2
−
2
−
kx
k
k
S̃x
k
kS̃ 2 xN k2C
C
C
32 N 6γ
N 2γ(α−1)
2N 2γ(2α−1)
`α−1
B N := −2 γ(α−1) hC 1/2 z N , SxN i
N
`3
`(2α−1)
1/2
H N := − 3γ hx, C 1/2 z N iC − γ(2α−1) hS̃CN z N , S̃xN iC
4N
N
`(3α−1)
1/2
− γ(3α−1) hS̃CN z N , S̃ 2 xN iC ,
N
I N := −

(B.22)
(B.23)

(B.24)

γ(α−1) N
p can be therefore expressed as
and iN
1 is defind in (B.4). The j-th component of N
follows:
i
i
h
h
γ(α−1) j,N
γ(α−1)
QN
γ(α−1)
RN
j,N
j,N
+ T0j (B.25)
=N
Ex (1 ∧ e )λj z
N
p = N
Ex (1 ∧ e )λj z

where T0j := hT0 , ϕj i and
T0 := N

γ(α−1)

h
Ex

QN

(1 ∧ e

RN

) − (1 ∧ e


i
1/2 N
) CN z .

(B.26)

We now decompose RN into a component which depends on z j,N , RjN , and a component that
N
does not depend on z j,N , Rj,⊥
:
N
RN = RjN + Rj,⊥
,

with
RjN := (i1 )j + (B N )j + (H N )j ,
N
N
N
having denoted by (iN
and H N , respectively, that
1 )j , (B )j and (H )j the part of i1 , B
depend on z j,N . That is,
`4
2
(i1 )j = − 4γ z j,N ;
8N
as for (H N )j , it suffices to notice that

(B N )j + (H N )j = ZjN ,
and the expression for ZjN is detailed in (B.16) (just set α = 2 in (B.16)). With this notation,
from (B.25), we further write
h
i
h
i
N
N
N
N
N γ(α−1) Ex (1 ∧ eQ )λj z j = N γ(α−1) Ex (1 ∧ e[R −(i1 )j )−Hj ] λj z j + T0j + T1j
(B.27)
where, like before, T1j := hT1 , ϕj i and


i
h
N
N
N
N
1/2
CN z N .
T1 := N γ(α−1) `Ex (1 ∧ eR ) − 1 ∧ eR −(i1 )j −Hj
We recall that the notation Ex denotes expected value given x, where the expectation is
taken over all the sources of noise contained in the integrand. In order to further evaluate
the RHS of (B.27) we calculate the expected value of the integrand with respect to the law
j
j
of z j (we denote such expected value by Ez and use Ez− to denote expectation with respect
to z N \ z j,N ); to this end, we use the following lemma.
37

Lemma B.3. If G is a normally distributed random variable with G ∼ N (0, 1) then




µ
µ+δ 2 /2
δG+µ
= δe
Φ − − |δ| ,
E G 1∧e
|δ|
where Φ is the CDF of a standard Gaussian.
N
and δ = δjB , where
We apply the above lemma with µ = Rj,⊥

δjB

`α−1
:= −2 γ(α−1) λj (SxN )j .
N

We therefore obtain
!
N
h
i
j
R
N
B 2
N
N
N
z
j,⊥
N γ(α−1) λj Ex (1 ∧ e[R −(i1 )j −Hj )] z j = N γ(α−1) λj Ex− δjB eRj,⊥ +(δj ) /2 Φ − B − δjB
δj
=

N
B 2
−2`α−1 λ2j (S N xN )j Ezx eRj,⊥ +(δj ) /2 Φ
N

B 2 /2

= −2`α−1 λ2j Ex (Sx)j eRj,⊥ +(δj )

−

N
Rj,⊥

δjB

!
−

δjB

j
j
N <0} + T2 + T3
1{Rj,⊥

N

= −2`α−1 λ2j Ex (Sx)j eQ̃ 1{Q̃<0} + T2j + T3j + T4j
N

= −2`α−1 λ2j Ex (Sx)j eQ 1{Q<0} + T2j + T3j + T4j + T5j ,
"
N +(δ B )2 /2
Rj,⊥
j

T2j := −2`α−1 λ2j Ex (Sx)j e

Φ −

N
Rj,⊥

δjB

!
− δjB

−Φ −

N
Rj,⊥

!#

δjB

and
"
T3j

N +(δ B )2 /2
j Rj,⊥
j

:= −2`α−1 λ2j Ex (Sx) e

Φ −

N
Rj,⊥

δjB

!

#
N <0}
− 1{Rj,⊥

h N
i
B 2
QN
N <0} − e
T4j := −2`α−1 λ2j Ex (Sx)j eRj,⊥ +(δj ) /2 1{Rj,⊥
1{QN <0}
h N
i
j
α−1 2
j
Q
Q
T5 := −2` λj Ex (Sx) e 1{QN <0} − e 1{Q<0} .

(B.28)

To prove the statement it suffices to show that
EπN

5
X

kTn k2 → 0 as N → ∞

n=0

These calculations are a a bit lengthy, so we gather the proof of the above in Lemma B.4 below. Assuming for the moment that the above is true, the proof is concluded after recognising
that
−2`2 λ2j (Sx)j EeQ 1{Q<0} = −2`2 (S̃ N xN )j νp .
• Proof of (iii). By acting as in the proof of [22, Lemma 4.5 and Corollary 4.6] we see
that (B.20) is a consequence of (B.2), (B.12) and of the following limit:
EπN eN
?

2

−→ 0 as N → ∞.
38

The above follows from the definition (B.3), Lemma B.1, Assumption 5.2 and [22, equation
(4.7)].
• Proof of (iv). One could show this with the same procedure as in (iii). However, in
this case things are easier, indeed we can write
N

EπN Ex (1 ∧ eQ ) − (1 ∧ e0 )

2

≤ EπN QN

2

→ 0.

The above limit follows simply by the assumption that γ > 1/6 and c1 = c2 = c3 = 0.



Lemma B.4. If Assumption 5.1, Assumption 5.2, Assumption 5.4 and Condition 5.3 hold,
then
N
X
2 N →∞
EπN kTn k2 = EπN
Tnj −→ 0, for all n ∈ {0, 1, . . ., 5},
i=1

where Tn =

PN

i=1 hTn , ϕi iϕi

and the terms T0j , . . ., T5j have been introduced in (B.26)- (B.28).

Proof of Lemma B.4. We consecutively prove the above limit for the terms T0j , . . ., T5j .
• Using the Lipshitzianity of the function u → 1 ∧ eu , the result for T0 follows from the
definition of RN , equation (B.21), and Lemma B.1, statements iii) and v). The result for T1
can be obtained similarly.
• Term T2 : we using the the lipshitzianity of the function Φ and observe that the following
holds
N

EπN ec(Rj,⊥ +δ

2 /2)

N

. EπN ecR . 1 for all c > 0.

The above can be obtained with a reasoning similar to the one detailed in [18, page 916 and
(5.20)], using (i) of Assumption 5.4 . Using the above observations and applying the Hoelder
inequality with the exponent r appearing in (ii) of Assumption 5.4, one then gets,
EπN kT2 k2 .

EπN

4
N
X
λ6j |(Sx)j |

!r !1/r
.

N 2γ(α−1)

j=1

Therefore the term T2 goes to zero by Assumption 5.4.
• The term T3 can be treated with calculations completely analogous to those in [18,
Lemma 5.8]. As a result of such calculations, using the fact that the noise z j,N is always
independent on the current position x, and recalling equation (B.16), we obtain that for any
r, q > 1 (to be later appropriately chosen), the following bound holds:


"
!#2/q r 1/r
N


j,N 2
j
X
z
(λj |(Sx) | + 1)
2
4
j 2
N


EπN kT3 k . EπN
λj (Sx)
Ex N γ(α−1)
· 1+
+ Zj


|R | N
+1
N 4γ
j=1

(B.29)
Now set

DN :=

1
Ex

(|RN | N γ(α−1) + 1)2
39

r/q
,

so that
(
EπN kT3 k2 .

EπN DN

( N
X

j 2

λ4j (Sx)

h

λ4j (Sx)j

2

h

Ex λj (Sx)j + 1

2 i1/q

)r )1/r
(I)

j=1

(
+

EπN DN

( N
X

Ex λj (Sx)j + 1

2

ZjN

2

i1/q

)r )1/r
. (II)

j=1

Notice that by the bounded convergence theorem, we have
r/q

1
.
EπN DN .
N γ(α−1)

(B.30)

With this observation it is easy to show that the term (I) tends to zero. It is less easy to
show that (II) tends to zero, so for this term we detail calculations a bit more.

 1/(2r)

 "
#2 1/q 2r

N
X

xj,N z j,N

4
j 2
j
2 1/(2r) 

λj (Sx)
Ex (λj (Sx) + 1)
(II) . (EπN DN )
EπN

3γ


N λj
 j=1


2 1/(2r) 
+ (EπN DN
)
EπN


+

2 1/(2r)
(EπN DN
)

+



"

2
λ4j (Sx)j Ex (λj (Sx)j + 1)


 j=1



N
X














j

λj |(Sx) | z
N γ(α−1)

j

j,N

λj (S S̃x) z

4
j 2
j
EπN

λj (Sx) Ex (λj (Sx) + 1)


N γ(2α−1)

 j=1


2 1/(2r)
(EπN DN
)



N
X

 1/(2r)
#2 1/q 2r







j,N

N
λj (S S̃ 2 x)j z j,N
X

j 2
4
j
EπN
λj (Sx) Ex (λj (Sx) + 1)


N γ(3α−1)

 j=1

 1/(2r)
2 1/q 2r

 







 1/(2r)
2 1/q 2r

 



.





We denote by (II)1 to (II)4 the terms in line 1 to 4 of the above array of equations and the
scond factor in line i we denote by (II)ib , so e.g.
2 1/(2r)
)
((II)1b )1/(2r) ,
(II)1 = (EπN DN

where
(II)1b := EπN



N
X



"

2
λ4j (Sx)j Ex (λj (Sx)j + 1)


 j=1

j,N j,N

x z
N 3γ λj


#2 1/q 2r


.



To streamline the presentation we have written the calculations leading to the above four
addends in a way that it looks like the choice of q should be the same for the four terms
above. However, acting appropriately in the computations that give (B.29), one can see that
the q does not need to be the same for each one of the above addends. We show how to
40

study (II)1 and (II)3 , the other terms can be done with similar tricks. Starting from (II)1 ,
because of (B.30), we just need to prove that (II)1b is bounded. We will do slightly better
in what follows. Recall that by assumption
!
2 p
N
N
j 2p
X
X
λ2p
λ2j |(Sx)j |
j |(Sx) |
EπN
< ∞.
(B.31)
≤ EπN
N 2pγ(α−1)
N 2γ(α−1)
j=1
j=1
Choosing q = 2 in the definition of (II)1b and recalling xj,N ∼ λj ρj , we get
!
3 2r
N
N
j 6r
X
X
λ2j λ3j |(Sx)j |
λ2j λ6r
j |(Sx) |
N
N
(II)1b = Eπ
. Eπ
−→ 0,
N 3γ
N 6rγ
j=1
j=1
where in the last inequality we have used the weighted Jentsen’s inequality (relying on the
fact that {λ2j }j is summable) and the convergence of the RHS to zero follows from (B.31).
The term (II)2b can be dealt with analogously, choosing q = 4 (this time when applying the
2
weighted Jentsen’s inequality one should rely on the fact that the sequence {λ4j |(Sx)j | }j is
summable for every x ∈ H). Finally, to deal with (II)3b , we use the fact that the sequence
{(S̃ 2 x)j }j is, by assumpion, bounded for every x ∈ H. Therefore, choosing q = 2 we have:

2r
j
N
λ
(S
S̃x)
X
j
2

(II)3b = EπN 
λ4j (Sx)j λj (Sx)j
γ(2α−1)
N
j=1
= EπN

3
N
X
λ4j |(Sx)j |
j=1

≤ EπN

N γ(2α−1)

2
N
X
λ2j |(Sx)j |
j=1

N γ(2α−1)

!2r
(S̃ 2 x)j
!2r
λ2j (Sx)j

.

Because 2αγ − γ > 2γ(α − 1), the RHS of the above tends to zero by using (B.31). The
term (II)4b can be dealt with in a completely analogous manner.
• The terms T4 and T5 can be studied similarly to what has been done in [14], see calcui,N
lations from equation (8.31), in particular the terms ei,N

3,k , e5,k .
References
[1] Erich Berger. Asymptotic behaviour of a class of stochastic approximation procedures. Probability theory
and related fields, 71(4):517–552, 1986.
[2] Joris Bierkens. Non-reversible metropolis-hastings. Statistics and Computing, pages 1–16, 2015.
[3] A. Bouchard-Côte, A. Doucet and S.J. Vollmer. The bouncy particle sampler: A non-reversible rejectionfree Markov Chain Monte Carlo method. submitted, 2015.
[4] O. F. Christensen, G. O. Roberts, and J. S. Rosenthal. Scaling limits for the transient phase of local
Metropolis-Hastings algorithms. J. R. Stat. Soc. Ser. B Stat. Methodol., 67(2):253 – 268, 2005.
[5] Persi Diaconis, Susan Holmes, and Radford M Neal. Analysis of a nonreversible Markov chain sampler.
Annals of Applied Probability, pages 726–752, 2000.
[6] AB Duncan, T Lelievre, and GA Pavliotis. Variance reduction using nonreversible langevin samplers.
Journal of Statistical Physics, 163(3):457–491, 2016.
[7] AB Duncan, GA Pavliotis, and K Zygalakis. Nonreversible langevin samplers: Splitting schemes, analysis and implementation. submitted, 2017.
41

[8] Aryeh Dvoretzky et al. Asymptotic normality for sums of dependent random variables. In Proc. 6th
Berkeley Symp. Math. Statist. Probab, volume 2, pages 513–535, 1972.
[9] Alan M Horowitz. A generalized guided Monte Carlo algorithm. Physics Letters B, 268(2):247–252,
1991.
[10] Chii-Ruey Hwang, Shu-Yin Hwang-Ma, and Shuenn-Jyi Sheu. Accelerating diffusions. The Annals of
Applied Probability, 15(2):1433–1444, 2005.
[11] B. Jourdain, T.Lelievre, and B. Miasojedow. Optimal scaling for the transient phase of Metropolis
Hastings algorithms: the longtime behavior. Bernoulli, 20(4):1930 – 1978, 2014.
[12] B. Jourdain, T. Lelievre, and B. Miasojedow. Optimal scaling for the transient phase of the random
walk Metropolis algorithm: the mean field limit. Ann. Appl. Probab., 25(4):2263 – 2300, 2015.
[13] E. P. Bernard, W. Krauth, D. B. Wilson. Event-chain algorithms for hard-sphere systems. Physical
Review E, (2009)
[14] J. Kuntz, Michela Ottobre, and Andrew M Stuart. Diffusion limit for the random walk metropolis
algorithm out of stationarity. arXiv preprint arXiv:1405.4896, 2014.
[15] J. Kuntz, Michela Ottobre, and Andrew M Stuart. Non-stationary phase of the MALA algorithm. arXiv
preprint arXiv:1608.08379, 2016.
[16] Jianfeng Lu and Konstantinos Spiliopoulos. Multiscale integrators for stochastic differential equations
and irreversible langevin samplers. arXiv preprint arXiv:1606.09539, 2016.
[17] Y.–A Ma, T. Chen and E. Fox. A complete recipe for stochastic-gradient MCMC. Advances in Neural
Information Processing systems 28, pages 2899–2907, 2015.
[18] Jonathan C Mattingly, Natesh S Pillai, and Andrew M Stuart. Diffusion limits of the random walk
metropolis algorithm in high dimensions. The Annals of Applied Probability, 22(3):881–930, 2012.
[19] P. Monmarche. Piecewise deterministic simulated annealing, ALEA, Lat. Am. J. Probab. Math. Stat.
13 (1), 357–398 (2016)
[20] M. Ottobre. Markov Chain Monte Carlo and Irreversibility. Reports on Math. Phys. (2016)
[21] Michela Ottobre, Natesh S Pillai, Frank J Pinski, and Andrew M Stuart. A function space hmc algorithm
with second order langevin diffusion limit. Bernoulli, 22(1):60–106, 2016.
[22] Natesh S Pillai, Andrew M Stuart, and Alexandre H Thiéry. Optimal scaling and diffusion limits for
the langevin algorithm in high dimensions. The Annals of Applied Probability, 22(6):2320–2356, 2012.
[23] R. Poncet. Generalized and hybrid Metropolis-Hastings overdamped Langevin algorithms. arXiv
1701.05833
[24] Luc Rey-Bellet and Konstantinos Spiliopoulos. Irreversible langevin samplers and variance reduction: a
large deviations approach. Nonlinearity, 28(7):2081–2103, 2015.
[25] Luc Rey-Bellet and Konstantinos Spiliopoulos. Variance reduction for irreversible langevin samplers and
diffusion on graphs. Electronic Communications in Probability, 20, 2015.
[26] Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling of discrete approximations to langevin
diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268,
1998.
[27] Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling for various metropolis-hastings algorithms.
Statistical science, 16(4):351–367, 2001.
[28] Gareth O Roberts and Richard L Tweedie. Exponential convergence of langevin distributions and their
discrete approximations. Bernoulli, pages 341–363, 1996.
[29] P.J. Rossky, J.D. Doll and H.L. Friedman, Brownian Dynamics as smart Monte Carlo Simulation,
J.Chem. Phys., 69 (10), 4628–4633, 1978.
Department of Mathematics, Heriot Watt University,Edinburgh, EH14 4AS, Scotland
E-mail address: m.ottobre@hw.ac.uk
Department of Statistics, Harvard University, Cambridge, MA, 02138-2901,USA
E-mail address: pillai@fas.harvard.edu
Department of Mathematics and Statistics, Boston University, Boston, MA, 02215, USA
E-mail address: kspiliop@math.bu.edu
42

