High-dimensional Bayesian inference via the Unadjusted
Langevin Algorithm

arXiv:1605.01559v3 [] 21 Feb 2018

Alain Durmus1 and Éric Moulines2
1

CMLA - École normale supérieure Paris-Saclay, CNRS, Université Paris-Saclay,
94235 Cachan, France.
2
Centre de Mathématiques Appliquées, UMR 7641, Ecole Polytechnique, France.
February 22, 2018

Keywords: total variation distance, Langevin diffusion, Markov Chain Monte Carlo,
Metropolis Adjusted Langevin Algorithm, Rate of convergence
AMS subject classification (2010): primary 65C05, 60F05, 62L10; secondary 65C40,
60J05,93E35
Abstract
: We consider in this paper the problem of sampling a high-dimensional probability distribution π having a density w.r.t. theR Lebesgue measure on Rd , known up
to a normalization factor x 7→ π(x) = e−U(x) / Rd e−U(y) dy. Such problem naturally
occurs for example in Bayesian inference and machine learning. Under the assumption that U is continuously differentiable, ∇U is globally Lipschitz and U is strongly
convex, we obtain non-asymptotic bounds for the convergence to stationarity in
Wasserstein distance of order 2 and total variation distance of the sampling method
based on the Euler discretization of the Langevin stochastic differential equation,
for both constant and decreasing step sizes. The dependence on the dimension of
the state space of the obtained bounds is studied to demonstrate the applicability
of this method. The convergence of an appropriately weighted empirical measure
is also investigated and bounds for the mean square error and exponential deviation inequality are reported for functions which are measurable and bounded. An
illustration to Bayesian inference for binary regression is presented.

1

Introduction

There has been recently an increasing interest in Bayesian inference applied to highdimensional models often motivated by machine learning applications. Rather than
1
2

Email: alain.durmus@cmla.ens-cachan.fr
eric.moulines@polytechnique.edu

1

obtaining a point estimate, Bayesian methods attempt to sample the full posterior distribution over the parameters and possibly latent variables which provides a way to
assert uncertainty in the model and prevents from overfitting [28], [37].
The problem can be formulated as follows. RWe aim at sampling a posterior distribution π on Rd , d ≥ 1, with density x 7→ e−U (x) / Rd e−U (y) dy w.r.t. the Lebesgue measure,
where U is continuously differentiable. The Langevin stochastic differential equation
associated with π is defined by:
√
(1)
dYt = −∇U (Yt )dt + 2dBt ,
where (Bt )t≥0 is a d-dimensional Brownian motion defined on the filtered probability
space (Ω, F, (Ft )t≥0 , P), satisfying the usual conditions. Under mild technical conditions,
the Langevin diffusion admits π as its unique invariant distribution.
We study the sampling method based on the Euler-Maruyama discretization of
(1). This scheme defines the (possibly) non-homogeneous, discrete-time Markov chain
(Xk )k≥0 given by
p
Xk+1 = Xk − γk+1 ∇U (Xk ) + 2γk+1 Zk+1 ,
(2)

where (Zk )k≥1 is an i.i.d. sequence of d-dimensional standard Gaussian random variables
and (γk )k≥1 is a sequence of step sizes, which can either be held constant or be chosen
to decrease to 0. This algorithm has been first proposed by [14] and [30] for molecular
dynamics applications. Then it has been popularized in machine learning by [17], [18]
and computational statistics by [28] and [32]. Following [32], in the sequel this method
will be refer to as the unadjusted Langevin algorithm (ULA). When the step sizes are held
constant, under appropriate conditions on U , the homogeneous Markov chain (Xk )k≥0
has a unique stationary distribution πγ , which in most of the cases differs from the
distribution π. It has been proposed in [33] and [32] to use a Metropolis-Hastings step
at each iteration to enforce reversibility w.r.t. π. This algorithm is referred to as the
Metropolis adjusted Langevin algorithm (MALA).
The ULA algorithm has already been studied in depth for constant step sizes in
[35], [32] and [27]. In particular, [35, Theorem 4] gives an asymptotic expansion for the
weakPerror between π and πγ . For sequence of step sizes such that limk→+∞ γk = 0
and ∞
k=1 γk = ∞, weak convergence of the weighted empirical distribution of the ULA
algorithm has been established in [23], [24] and [25].
Contrary to these previously reported works, we focus in this paper on non-asymptotic
results. More precisely, we obtain explicit bounds between the distribution of the n-th
iterate of the Markov chain defined in (2) and the target distribution π in Wasserstein
and total variation distance for nonincreasing step sizes. When the sequence of step sizes
is held constant γk = γ for all k ≥ 0, both fixed horizon (the total computational budget
is fixed and the step size is chosen to minimize the upper bound on the Wasserstein or
total variation distance) and fixed precision (for a fixed target precision, the number of
iterations and the step size are optimized simultaneously to meet this constraint) strategies are studied. In addition,
Pquantitative estimates between π and πγ are obtained.
When limk→+∞ γk = 0 and ∞
k=1 γk = ∞, we show that the marginal distribution of
2

the non-homogeneous Markov chain (Xk )k≥0 converges to the target distribution π and
provide explicit convergence bounds. A special attention is paid to the dependency of
the proposed upper bounds on the dimension of the state space, since we are particularly
interested in the application of these methods to sampling in high-dimension.
Compared to [10] and [11], our contribution is twofold. First, we report sharper
convergence bounds for constant and non-increasing step sizes. In particular, it is shown
that the number of iterations required to reach the target distribution π with a precision
ǫ > 0 in total variation is upper bounded (up to logarithmic factors) by (ǫd)−1 , where
d is the dimension, under appropriate regularity assumption (in [11], the upper bound
was shown to be d−1 ǫ−2 ). We show that our result is optimal (up to logarithmic factors
again) in the case of the standard d-dimensional Gaussian distribution, for which explicit
computation can be done. Besides, we establish computable bound for the mean square
error associated with ULA for measurable bounded functions, which is an important
result in Bayesian statistics to estimate credibility regions. For that purpose, we study
the convergence of the Euler-Maruyama discretization towards its stationary distribution
in total variation using a discrete time version of reflection coupling introduced in [4].
The paper is organized as follows. In Section 2, we study the convergence in the
Wasserstein distance of order 2 of the Euler discretization for constant and decreasing
step sizes. In Section 3, we give non asymptotic bounds in total variation distance
between the Euler discretization and π. This study is completed in Section 4 by nonasymptotic bounds of convergence of the weighted empirical measure applied to bounded
and measurable functions. Our claims are supported in a Bayesian inference for a binary
regression model in Section 5. The proofs are given in Section 6. Finally in Section 8,
some results of independent interest, used in the proofs, on functional autoregressive
models are gathered. Some technical proofs and derivations are postponed and carried
out in a supplementary paper [12].

Notations and conventions
Denote by B(Rd ) the Borel σ-field of Rd , F(Rd ) the set of all Borel measurable functions
on Rd and for f ∈ F(Rd ), kf k∞ = supx∈Rd |f (x)|. For µ a probability measure on
(Rd , B(Rd )) and f ∈ F(Rd ) a µ-integrable function, denote by µ(f ) the integral of f
w.r.t. µ. We say that ζ is a transference plan of µ and ν if it is a probability measure
on (Rd × Rd , B(Rd × Rd )) such that for all measurable set A of Rd , ζ(A × Rd ) = µ(A)
and ζ(Rd × A) = ν(A). We denote by Π(µ, ν) the set of transference plans of µ and ν.
Furthermore, we say that a couple of Rd -random variables (X, Y ) is a coupling of µ and
ν if there exists ζ ∈ Π(µ, ν) such that (X, Y ) are distributed according to ζ. For two
probability measures µ and ν, we define the Wasserstein distance of order p ≥ 1 as
def

Wp (µ, ν) =



inf

Z

1/p
.
kx − yk dζ(x, y)
p

ζ∈Π(µ,ν) Rd ×Rd

By [36, Theorem 4.1], for all µ, ν probability measures on Rd , there exists a transference plan ζ ⋆ ∈ Π(µ, ν) such that for any coupling (X, Y ) distributed according to ζ ⋆ ,
3

Wp (µ, ν) = E[kX − Y kp ]1/p . This kind of transference plan (respectively coupling) will
be called an optimal transference plan (respectively optimal coupling) associated with
Wp . We denote
by Pp (Rd ) the set of probability measures with finite p-moment: for all
R
µ ∈ Pp (Rd ), Rd kxkp dµ(x) < +∞. By [36, Theorem 6.16], Pp (Rd ) equipped with the
Wasserstein distance Wp of order p is a complete separable metric space.
Let f : Rd → R be a Lipschitz function, namely there exists C ≥ 0 such that for all
x, y ∈ Rd , |f (x) − f (y)| ≤ C kx − yk. Then we denote
kf kLip = inf{|f (x) − f (y)| kx − yk−1 | x, y ∈ Rd , x 6= y} .
The Monge-Kantorovich theorem (see [36, Theorem 5.9]) implies that for all µ, ν probabilities measure on Rd ,

Z
Z
d
f (x)dν(x) | f : R → R ; kf kLip ≤ 1 .
f (x)dµ(x) −
W1 (µ, ν) = sup
Rd

Rd

Denote by Fb (Rd ) the set of all bounded Borel measurable functions on Rd . For f ∈
Fb (Rd ) set osc(f ) = supx,y∈Rd |f (x) − f (y)|. For two probability measures µ and ν on
Rd , the total variation distance distance between µ and ν is defined by kµ − νkTV =
supA∈B(Rd ) |µ(A) − ν(A)|. By the Monge-Kantorovich theorem the total variation distance
between µ and ν can be written on the form:
Z
1Dc (x, y)dζ(x, y) ,
kµ − νkTV = inf
ζ∈Π(µ,ν) Rd ×Rd

where D = {(x, y) ∈ Rd × Rd |x = y}. For all x ∈ Rd and M > 0, we denote by B(x, M ),
the ball centered at x of radius M . For a subset A ⊂ Rd , denote by Ac the complementary
of A. Let n ∈ N∗ and M be a n × n-matrix, then denote by M T the transpose of M
and kM k the operator norm associated with M defined by kM k = supkxk=1 kM xk.
Define the Frobenius norm associated with M by kM k2F = Tr(M T M ). Let n, m ∈ N∗
and F : Rn → Rm be a twice continuously differentiable function. Denote by ∇F and
~ the vector
∇2 F the Jacobian and the Hessian of F respectively. Denote also by ∆F
d
m
~ (x) is the vector of R such that for all
Laplacian of F defined by: for all x ∈ R , ∆F
P
~
i ∈ {1, · · · , m}, the i-th component P
of ∆F (x) is Q
equals to dj=1 (∂ 2 Fi /∂x2j )(x). In the
sequel, we take the convention that np = 0 and np = 1 for n, p ∈ N, n < p.

2

Non-asymptotic bounds in Wasserstein distance of order
2 for ULA

Consider the following assumption on the potential U :
H 1. The function U is continuously differentiable on Rd and gradient Lipschitz: there
exists L ≥ 0 such that for all x, y ∈ Rd , k∇U (x) − ∇U (y)k ≤ L kx − yk.

4

Under H1, for all x ∈ Rd by [21, Theorem 2.5, Theorem 2.9 Chapter 5] there exists
a unique strong solution (Yt )t≥0 to (1) with Y0 = x. Denote by (Pt )t≥0 the semi-group
associated with (1). It is well-known that π is its (unique) invariant probability. To get
geometric convergence of (Pt )t≥0 to π in Wasserstein distance of order 2, we make the
following additional assumption on the potential U .
H 2. U is strongly convex, i.e. there exists m > 0 such that for all x, y ∈ Rd ,
U (y) ≥ U (x) + h∇U (x), y − xi + (m/2) kx − yk2 .
Under H2, [29, Theorem 2.1.8] shows that U has a unique minimizer x⋆ ∈ Rd . If in
addition H1 holds, then [29, Theorem 2.1.12, Theorem 2.1.9] show that for all x, y ∈ Rd :
h∇U (y) − ∇U (x), y − xi ≥

1
κ
ky − xk2 +
k∇U (y) − ∇U (x)k2 ,
2
m+L

(3)

where

2mL
.
(4)
m+L
We briefly summarize some background material on the stability and the convergence in
W2 of the overdamped Langevin diffusion under H1 and H2. Most of the statements in
Proposition 1 are known and are recalled here for ease of references; see e.g. [6].
κ=

Proposition 1. Assume H1 and H2.
(i) For all t ≥ 0 and x ∈ Rd ,
Z
ky − x⋆ k2 Pt (x, dy) ≤ kx − x⋆ k2 e−2mt + (d/m)(1 − e−2mt ) .
Rd

(ii) The stationary distribution π satisfies

R

Rd

kx − x⋆ k2 π(dx) ≤ d/m.

(iii) For any x, y ∈ Rd and t > 0, W2 (δx Pt , δy Pt ) ≤ e−mt kx − yk.

(iv) For any x ∈ Rd and t > 0, W2 (δx Pt , π) ≤ e−mt kx − x⋆ k + (d/m)1/2 .

Proof. The proof is given in the supplementary document [12, Section 1.1].

Note that the convergence rate in Proposition 1-(iv) is independent from the dimension. Let (γk )k≥1 be a sequence of positive and non-increasing step sizes and for n, ℓ ∈ N,
denote by
ℓ
X
def
γk ,
Γn = Γ1,n .
(5)
Γn,ℓ =
k=n

For γ > 0, consider the Markov kernel Rγ given for all A ∈ B(Rd ) and x ∈ Rd by
Z


Rγ (x, A) = (4πγ)−d/2 exp −(4γ)−1 ky − x + γ∇U (x)k2 dy .
A

5

(6)

The process (Xk )k≥0 given in (2) is an inhomogeneous Markov chain with respect to the
family of Markov kernels (Rγk )k≥1 . For ℓ, n ∈ N∗ , ℓ ≥ n, define
Qn,ℓ
γ = Rγn · · · Rγℓ ,

Qnγ = Q1,n
γ

(7)

with the convention that for n, ℓ ∈ N, n < ℓ, Qℓ,n
γ is the identity operator.
∗
We first derive a Foster-Lyapunov drift condition for Qn,ℓ
γ , ℓ, n ∈ N , ℓ ≥ n.
Proposition 2. Assume H1 and H2.
(i) Let (γk )k≥1 be a nonincreasing sequence with γ1 ≤ 2/(m + L). Let x⋆ be the unique
minimizer of U . Then for all x ∈ Rd and n, ℓ ∈ N∗ ,
Z
ky − x⋆ k2 Qn,ℓ
γ (x, dy) ≤ ̺n,ℓ (x) ,
Rd

where ̺n,ℓ (x) is given by
̺n,ℓ (x) =

ℓ
Y

(1 − κγk ) kx − x⋆ k2 + 2dκ−1

k=n

(

1 − κ−1

ℓ
Y

i=n

)

(1 − κγi )

,

(8)

and κ is defined in (4).
(ii) For any γ ∈ (0, 2/(m + L)], Rγ has a unique stationary distribution πγ and
Z
kx − x⋆ k2 πγ (dx) ≤ 2dκ−1 .
Rd

Proof. The proof is postponed to [12, Section 1.2].
We now proceed to establish the contraction property of the sequence (Qnγ )n≥1 in
W2 . This result implies the geometric convergence of the sequence (δx Rγn )n≥1 to πγ in
W2 for all x ∈ Rd . Note that the convergence rate is again independent of the dimension.
Proposition 3. Assume H1 and H2. Then,
(i) Let (γk )k≥1 be a nonincreasing sequence with γ1 ≤ 2/(m + L). For all x, y ∈ Rd
and ℓ ≥ n ≥ 1,
n,ℓ
W2 (δx Qn,ℓ
γ , δy Qγ ) ≤

(

ℓ
Y

(1 − κγk ) kx − yk2

k=n

(ii) For any γ ∈ (0, 2/(m + L)), for all x ∈ Rd and n ≥ 1,

)1/2

.

n
o1/2
W2 (δx Rγn , πγ ) ≤ (1 − κγ)n/2 kx − x⋆ k2 + 2κ−1 d
.
6

Proof. The proof is postponed to [12, Section 1.3]
Corollary 4. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 2/(m + L). Then for all Lipschitz functions f : Rd → R and ℓ ≥ n ≥ 1, Qn,ℓ
γ f is a
Qℓ
1/2 kf k
Lipschitz function with kQn,ℓ
(1
−
κγ
)
.
f
k
≤
γ
Lip
Lip
k
k=n
Proof. The proof follows from Proposition 3-(i) using

n,ℓ
n,ℓ
n,ℓ
Qn,ℓ
γ f (y) − Qγ f (z) ≤ kf kLip W2 (δy Qγ , δz Qγ ) .

We now proceed to establish explicit bounds for W2 (δx Qnγ , π), with x ∈ Rd . Since π
is invariant for Pt for all t ≥ 0, it suffices to get some bounds on W2 (δx Qnγ , ν0 PΓn ), with
ν0 ∈ P2 (Rd ) and take ν0 = π. To do so, we construct a coupling between the diffusion
and the linear interpolation of the Euler discretization. An obvious candidate is the
synchronous coupling (Yt , Y t )t≥0 defined for all n ≥ 0 and t ∈ [Γn , Γn+1 ) by
(
√
Rt
Yt = YΓn − Γn ∇U (Ys )ds + 2(Bt − BΓn )
√
(9)
Ȳt = ȲΓn − ∇U (ȲΓn )(t − Γn ) + 2(Bt − BΓn ) ,
with Y0 is distributed according to ν0 , Ȳ0 = x and (Γn )n≥1 is given in (5). Therefore since
for all n ≥ 0, W22 (δx Qnγ , ν0 PΓn ) ≤ E[kYΓn − ȲΓn k2 ], taking ν0 = π, we derive an explicit
bound on the Wasserstein distance between the sequence of distributions (δx Qkγ )k≥0 and
the stationary measure π of the Langevin diffusion (1).
Theorem 5. Assume H1 and H2. Let (γk )k≥1 be a nonincreasing sequence with γ1 ≤
1/(m + L). Then for all x ∈ Rd and n ≥ 1,
n
o
⋆ 2
(2)
W22 (δx Qnγ , π) ≤ u(1)
n (γ) kx − x k + d/m + un (γ) ,

where

u(1)
n (γ)

=2

n
Y

(1 − κγk /2)

(10)

k=1

κ is defined in (4) and
"
#

 Y
n
n
X
 −1
dL2 γi dL2 γi2
2
(2)
2
γi κ + γi
un (γ) = L
2d +
(1 − κγk /2) . (11)
+
m
6
i=1

k=i+1

Proof. Let x ∈ Rd , n ≥ 1 and ζ0 = π ⊗ δx . Let (Yt , Y t )t≥0 with (Y0 , Y 0 ) distributed
according to ζ0 and defined by (9). Byhdefinition of W
i 2 and since for all t ≥ 0, π
2
n
2
is invariant for Pt , W2 (µ0 Q , π) ≤ Eζ0 YΓn − Y Γn
. Lemma 20 with ǫ = κ/4,
Proposition 1-(ii) imply, using a straightforward induction, that for all n ≥ 0
Z
i
h
2
(1)
ky − xk2 π(dy) + An (γ) ,
(12)
≤ un (γ)
Eζ0 YΓn − Y Γn
Rd

7

(1)

where (un (γ))n≥1 is given by (10), and
def

An (γ) = L

2

n
X
i=1

γi2



−1

κ

+ γi (2d +

dL2 γi2 /6)

k=i+1

+ L4

n
X
i=1

with

n
Y

(1 − κγk /2)


δ̃i γi3 κ−1 + γi

n
Y

(1 − κγk /2) (13)

k=i+1

i
h
δ̃i = e−2mΓi−1 Eζ0 kY0 − x⋆ k2 + (1 − e−2mΓi−1 )(d/m) ≤ d/m .

Since Y0 is distributed according to π, Proposition 1-(ii) shows that for all i ∈ {1, · · · , n},
δ̃i ≤ d/m .

2

(14)
x⋆ k2 +kx⋆

2

In addition since
for all y ∈ Rd , kx − yk ≤ 2(kx −
− yk ), using Proposition 1R
(ii), we get Rd ky − xk2 π(dy) ≤ kx − x⋆ k2 + d/m. Plugging this result, (14) and (13)
in (12) completes the proof.
Corollary 6. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 1/(m + L). Assume that limk→∞ γk = 0 and limn→+∞ Γn = +∞. Then for all
x ∈ Rd , limn→∞ W2 (δx Qnγ , π) = 0.
We preface the proof by a technical lemma which is established in the supplementary
document Section 2.1.
Lemma 7. Let (γk )k≥1 be a sequence of nonincreasing real numbers, ̟ > 0 and γ1 <
̟−1 . Then for all n ≥ 0, j ≥ 1 and ℓ ∈ {1, . . . , n + 1},
n+1
Y
X n+1

i=1 k=i+1

(1 −

̟γk ) γij

≤

n+1
Y
k=ℓ

(1 − ̟γk )

ℓ−1
X

γij +

i=1

(1)

γℓj−1
.
̟
(2)

Proof of Corollary 6. By Theorem 5, it suffices to show that un (γ) and un (γ), defined
by (10) and (11) respectively, goes to 0 as n → +∞. Using the bound 1 + t ≤ et
(1)
for t ∈ R, and limn→+∞ Γn = +∞, we have limn→+∞ un (γ) = 0. Since (γk )k≥0
(2)
is nonincreasing,
Pn Qnnote that to show2 that limn→+∞ un (γ) = 0, it suffices to prove
limn→+∞ i=1 k=i+1 (1 − κγk /2)γi = 0. But since (γk )k≥1 is nonincreasing, there
exists c ≥ 0 such that cΓn ≤ n − 1 and by Lemma 7 applied with ℓ = ⌊cΓn ⌋ the integer
part of cΓn :
n
n
Y
X

i=1 k=i+1

(1 −

κγk /2) γi2

−1

≤ 2κ

−1

γ⌊cΓn ⌋ + exp −2

⌊cΓn ⌋−1
 X
γi .
κ(Γn − Γ⌊cΓn ⌋ )

(15)

i=1

Since limk→+∞ γk = 0, by the Cesáro theorem, limn→+∞ n−1 Γn = 0. Therefore since
limn→+∞ Γn = +∞, limn→+∞ (Γn )−1 Γ⌊cΓn ⌋ = 0, and the conclusion follows from comP⌊cΓn ⌋−1
γi ≤
bining in (15), this limit, limk→+∞ γk = 0, limn→+∞ Γn = +∞ and
i=1
cγ1 Γn .
8

In the case of constant step sizes γk = γ for all k ≥ 1, we can deduce from Theorem 5,
a bound between π and the stationary distribution πγ of Rγ .
Corollary 8. Assume H 1 and H 2. Let (γk )k≥1 be a constant sequence γk = γ for all
k ≥ 1 with γ ≤ 1/(m + L). Then

W22 (π, πγ ) ≤ 2κ−1 L2 γ κ−1 + γ (2d + dL2 γ/m + dL2 γ 2 /6) .

Proof. Since by Proposition 3, for all x ∈ Rd , (δx Rγn )n≥0 converges to πγ as n → ∞
in (P2 (Rd ), W2 ), the proof then follows from Theorem 5 and Lemma 7 applied with
ℓ = 1.
We can improve the bound provided by Theorem 5 under additional regularity assumptions on the potential U .
H 3. The potential U is three times continuously differentiable and there exists L̃ such
that for all x, y ∈ Rd , ∇2 U (x) − ∇2 U (y) ≤ L̃ kx − yk.
Note that under H1 and H3, we have that for all x, y ∈ Rd ,
~
∇2 U (x)y ≤ L kyk , ∆(∇U
)(x)

2

≤ d2 L̃2 .

(16)

Theorem 9. Assume H 1, H 2 and H 3. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 1/(m + L). Then for all x ∈ Rd and n ≥ 1,
n
o
⋆ 2
W22 (δx Qnγ , π) ≤ u(1)
(γ)
kx
−
x
k
+
d/m
+ u(3)
n
n (γ) ,
(1)

where un is given by (10), κ in (4) and
(
"
!
)
n


X
dL̃2
4L4
2
−1
3
(3)
4 γi
4
−1
dγi 2L + κ
un (γ) =
+ γi L
+ γi L +
+m
3
3m
6
i=1
#
n 
Y
κγk 
1−
×
.
2

(17)

k=i+1

Proof. The proof of is along the same lines than Theorem 5, using Lemma 21 in place
of Lemma 20.
In the case of constant step sizes γk = γ for all k ≥ 1, we can deduce from Theorem 9,
a sharper bound between π and the stationary distribution πγ of Rγ .
Corollary 10. Assume H 1, H 2 and H 3. Let (γk )k≥1 be a constant sequence γk = γ
for all k ≥ 1 with γ ≤ 1/(m + L). Then
!
)
(
4
2
4L
d
L̃
+ γL4 (γ/6 + m−1 ) .
+ γL4 +
W22 (π, πγ ) ≤ 2κ−1 dγ 2 2L2 + κ−1
3
3m
9

Proof. The proof follows the same line as the proof of Corollary 8 and is omitted.
Using Proposition 3-(ii) and Corollary 6 or Corollary 10, given ε > 0, we determine the smallest number of iterations nε and an associated step-size γε satisfying
W2 (δx⋆ Rγnε , π) ≤ ε for all n ≥ nε . The dependencies on the dimension d and the
precision ε of nε based on Theorem 5 and Theorem 9 are reported in Table 1. Under
H 1 and H 2, the complexity we get is the same as the one established in [11] for the
total variation distance. However if in addition H 3 holds, we improve this complexity
with respect to the precision parameter ε as well. Finally, in the case where L̃ = 0 (for
example for non-degenerate d-dimensional Gaussian distributions), then the dependency
on the dimension of the number of iterations nε given by Theorem 9 turns out to be
only of order O(d1/2 log(d)).
Details and further discussions are included in the supplementary paper [12, Section 6-Section 7-Section 8]. In particular, the dependencies of the obtained bounds with
respect to the constants m and L which appear in H1, H2 are evidenced.
In a recent work [9] (based on a previous version of this paper), an improvement of
the proof of Theorem 5 has been proposed for constant step size. Whereas the constants
are sharper, the dependence with respect to the dimension and the precision parameter
ε > 0 is the same (first line of Table 1).
Parameter
Theorem 5 and Proposition 3-(ii)
Theorem 9 and Proposition 3-(ii)

d
O(d log(d))
O(d log(d))

ε
O(ε−2 |log(ε)|)
O(ε−1 |log(ε)|)

Table 1: Dependencies of the number of iterations nε to get W2 (δx⋆ Rγnεε , π) ≤ ε
For simplicity, consider sequences (γk )k≥1 defined for all k ≥ 1 by γk = γ1 k−α , for
(1)
(2)
γ1 < 1/(m + L) and α ∈ (0, 1]. Then for n ≥ 1, un ≤ O(e−κΓn /2 ), un ≤ dO(n−α ) and
(3)
un ≤ d2 O(n−2α ) (see [12, Sections 6 to 8] for details). Based on these upper bounds,
we can obtained the convergence rates provided by Theorem 5 and Theorem 9 for this
particular setting, see Table 2.
Theorem 5
Theorem 9

α ∈ (0, 1)
d O(n−α )
d2 O(n−2α )

α=1
d O(n−1 ) for γ1 > 2κ−1 see [12, Section 8]
d2 O(n−2 )

Table 2: Order of convergence of W22 (δx⋆ Qnγ , π) for γk = γ1 k−α

3

Quantitative bounds in total variation distance

We deal in this section with quantitative bounds in total variation distance. For Bayesian
inference application, total variation bounds are of importance for computing highest
posterior density (HPD) credible regions and intervals. For computing such bounds we
10

will use the results of Section 2 combined with the regularizing property of the semigroup
(Pt )t≥0 . Under H2, define χm for all t ≥ 0 by
p
χm (t) = (4/m)(e2mt − 1) .

The first key result consists in upper-bounding the total variation distance kµPt −νPt kTV
for µ, ν ∈ P1 (Rd ). To that purpose, we use the coupling by reflection; see [26] and [6,
Example 3.7]. This coupling is defined as the unique solution (Xt , Yt )t≥0 of the SDE:
(
√
dXt = −∇U (Xt )dt + 2dBtd
√
where et = e(Xt − Yt )
dYt = −∇U (Yt )dt + 2(Id −2et eTt )dBtd ,

with X0 = x, Y0 = y, e(z) = z/ kzk for z 6= 0 and e(0) = 0 otherwise. Define the
coupling time Tc = inf{s ≥ 0 R| Xs = Ys }. By construction Xt = Yt for t ≥ Tc . By
t
Levy’s characterization, B̃td = 0 (Id −2es eTs )dBsd is a d-dimensional Brownian motion,
therefore (Xt )t≥0 and (Yt )t≥0 are weak solutions to (1) started at x and y respectively.
Then by Lindvall’s inequality, for all t > 0 we have kPt (x, ·) − Pt (y, ·)kTV ≤ P (Xt 6= Yt ).
Theorem 11. Assume H1 and H2.
(i) For any x, y ∈ Rd and t > 0, it holds
kPt (x, ·) − Pt (y, ·)kTV ≤ 1 − 2Φ{− kx − yk /χm (t)} ,
where Φ is the cumulative distribution function of the standard normal distribution.
(ii) For any µ, ν ∈ P1 (Rd ) and t > 0,
kµPt − νPt kTV ≤ p

W1 (µ, ν)
(2π/m)(e2mt − 1)

.

(iii) For any x ∈ Rd and t ≥ 0,

(d/m)1/2 + kx − x⋆ k
kπ − δx Pt kTV ≤ p
.
(2π/m)(e2mt − 1)

(i) We now bound for t > 0, P (Xt 6= Yt ). For t < Tc , denoting by B1t =
T
d
0 1{s<Tc } es dBs , we get
√
d{Xt − Yt } = − {∇U (Xt ) − ∇U (Yt )} dt + 2 2et dB1t .

Proof.
Rt

By Itô’s formula and H2, we have for t < Tc ,
Z t
√
h∇U (Xs ) − ∇U (Ys ), es i ds + 2 2B1t
kXt − Yt k = kx − yk −
0
Z t
√
kXs − Ys k ds + 2 2B1t .
≤ kx − yk − m
0

11

By Grönwall’s inequality, we have
−mt

kXt − Yt k ≤ e

√ 1
√ Z t 1 −m(t−s)
kx − yk + 2 2Bt − m2 2
Bs e
ds .
0

Therefore by integration by part, kXt − Yt k ≤ Ut where (Ut )t∈(0,Tc ) is the one-dimensional
Ornstein-Uhlenbeck process defined by
Z 8t
√ Z t m(s−t) 1
−mt
−mt
Ut = e
kx − yk + 2 2
e
dBs = e
kx − yk +
em(s−t) dB̃s1
0

0

Therefore, for all x, y ∈ Rd and t ≥ 0, we get


P(Tc > t) ≤ P min Ut > 0 .
0≤s≤t

Finally the proof follows from [2, Formula 2.0.2, page 542].
(ii) Let µ, ν ∈ P1 (Rd ) and ξ ∈ Π(µ, ν) be an optimal transference plan for (µ, ν)
w.r.t. W1 . Since for all s > 0, 1/2 − Φ(−s) ≤ (2π)−1/2 s, (i) implies that for all x, y ∈ Rd
and t > 0,
Z
kx − yk
p
kµPt − νPt kTV ≤
dξ(x, y) ,
(2π/m)(e2mt − 1)
Rd ×Rd
which is the desired result.
(iii) The proof is a straightforward consequence of (ii) and Proposition 1-(iv).
Since for all s > 0, s ≤ es − 1, note that Theorem 11-(ii) implies that for all t > 0
and µ, ν ∈ P1 (Rd ),
kµPt − νPt kTV ≤ (4πt)−1/2 W1 (µ, ν) .
(18)

Therefore for all bounded measurable function f , Pt f is a Lipschitz function for all t > 0
with Lipshitz constant
kPt f kLip ≤ (4πt)−1/2 osc(f ) .
(19)
For all n, ℓ ≥ 1, n < ℓ and (γk )k≥1 a nonincreasing sequence denote by


ℓ
Y

−1
−1
Λn,ℓ (γ) = κ
(1 − κγj ) − 1 ,
Λℓ (γ) = Λ1,ℓ (γ) .



(20)

j=n

We will now study the contraction of Qn,ℓ
γ in total variation for non-increasing sequences
(γk )k≥1 . Strikingly, we are able to derive results which closely parallel Theorem 11. The
proof is nevertheless completely different because the reflection coupling is no longer
applicable in discrete time. We use a coupling construction inspired by the method of
[4, Section 3.3] for Gaussian random walks. This construction has been used in [13]
to establish convergence of homogeneous Markov chain in Wasserstein distances using
different method of proof. To avoid disrupting the stream of the proof, this construction
is introduced in Section 8.
12

Theorem 12. Assume H1 and H2.
(i) Let (γk )k≥1 be a nonincreasing sequence satisfying γ1 ≤ 2/(m + L). Then for all
x, y ∈ Rd and n, ℓ ∈ N∗ , n < ℓ, we have
n,ℓ
1/2
kδx Qn,ℓ
}.
γ − δy Qγ kTV ≤ 1 − 2Φ{− kx − yk /{8 Λn,ℓ (γ)}

(ii) Let (γk )k≥1 be a nonincreasing sequence satisfying γ1 ≤ 2/(m + L). For all µ, ν ∈
P1 (Rd ) and ℓ, n ∈ N∗ , n < ℓ,
n,ℓ
−1/2
kµQn,ℓ
W1 (µ, ν) .
γ − νQγ kTV ≤ {4πΛn,ℓ (γ)}

(iii) Let γ ∈ (0, 2/(m + L)]. Then for any x ∈ Rd and n ≥ 1,
o
n
kπγ − δx Rγn kTV ≤ {4πκ(1− (1− κγ)n/2 )}−1/2 (1− κγ)n/2 kx − x⋆ k + (2κ−1 d)1/2 .
Proof.

(i) By (3) for all x, y and k ≥ 1, we have
kx − γk ∇U (x) − y + γk ∇U (y)k ≤ (1 − κγk )1/2 kx − yk .

Let n, ℓ ≥ 1, n < ℓ, then applying Theorem 23 in Section 8, we get


n,ℓ
1/2
,
kδx Qn,ℓ
γ − δy Qγ kTV ≤ 1 − 2Φ − kx − yk /{8 Λn,ℓ (γ)}
(ii) Let f ∈ Fb (Rd ) and ℓ > n ≥ 1. For all x, y ∈ Rd by definition of the total
variation distance and (i), we have
n,ℓ
n,ℓ
n,ℓ
Qn,ℓ
γ f (x) − Qγ f (y) ≤ osc(f )kδx Qγ − δy Qγ kTV
o
n

,
≤ osc(f ) 1 − 2Φ − kx − yk /{8 Λn,ℓ (γ)}1/2

Using that for all s > 0, 1/2 − Φ(−s) ≤ (2π)−1/2 s concludes the proof.

(iii) The proof follows from (iii), the bound for all s > 0, 1/2 − Φ(−s) ≤ (2π)−1/2 s
and Proposition 2-(ii).
We can combine this result and Theorem 5 or Theorem 9 to get explicit bounds in
total variation between the Euler-Maruyama discretization and the target distribution
π. To that purpose, we use the following decomposition, for all nonincreasing sequence
(γk )k≥1 , initial point x ∈ Rd and ℓ ≥ 0:
kπ − δx Qℓγ kTV ≤ kπ − δx PΓℓ kTV + kδx PΓℓ − δx Qℓγ kTV .

13

(21)

The first term is dealt with Theorem 11-(iii). It remains to bound the second term
in (21). Since we will use Theorem 5 and Theorem 9, we have two different results
depending on the assumptions on U again. Define for all x ∈ Rd and n, p ∈ N,
2
ϑ(1)
n,p (x) = L

n
X

γi2

i=1

ϑ(2)
n,p (x)

=

n
X

γi3

i=1

n
Y

(1 − κγk /2)

k=i+1

n
Y



κ−1 + γi (2d + dL2 γi2 /6)

(22)



+L2 γi δi,n,p (x) κ−1 + γi


(1 − κγk /2) L4 δi,n,p (x)(4κ−1 /3 + γn+1 )

(23)

k=i+1

n
oi
2
+d 2L2 + 4κ−1 (dL̃2 /12 + γn+1 L4 /4) + γn+1
L4 /6
,

where ̺n,p (x) is given by (8) and

δi,n,p (x) = e−2mΓi−1 ̺n,p (x) + (1 − e−2mΓi−1 )(d/m) ,
Theorem 13. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 1/(m + L). Then for all x ∈ Rd and ℓ, n ∈ N∗ , ℓ > n,
kδx PΓℓ − δx Qℓγ kTV ≤ (ϑn (x)/(4πΓn+1,ℓ ))1/2
+ 2−3/2 L

ℓ
X


(γk3 L2 /3)̺1,k−1 (x) + dγk2

k=n+1

!1/2

. (24)

(2)

where ̺1,n (x) is defined by (8), ϑn (x) is equal to ϑn,0 (x) given by (23), if H3 holds, and
(1)

to ϑn,0 (x) given by (22) otherwise.
Proof. Applying Lemma 20 or Lemma 21, we get that for all x ∈ Rd
( (1)
ϑn,0 (x) (H1, H2) ,
1/2
n
W1 (δx Qγ , δx PΓn ) ≤ {ϑn (x)}
, ϑn (x) =
(2)
ϑn,0 (x) (H1, H2, H3) ,

(25)

By the triangle inequality, we get
δx PΓℓ − δx Qℓγ

TV

≤



δx PΓn − δx Qnγ PΓn+1,ℓ

Using (18) and (25), we have

δx PΓn − δx Q1,n
PΓn+1,ℓ
γ

TV

n
o
+ δx Qnγ Qγn+1,ℓ − PΓn+1,ℓ

TV

14

≤ (ϑn (x)/(4πΓn+1,ℓ ))1/2 .

TV

. (26)

(27)

For the second term, by [10, Equation 11] (note that we have a different convention for
the total variation distance) and the Pinsker inequality, we have
n
o 2
n+1,ℓ
Q
−
P
δx Q1,n
Γ
γ
γ
n+1,ℓ
TV

−3

≤2

L

2


Z
ℓ
X
3
(γk /3)

k=n+1

⋆

Rd

2

+

dγk2

(γk3 L2 /3)̺1,k−1 (x) + dγk2

.

k∇U (z) − ∇U (x )k

Qγk−1 (x, dz)



.

By H1 and Proposition 2, we get
n
o
n+1,ℓ
Q
−
P
δx Q1,n
Γ
γ
γ
n+1,ℓ

2
TV

≤ 2−3 L2

ℓ
X


k=n+1

Combining the last inequality and (27) in (26) concludes the proof.
Consider the case of decreasing step sizes defined: γk = γ1 k−α for k ≥ 1 and α ∈
(0, 1]. Under H 1 and H 2 , choosing n = ℓ − ⌊ℓα ⌋ in the bound given by Theorem 13
and using Table 2 implies that kδx⋆ Qℓγ − πkTV = d1/2 O(ℓ−α/2 ). Note that for α = 1,


this rate is true only for γ1 > 2κ−1 . If in addition H 3 holds, choosing n = ℓ − ℓα/2
in the bound given by (24) and using Table 2, (21) and Theorem 11-(iii) implies that
kδx⋆ Qℓγ − πkTV = dO(ℓ−(3/4)α ). Note that these rates are explicit compared to the ones
obtained in [11, Proposition 4]. These conclusions and the dependency on the dimension
are summarized in Table 3.
H1, H2
γk = γ1

k−α ,

α ∈ (0, 1]

d1/2 O(ℓ−α/2 )

H1, H2 and H3
dO(ℓ−3α/4 )

Table 3: Order of convergence of kδx⋆ Qℓγ − πkTV for γk = γ1 k−α based on Theorem 13
 −1 
When γk = γ ∈ (0,
1/(m
+
L))
for
all
k
≥
1,
under
H
1
and
H
2,
for
ℓ
>
γ

−1
choosing n = ℓ − γ
implies that (see the supplementary document [12, Section 2.2])

1/2
kδx Rγℓ − δx Pℓγ kTV ≤ (4π)−1/2 γD1 (γ, d) + γ 3 D2 (γ)D3 (γ, d, x)
+ D4 (γ, d, x) , (28)

where




D1 (γ, d) = 2L2 κ−1 κ−1 + γ 2d + L2 γ 2 /6 , D2 (γ) = L4 κ−1 + γ
o
n


−1
D3 (γ, d, x) = (ℓ − γ −1 )e−mγ(ℓ−⌈γ ⌉−1) kx − x⋆ k2 + 2d(κγm)−1

D4 (γ, d, x) = 2−3/2 L [dγ(1 + γ)
n
oi1/2
−1
+(L2 γ 3 /3) (1 + γ −1 )(1 − κγ)ℓ−⌈γ ⌉ kx − x⋆ k2 + 2(1 + γ)κ−1 d
.

(29)

Using this bound and Theorem 11-(iii), the minimal number of iterations ℓε > 0 for a
target precision ε > 0 to get kδx⋆ Rγℓεε − πkTV ≤ ε is of order d log(d)O(|log(ε)| ε−2 ) (the
proper choice of the step size γε is given in Table 5). In addition, letting ℓ go to infinity
in (28) we get the following result.
15

Corollary 14. Assume H1 and H2. Let γ ∈ (0, 1/(m + L)]. Then it holds

1/2
kπγ − πkTV ≤ 2−3/2 L dγ(1 + γ) + 2(L2 γ 3 /3)(1 + γ)κ−1 d

1/2
+ (4π)−1/2 γD1 (γ, d) + 2dγ 2 D2 (γ)(κm)−1
,

where D1 (γ) and D2 (γ) are given in (29).

Note that the bound provided by Corollary 14 is of order d1/2 O(γ 1/2 ).
If in addition H3 holds, for constant step sizes, we can improve with respect to the
step size γ, the bounds given by Corollary 14.
Theorem 15. Assume H1, H2 and H3. Let γ ∈ (0, 1/(m + L)]. Then it holds

1/2
kπγ − πkTV ≤ (4π)−1/2 γ 2 E1 (γ, d) + 2dγ 2 E2 (γ)/(κm)



1/2
+ (4π)−1/2 log γ −1 / log(2) γ 2 E1 (γ, d) + γ 2 E2 (γ)(2κ−1 d + d/m)

+ 2−3/2 L 2dγ 3 L2 /(3κ) + dγ 2

1/2

,

where E1 (γ, d) and E2 (γ) are defined by
n
o
E1 (γ, d) = 2dκ−1 2L2 + 4κ−1 (dL̃2 /12 + γL4 /4) + γ 2 L4 /6
E2 (γ) = L4 (4κ−1 /3 + γ) .

Proof. The proof is postponed to Section 7.
Note that the bound provided by Theorem 15 is of order dO(γ |log(γ)|), improving
the dependency given by Corollary 14, with respect to the step size γ. Furthermore when
L̃ = 0, this bound given by Theorem 15 is of order d1/2 O(γ |log(γ)|) and is sharp up
to a logarithmic factor. Indeed, assume that π is the d-dimensional standard Gaussian
distribution. In such case, the ULA sequence (Xk )k≥0 is the autoregressive process
√
given for all k ≥ 0 by Xk+1 = (1 − γ)Xk + 2γZk+1 . For γ ∈ (0, 1), this sequence
has a stationary distribution πγ , which is a d-dimensional Gaussian distribution with
zero-mean and covariance matrix σγ2 Id , with σγ2 = (1 − γ/2)−1 . Therefore, using [22,
Lemma 4.9] (or the Pinsker inequality), we get the following upper bound: kπ −πγ kTV ≤
Cd1/2 |σγ2 − 1| = Cd1/2 γ/2, where C is a universal constant.
We can also for a precision target ε > 0 choose γε > 0 and the number of iterations
nε > 0 to get kδx Rγnεε −πkTV ≤ ε. By Theorem 11-(iii), Theorem 12-(iii) and Theorem 15,
the minimal number of iterations ℓε is of order d log2 (d)O(ε−1 log2 (ε)) for a well chosen
step size γε . This result improves the conclusion of [11] and 14 with respect to the
precision parameter ε, which provides an upper bound of the number of iterations of order
d log(d)O(ε−2 log2 (ε)). We can also compare our reported upper bound with the one
obtained for the d-dimensional standard Gaussian distribution. If the initial distribution
is the Dirac mass at zero (the minimum of the potential U (x) = kxk2 /2) and γ ∈ (0, 1),
the distribution of the ULA sequence after n iterations is zero-mean Gaussian with
16

covariance (1 − (1 − γ)2(n+1) )/(1 − γ/2) Id . If we use [22, Lemma 4.9] again, we get for
γ ∈ (0, 1),
kδ0 Rγn − πkTV ≤ Cd1/2 γ|1 − 2γ −1 (1 − γ)2(n+1) | ,
where C is a universal constant. To get an ε precision we need to choose γε = d−1/2 ε/(2C)
and then nε = ⌈(1/2) log(γε /4)/ log(1 − γε )⌉ = d1/2 log(d)O(ε−1 | log(ε)|). On the other
hand since L̃ = 0, based on the bound given by Theorem 15, the minimal number of
iterations to get kδx Rγnεε − πkTV ≤ ε is of order d1/2 log2 (d)O(ε−1 log2 (ε)). It follows that
our upper bound for the stepsize and the optimal number of iterations is again sharp
up to a logarithmic factor in the dimension and the precision. The discussions on the
bounds for constant sequences of step sizes are summarized in Table 4 and Table 5.
H1, H2
kπ − πγ kTV

d1/2 O(γ 1/2 )

H1, H2 and H3
dO(γ |log(γ)|)

Table 4: Order of the bound between π and πγ in total variation function of the step
size γ > 0 and the dimension d.

γε
nε

H1, H2 and H3
d−1 log−1 (d)O(ε log−1 (ε) )
d log2 (d)O(ε−1 log2 (ε))

H1, H2
d−1 O(ε2 )
d log(d)O(ε−2 |log(ε)|)

Table 5: Order of the step size γε > 0 and the number of iterations nε ∈ N∗ to get
kδx⋆ Rγnεε − πkTV ≤ ε for ε > 0.

4

Mean square error and concentration for bounded measurable functions

R
The result of the previous section allows us to study the approximation of Rd f (y)π(dy)
by the weighted average estimator π̂nN (f ) defined, for a measurable and bounded function
f : Rd → R, N, n ∈ N, n ≥ 1 by
π̂nN (f ) =

N
+n
X

N
f (Xk ) ,
ωk,n

N
= γk+1 Γ−1
ωk,n
N +2,N +n+1 .

(30)

k=N +1

We derive a bound on the mean-square error, defined as MSEN,n
= Ex
f
The

MSEN,n
f

h

π̂nN (f ) − π(f )

can be decomposed as the sum of the squared bias and variance:

MSEN,n
= Ex [π̂nN (f )] − π(f )
f

17

2


+ Varx π̂nN (f )

.

2

i

.

We first obtain a bound for the bias. By the Jensen inequality and because f is bounded,
we have:
N
+n
X
2
N
Ex [π̂nN (f )] − π(f ) ≤ osc(f )2
kδx Qkγ − πk2TV ;
ωk,n
k=N +1

Using the results of Section 3, we can deduce different bounds for the bias, depending
on the assumptions on U and the sequence of step sizes (γk )k≥1 . It is now required to
derive a bound for the variance. Our main tool is the Gaussian Poincaré inequality [3,
Theorem 3.20] which can be applied to Rγ (y, ·) defined by (6), noticing that Rγ (y, ·)
is a Gaussian distribution with mean y − γ∇U (y) and covariance matrix 2γ Id : for all
Lipschitz function g : Rd → R
Rγ {g(·) − Rγ g(y)}2 (y) ≤ 2γ kgk2Lip .

(31)

Lemma 16. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with γ1 ≤
2/(m + L). Let g : Rd → R be a Lipschitz function. Then for all n, p ≥ 1, n ≤ p and
y ∈ Rd
Z

2
n,p
Qn,p
≤ 2κ−1 kgk2Lip ,
0≤
γ (y, dz) g(z) − Qγ g(y)
Rd

where

Qn,p
γ

is given by (7).

Proof. By decomposing g(Xp ) − EGy n [g(Xp )] =
and using EGy k [g(Xp )] = Qγk+1,p g(Xk ), we get
VarGy n

{g(Xp )} =
=

p
X

k=n+1
p
X

EGy n
EGy n

k=n+1





G
Ey k−1

Rγk

n

Pp

Gk
k=n+1 {Ey

G

[g(Xp )] − Ey k−1 [g(Xp )]},


2 
Gk−1
Gk
Ey [g(Xp )] − Ey
[g(Xp )]
Qγk+1,p g(·)

−

o2

Rγk Qγk+1,p g(Xk−1 )


(Xk−1 ) .

P
Equation (31) implies VarGy n {g(Xp )} ≤ 2 pk=n+1 γk kQγk+1,p gk2Lip . The proof follows
from Corollary 4 and Lemma 7, using the bound (1 − t)1/2 ≤ 1 − t/2 for t ∈ [0, 1].
Theorem 17. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 2/(m + L). Then for all N ≥ 0, n ≥ 1, x ∈ Rd and f ∈ Fb (Rd ), we get
(4)

Varx {π̂nN (f )} ≤ osc(f )2 {2γ1 Γ−1
N +2,N +n+1 + uN,n (γ)}
(4)

uN,n (γ) =

NX
+n−1
k=N

γk+1

( N +n
X

N
ωi,n

(πΛk+2,i (γ))1/2
i=k+2

for n1 , n2 ∈ N, Λn1 ,n2 (γ) is given by (20).

18

)2

+ κ−1

(

N
+n
X

i=N +1

N
ωi,n

(4πΛN +1,i (γ))1/2

)2

,

Proof. Let N ≥ 0, n ≥ 1, x ∈ Rd and f ∈ Fb (Rd ). We decompose π̂nN (f ) − Ex [π̂nN (f )] as
the sum of martingale increments, w.r.t. (Gn )n≥0 , the natural filtration associated with
Euler approximation (Xn )n≥0 , and we get


NX
+n−1
 N

 N 2
Gk+1  N
Gk
Varx π̂n (f ) =
Ex Ex
π̂n (f ) − Ex π̂n (f )
k=N

+ Ex

h

EGxN





π̂nN (f )

−

2
Ex [π̂nN (f )]

i

. (32)





G
The martingale increment Exk+1 π̂nN (f ) − EGxk π̂nN (f ) has a simple expression. For
k = N + n − 1, . . . , N + 1, define (backward in time) the function
N
N
ΦN
n,k : xk 7→ ωk,n f (xk ) + Rγk+1 Φn,k+1 (xk ) ,

(33)

N
N
where ΦN
n,N +n : xN +n 7→ Φn,N +n (xN +n ) = ωN +n,n f (xN +n ). Denote finally
N
ΨN
n : xN 7→ RγN+1 Φn,N +1 (xN ) .

(34)

Note that for k ∈ {N, . . . , N + n − 1}, by the Markov property,



Gk+1  N
N
π̂n (f ) − EGxk π̂nN (f ) ,
ΦN
(35)
n,k+1 (Xk+1 ) − Rγk+1 Φn,k+1 (Xk ) = Ex


GN π̂ N (f ) . With these notations, (32) may be equivalently expressed
and ΨN
n (XN ) = Ex
n
as
NX
+n−1
i
h

 N
2
N
(X
)
(·)
−
R
Φ
(X
)
Ex Rγk+1 ΦN
Varx π̂n (f ) =
γ
k
k
n,k+1
k+1 n,k+1
k=N


+ Varx ΨN
n (XN )

.

(36)

Let k ∈ {N, . . . , N + n − 1}. We cannot directly apply the Poincaré inequality (31) since
the function ΦN
n,k , defined in (33), is not Lipschitz. However, Theorem 12-(ii) shows that
for all ℓ, n ∈ N∗ , n < ℓ, Qn,ℓ
γ f is a Lipschitz function with
Qn,ℓ
γ f

Lip

≤ osc(f )/{4πΛn,ℓ (γ)}1/2 .

N
N
N
Using (33), we may decompose ΦN
n,k = ωk+1,n f + Φ̃n,k , where Φ̃n,k =
which is Lipshitz with constant

Φ̃N
n,k

Lip

≤

N
+n
X

N
ωi,n
Qγk+2,i f

i=k+2

Using the inequality (a

+ b)2

≤

2a2

Lip

+


N
Rγk+1 ΦN
n,k+1 (·) − Rγk+1 Φn,k+1 (y)

≤ osc(f )

2b2 ,

N
+n
X

i=k+2

(37)
PN +n

k+2,i
N
f
i=k+2 ωi,n Qγ

N
ωi,n
/ {4πΛk+2,i (γ)}1/2 .

(38)

(31), we finally get for any y ∈ Rd

2

N
)2 osc(f )2
(y) ≤ 2(ωk+1,n
( N +n
)2
X
2
N
1/2
. (39)
+ γk+1 osc(f )
ωi,n /{πΛk+2,i (γ)}
i=k+2

19

3.5

3.0

3.0

Γn u0(4)
,n (γ)

Γn u0(4)
,n (γ)

3.5

(Γn u0(4)
,n (γ))n for γk =1/2 for all k ≥1

2.5

2.0

1.5

2.0

−α
(Γn u0(4)
,n (γ))n for γk =(1 + k)

1.5

1.00

1000

2000

3000

Iteration

n

4000

1.00

5000

3.5

3.5

3.0

3.0

2.5

2.0

−α
(Γn u0(4)
,n (γ))n for γk =(1 + k)

1.5

1.00

1000

2000

3000

Iteration

n

1000

2000

3000

Iteration

Γn u0(4)
,n (γ)

Γn u0(4)
,n (γ)

2.5

and

α =1/2

4000

α =1/4

4000

n

5000

2.5

2.0

−α
(Γn u0(4)
,n (γ))n for γk =(1 + k)

1.5

1.00

5000

and

1000

2000

3000

Iteration

n

and

α =3/4

4000

5000

(4)

Figure 1: Plots of (u0,n )n≥1 Γn for four sequences of step-sizes (γk )k≥1 , γk = (1 + k)−α
for α = 0, 1/4, 1/2, 3/4

N
N
It remains to control Varx ΨN
n (XN ) , where Ψn is defined in (34). Using (37), Ψn is
a Lipschitz function with Lipschitz constant bounded by:
ΨN
n Lip

≤

N
+n
X

i=N +1

N
ωi,n

QγN +1,i f Lip

≤ osc(f )

N
+n
X

N
ωi,n
/{4πΛN +1,i (γ)}1/2 .

(40)

i=N +1

By Lemma 16, we have the following result which is the counterpart of Corollary 3:
for all y ∈ Rd ,


−1
Vary ΨN
kf k2∞
n (XN ) ≤ 2κ

(

N
+n
X

N
ωi,n
/(πΛN +1,i )1/2

i=N +1

)2

.

(41)

Finally, the proof follows from combining (39) and (41) in (36).
(4)

To illustrate the result Theorem 17, we study numerically the behaviour (uN,n )n≥1
for κ = 1 N = 0, and four different non-increasing sequences of step-sizes (γk )k≥1 ,
γk = (1 + k)−α for α = 1/4, 1/2, 3/4 and γk = 1/2 for k ≥ 1. These results are
(4)
gathered in Figure 1, where it can be observed that (Γn u0,n (γ))n≥1 converges to a limit
as n → +∞.
We now establish an exponential deviation inequality for π̂nN (f ) − Ex [π̂nN (f )] given
by (30) for a bounded measurable function f .

20

Theorem 18. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 2/(m + L). Let (Xn )n≥0 be given by (2) and started at x ∈ Rd . Then for all N ≥ 0,
n ≥ 1, r > 0, and functions f ∈ Fb (Rd ):

where



2 (5)
−1 2
Px π̂nN (f ) ≥ Ex [π̂nN (f )] + r ≤ e−{r−osc(f )(ΓN+2,N+n+1 ) } /{2osc(f ) uN,n (γ)} ,
(5)

uN,n (γ) =

NX
+n−1
k=N

γk+1

( N +n
X

N
ωi,n

(πΛk+2,i )1/2
i=k+2

)2

+ κ−1

(

N
+n
X

N
ωi,n

i=N +1

(πΛN +1,i )1/2

)2

.

Proof. The proof is postponed in the supplementary document to Section 3.

5

Numerical experiments

Consider a binary regression set-up in which the binary observations (responses) {Yi }pi=1
β T Xi )}pi=1 ,
are conditionally independent Bernoulli random variables with success probability {̺(β
p
where ̺ is the logistic function defined for z ∈ R by ̺(z) = ez /(1 + ez ) and {Xi }i=1 and
β are d dimensional vectors of known covariates and unknown regression coefficients,
respectively. The prior distribution for the parameter β is a zero-mean Gaussian distribution with covariance matrix Σβ . The density of the posterior distribution of β is up
to a proportionality constant given by
!
p n
o
X
p
−1 T −1
T
β T Xi
β |{(Xi , Yi )}i=1 ) ∝ exp
) − 2 β Σβ β .
Yiβ Xi − log(1 + e
πβ (β
i=1

Bayesian inference for the logistic regression model has long been recognized as a numerically involved problem, due to the analytically inconvenient form of the likelihood function. Several algorithms have been proposed, trying to mimick the data-augmentation
(DA) approach of [1] for probit regression; see [20], [15] and [16]. Recently, a very
promising DA algorithm has been proposed in [31], using the Polya-Gamma distribution
in the DA part. This algorithm has been shown to be uniformly ergodic for the total
variation by [7, Proposition 1], which provides an explicit expression for the ergodicity
constant. This constant is exponentially small in the dimension of the parameter space
and the number of samples. Moreover, the complexity of the augmentation step is cubic
in the dimension, which prevents from using this algorithm when the dimension of the
regressor is large.
We apply ULA to sample from the posterior distribution πβ (·|{(Xi , Yi )}pi=1 ). The
gradient of its log-density may be expressed as
β |{Xi , Yi }pi=1 )}
∇ log{πβ (β

p 
X
=
Yi Xi −
i=1

21

Xi
T
1 + e−ββ Xi



− Σβ−1β ,

Therefore − log πβ (·|{XiP
, Yi }pi=1 ) is strongly convex H 2 with m = λ−1
max (Σβ ) and satisfies H 1 with L = (1/4) pi=1 XiT Xi + λ−1
(Σ
),
where
λ
(Σ
)
and
λmax (Σβ ) denote
min
β
β
min
the minimal and maximal eigenvalues of Σβ , respectively. We first compare the histograms produced by ULA and the Pòlya-Gamma Gibbs sampling from [31]. For that
purpose, we take d P
= 5, p = 100, generate synthetic data (Yi )1≤i≤p and (Xi )1≤i≤p , and
set Σβ−1 = (dp)−1 ( pi=1 XiT Xi ) Id . We produce 108 samples from the Pólya-Gamma
sampler using the R package BayesLogit [38]. Next, we make 103 runs of the Euler approximation scheme with n = 106 effective iterations, with a constant sequence (γk )k≥1 ,
γk = 10(κn1/2 )−1 for all k ≥ 0 and a burn-in period N = n1/2 . The histogram of
the Pólya-Gamma Gibbs sampler for first component, the corresponding mean of the
obtained histograms for ULA and the 0.95 quantiles are displayed in Figure 2. The
same procedure is also applied with the decreasing step size sequence (γk )k≥1 defined by
γk = γ1 k−1/2 , with γ1 = 10(κ log(n)1/2 )−1 and for the burn in period N = log(n), see also
Figure 2. In addition, we also compare MALA and ULA on five real data sets, which are
1.4

quantile at 95%
Mean of ULA

1.2

Polya−Gamma

1

quantile at 95%
Mean of ULA
Polya−Gamma

1.4
1.2
1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2
0.4

0.6

0.8

1

1.2

1.4

1.6

0.4

1.8

0.6

0.8

1

1.2

1.4

1.6

1.8

Figure 2: Empirical distribution comparison between the Polya-Gamma Gibbs Sampler
and ULA. Left panel: constant step size γk = γ1 for all k ≥ 1; right panel: decreasing
step size γk = γ1 k−1/2 for all k ≥ 1
summarized in Table 6. Note that for the Australian credit data set, the ordinal covariates have been stratified by dummy variables. Furthermore, we normalized the
Ppdata sets
−1
T
where
Σ
=
p
and consider the Zellner prior setting Σ−1 = (π2 d/3)Σ−1
X
i=1 Xi Xi
X
; see [34], [19] and the references therein. Also, we apply a pre-conditioned version of
1/2
MALA and ULA, targeting the probability density π̃β (·) ∝ πβ (ΣX ·). Then, we ob1/2
tain samples from πβ by post-multiplying the obtained draws by ΣX . We compare
MALA and ULA for each data sets by estimating for each component i ∈ {1, . . . , d} the
marginal accuracy between their d marginal empirical distributions and the d marginal
posterior distributions, where the marginal accuracy between two probability measure
µ, ν on (R, B(R)) is defined by
MA(µ, ν) = 1 − (1/2)kµ − νkTV .
This quantity has already been considered in [5] and [8] to compare approximate samplers. To estimate the d marginal posterior distributions, we run 2 · 107 iterations of
the Polya-Gamma Gibbs sampler. Then 100 runs of MALA and ULA (106 iterations
22

❵❵❵
❵❵❵
❵❵Dimensions
Observations p
❵❵❵
Data set
❵❵❵

German credit 1
Heart disease 2
Australian credit3
Pima indian diabetes4
Musk5

Covariates d

1000
270
690
768
476

25
14
35
9
167

Table 6: Dimension of the data sets
0.98

0.985

0.97
0.98
0.975

0.96
0.95
ULA

MALA

ULA

MALA

ULA

MALA

0.98
0.98

0.97

0.97

0.96

0.96
0.95

ULA

MALA

0.95

0.98
0.96
0.94
0.92
ULA

MALA

Figure 3: Marginal accuracy across all the dimensions.
Upper left: German credit data set. Upper right: Australian credit data set. Lower left:
Heart disease data set. Lower right: Pima Indian diabetes data set. At the bottom:
Musk data set
per run) have been performed. For MALA, the step-size is chosen so that the acceptance probability at stationarity is approximately equal to 0.5 for all the data sets. For
ULA, we choose the same constant step-size than MALA. We display the boxplots of the
mean of the estimated marginal accuracy across all the dimensions in Figure 3. These
results all imply that ULA is an alternative to the Polya-Gibbs sampler and the MALA
algorithm.
1

http://archive.ics.uci.edu/ml/datasets/Statlog+(German+Credit+Data)
http://archive.ics.uci.edu/ml/datasets/Statlog+(Heart)
3
http://archive.ics.uci.edu/ml/datasets/Statlog+(Australian+Credit+Approval)
4
http://archive.ics.uci.edu/ml/datasets/Pima+Indians+Diabetes
5
https://archive.ics.uci.edu/ml/datasets/Musk+(Version+1)
2

23

6
6.1

Proofs
Proof of Theorem 5

We preface the proof by two technical Lemmata.
Lemma 19. Let (Yt )t≥0 be the solution of (1) started at x ∈ Rd . For all t ≥ 0 and
x ∈ Rd ,
i
h
Ex kYt − xk2 ≤ dt(2 + L2 t2 /3) + (3/2)t2 L2 kx − x⋆ k2 .

Proof. Let A be the generator associated with (Pt )t≥0 defined by (S2).
R t Denote for all
x, y ∈ Rd , Ṽx (y) = ky − xk2 . Note that the process (Ṽx (Yt ) − Ṽx (x) − 0 A Ṽx (Ys )ds)t≥0 ,
is a (Ft )t≥0 -martingale under Px . Denote for all t ≥ 0 and x ∈ Rd by ṽ(t, x) = Pt Ṽx (x).
Then we get,
∂ṽ(t, x)
= Pt A Ṽx (x) .
(42)
∂t
By H2, we have for all y ∈ Rd , h∇U (y) − ∇U (x), y − xi ≥ m kx − yk2 , which implies


A Ṽx (y) = 2 (− h∇U (y), y − xi + d) ≤ 2 −mṼx (y) + d − h∇U (x), y − xi .
Using (42), this inequality and that Ṽx is positive, we get


Z
∂ṽ(t, x)
h∇U (x), y − xi Pt (x, dy) .
= Pt A Ṽx (x) ≤ 2 d −
∂t
Rd

(43)

By the Cauchy-Schwarz inequality, ∇U (x⋆ ) = 0, (1) and the Jensen inequality, we have,
|Ex [h∇U (x), Yt − xi]| ≤ k∇U (x)k kEx [Yt − x]k

Z t
⋆
≤ k∇U (x)k Ex
{∇U (Ys ) − ∇U (x )} ds
0

√
≤ t k∇U (x) − ∇U (x⋆ )k

Z

h

t

⋆

Ex k∇U (Ys ) − ∇U (x )k

0

2

i

ds

1/2

.

Furthermore, by H1 and Proposition 1-(i), we have
Z

Rd

h∇U (x), y − xi Pt (x, dy) ≤
≤

√

2

⋆

tL kx − x k



√

2

⋆

tL kx − x k

Z

0

t

h

⋆ 2

Ex kYs − x k

i

ds

1/2

1/2
1 − e−2mt
2tm + e−2mt − 1
⋆ 2
kx − x k +
(d/m)
2m
2m

≤ L2 kx − x⋆ k (t kx − x⋆ k + t3/2 d1/2 ) ,

where we used for the last line that by the Taylor theorem with remainder term, for all
s ≥ 0, (1 − e−2ms )/(2m) ≤ s and (2ms + e−2ms − 1)/(2m) ≤ ms2 , and the inequality
24

√
√
a + b ≤ a + b. Plugging this upper bound in (43), and since 2 kx − x⋆ k t3/2 d1/2 ≤
t kx − x⋆ k2 + t2 d, we get
√

∂ṽ(t, x)
≤ 2d + 3L2 t kx − x⋆ k2 + L2 t2 d
∂t
Since ṽ(0, x) = 0, the proof is completed by integrating this result.
Let (Ft )t≥0 be the filtration associated with (Bt )t≥0 and (Y0 , Y 0 ).
Lemma 20. Assume H 1 and H 2. Let (γk )k≥1 be a nonincreasing sequence with γ1 ≤
1/(m + L). Let ζ0 ∈ P2 (Rd × Rd ), (Yt , Y t )t≥0 such that (Y0 , Y 0 ) is distributed according
to ζ0 and given by (9). Then almost surely for all n ≥ 0 and ǫ > 0,
YΓn+1 − Y Γn+1

2

2

≤ {1 − γn+1 (κ − 2ǫ)} YΓn − Y Γn
Z Γn+1
−1
k∇U (Ys ) − ∇U (YΓn )k2 ds ,
+ (2γn+1 + (2ǫ) )
Γn
i
h
2
2
FΓn
≤ {1 − γn+1 (κ − 2ǫ)} YΓn − Y Γn
E
YΓn+1 − Y Γn+1


2
2
+ L2 γn+1
(1/(4ǫ) + γn+1 ) 2d + L2 γn+1 kYΓn − x⋆ k2 + dL2 γn+1
/6 .

(44)

(45)

Proof. Let n ≥ 0 and ǫ > 0, and set Θn = YΓn − Y Γn . We first show (44). By definition
we have:
2

2

kΘn+1 k = kΘn k +

Z

Γn+1
Γn



2

∇U (Ys ) − ∇U (Y Γn ) ds

− 2γn+1 Θn , ∇U (YΓn ) − ∇U (Y Γn ) − 2

Z

Γn+1

hΘn , {∇U (Ys ) − ∇U (YΓn )}i ds . (46)

Γn

Young’s inequality and Jensen’s inequality imply
Z

Γn+1

Γn



2

∇U (Ys ) − ∇U (Y Γn ) ds

2
≤ 2γn+1
∇U (YΓn ) − ∇U (Y Γn )

+ 2γn+1

Z

Γn+1

Γn

2

k∇U (Ys ) − ∇U (YΓn )k2 ds .

Using (3), γ1 ≤ 1/(m + L) and (γk )k≥1 is nonincreasing, (46) becomes
2

2

kΘn+1 k ≤ {1 − γn+1 κ} kΘn k + 2γn+1
−2

Z

Z

Γn+1
Γn

k∇U (Ys ) − ∇U (YΓn )k2 ds

Γn+1

hΘn , {∇U (Ys ) − ∇U (YΓn )}i ds .

Γn

Using the inequality | ha, bi | ≤ ǫkak2 + (4ǫ)−1 kbk2 concludes the proof of (44).
25

(47)

We now prove (45). Note that (44) implies that
i
h
EFΓn kΘn+1 k2 ≤ {1 − γn+1 (κ − 2ǫ)} kΘn k2
Z Γn+1
i
h
−1
EFΓn k∇U (Ys ) − ∇U (YΓn )k2 ds . (48)
+ (2γn+1 + (2ǫ) )
Γn

By H1, the Markov property of (Yt )t≥0 and Lemma 19, we have
Z Γn+1
i
h
EFΓn k∇U (Ys ) − ∇U (YΓn )k2 ds
Γn


2
4
3
+ dL2 γn+1
/12 + (1/2)L2 γn+1
kYΓn − x⋆ k2 .
≤ L2 dγn+1
The proof is then concluded plugging this bound in (48) .

6.2

Proofs of Theorem 9

Lemma 21. Assume H 1, H 2 and H 3. Let (γk )k≥1 be a nonincreasing sequence with
γ1 ≤ 1/(m + L). and ζ0 ∈ P2 (Rd × Rd ). Let (Yt , Y t )t≥0 be defined by (9) such that
(Y0 , Y 0 ) is distributed according to ζ0 . Then for all n ≥ 0 and ǫ > 0, almost surely
h
i
2
2
EFΓn YΓn+1 − Y Γn+1
≤ {1 − γn+1 (κ − 2ǫ)} YΓn − Y Γn
n
2
3
d(2L2 + ǫ−1 (dL̃2 /12 + γn+1 L4 /4) + γn+1
L4 /6)
+ γn+1
o
+L4 (ǫ−1 /3 + γn+1 ) kYΓn − x⋆ k2 .

Proof. Let n ≥ 0 and ǫ > 0, and set Θn = YΓn − Y Γn . Using Itô’s formula, we have for
all s ∈ [Γn , Γn+1 ),
Z sn
o
~
∇2 U (Yu )∇U (Yu ) + ∆(∇U
)(Yu ) du
∇U (Ys ) − ∇U (YΓn ) =
Γn
√ Z s 2
∇ U (Yu )dBu .
(49)
+ 2
Γn

Rs
Since Θn is FΓn -measurable and ( 0 ∇2 U (Yu )dBu )s∈[0,Γn+1 ] is a (Fs )s∈[0,Γn+1 ] -martingale
under H1, by (49) we have:
EFΓn [hΘn , ∇U (Ys ) − ∇U (YΓn )i]

Z
FΓn
= Θn , E

s
Γn

n

o 
~
∇ U (Yu )∇U (Yu ) + ∆(∇U )(Yu ) du
2

Combining this equality and | ha, bi | ≤ ǫkak2 + (4ǫ)−1 kbk2 in (47) we have
i
h
EFΓn kΘn+1 k2 ≤ {1 − γn+1 (κ − 2ǫ)} kΘn k2 + +(2ǫ)−1 A

Z Γn+1
2γn+1 EFΓn
k∇U (Ys ) − ∇U (YΓn )k2 ds ,
Γn

26

(50)

where
A=

Z

Γn+1

E

Z

FΓn

Γn

s


~
∇ U (Yu )∇U (Yu ) + (1/2)∆(∇U )(Yu )du
2

Γn

2

ds .

We now separately bound the two last terms of the right hand side. By H1, the Markov
property of (Yt )t≥0 and Lemma 19, we have
Z Γn+1
i
h
EFΓn k∇U (Ys ) − ∇U (YΓn )k2 ds
Γn


2
4
3
+ dL2 γn+1
/12 + (1/2)L2 γn+1
kYΓn − x⋆ k2 . (51)
≤ L2 dγn+1

We now bound A. We get using Jensen’s inequality, Fubini’s theorem, ∇U (x⋆ ) = 0 and
(16)
Z s
Z Γn+1
h
i
2
du ds
EFΓn ∇2 U (Yu )∇U (Yu )
(s − Γn )
A≤2
Γn

−1

+2
≤2

Z

Γn+1
Γn

Z

Γn
Γn+1

Γn

(s − Γn )L

4

Z

(s − Γn )
s

Γn

Z

s

E

FΓn

Γn



~
∆(∇U
)(Yu )

2



du ds

h
i
3
EFΓn kYu − x⋆ k2 du ds + γn+1
d2 L̃2 /6 .

(52)

By Lemma 19-(i), the Markov property and for all t ≥ 0, 1 − e−t ≤ t, we have for all
s ∈ [Γn , Γn+1 ],
Z s
i
h
EFΓn kYu − x⋆ k2 du ≤ (2m)−1 (1 − e−2m(s−Γn ) ) kYΓn − x⋆ k2 + d(s − Γn )2 .
Γn

Using this inequality in (52) and for all t ≥ 0, 1 − e−t ≤ t , we get

3
4
3
A ≤ (2L4 γn+1
/3) kYΓn − x⋆ k2 + L4 dγn+1
/2 + γn+1
d2 L̃2 /6 .

Combining this bound and (51) in (50) concludes the proof.

7

Proof of Theorem 15

We preface the proof by a preliminary lemma. Define for all γ > 0, the function n :
R∗+ → N by



n(γ) = log γ −1 / log(2) .
(53)

Lemma 22. Assume H 1, H 2 and H 3. Let γ ∈ (0, 1/(m + L)). Then for all x ∈ Rd
and ℓ ∈ N∗ , ℓ > 2n(γ) ,
(2)

kδx Pℓγ − δx Rγℓ kTV ≤ (ϑℓ−2n(γ) ,0 (x)/(π2n(γ)+2 γ))1/2
−3/2

+2



3

2

L (γ L /3)̺1,ℓ−1 (x) + dγ

2 1/2

+

n(γ)
X
k=1

27

(2)

(ϑ2k−1 ,ℓ−2k (x)/(π2k+1 γ))1/2 .

(2)

where ̺1,ℓ−1 (x) is defined by (8) and for all n1 , n2 ∈ N, ϑn1 ,n2 is given by (23).

Proof. Let γ ∈ (0, 1/(m + L)) and ℓ ∈ N∗ . For ease of notation, let n = n(γ), and
assume that ℓ > 2n . Consider the following decomposition
o
n
n
δx Pℓγ − δx Rγℓ
δx P(ℓ−2n )γ − δx Rγℓ−2 P2n γ
≤
TV

TV

+ δx Rγℓ−1 {Pγ − Rγ }

TV

+

n
X

ℓ−2k

δx Rγ

k=1

n

k−1

P2k−1 γ − Rγ2

o

P2k−1 γ

TV

. (54)

We bound each term in the right hand side. First by (18) and Equation (25), we have
o
n
n
(2)
≤ (ϑℓ−2n ,0 (x)/(π2n+1 γ))1/2 ,
(55)
δx P(ℓ−2n )γ − δx Rγℓ−2 P2n γ
TV

(2)

where ϑn,0 (x) is given by (23). Similarly but using in addition Proposition 2, we have
for all k ∈ {1, · · · , n},
o
n
k−1
k
(2)
P2k−1 γ
≤ (ϑ2k−1 ,ℓ−2k (x)/(π2k+1 γ))1/2 ,
(56)
δx Rγℓ−2 P2k−1 γ − Rγ2
TV

(2)

where ϑ2k−1 ,ℓ−2k (x) is given by (23). For the last term, by [10, Equation 11] and the
Pinsker inequality, we have


Z
2
≤ 2−3 L2 (γ 3 /3)
k∇U (z)k2 Rγℓ−1 (x, dz) + dγ 2 .
δx Rγℓ−1 {Pγ − Rγ }
TV

Rd

By H1 and Proposition 2, we get
δx Rγℓ−1 {Rγ − Pγ }

2
TV


≤ 2−3 L2 (γ 3 L2 /3)̺1,ℓ−1 (x) + dγ 2

.

(57)

Combining (55), (56) and (57) in (54) concludes the proof.

Proof of Theorem 15. By Lemma 22, we get (see the supplementary document [12, Section 7] for the detailed calculation):
o1/2
n
n
o
kδx Pℓγ − δx Rγℓ kTV ≤ 2−3/2 L (γ 3 L2 /3) (1 − κγ)ℓ−1 kx − x⋆ k2 + 2κ−1 d + dγ 2

1/2
+ (4π)−1/2 γ 2 E1 (γ, d) + γ 3 E2 (γ)E3 (γ, d, x)

1/2
+ n(γ)(4π)−1/2 γ 2 E1 (γ, d) + γ 2 E2 (γ)E4 (γ, d, x)
, (58)

where

E3 (γ, d, x) = (ℓ − γ −1 )e−mγ(ℓ−2γ

E4 (γ, d, x) = e−mγ(ℓ−2γ

−1 −1)

−1 −1)

kx − x⋆ k2 + 2d/(κγm)

kx − x⋆ k2 + 2κ−1 d + d/m .

Letting ℓ go to infinity, using Theorem 11-(iii) and Theorem 12-(iii), we get the desired
conclusion.

28

8

Contraction in total variation for functional autoregressive models

In this section, we consider functional autoregressive models of the form: for all k ≥ 0
Xk+1 = hk+1 (Xk ) + σk+1 Zk+1 ,

(59)

where (Zk )k≥1 is a sequence of i.i.d. d dimensional standard Gaussian random variables,
(σk )k≥1 is a sequence of positive real numbers and (hk )k≥1 are a sequence of measurable
map from Rd to Rd which satisfies the following assumption:
AR 1. For all k ≥ 1, there exists ̟k ∈ R such that hk is 1 − ̟k -Lipschitz, i.e. for all
x, y ∈ Rd , khk (x) − hk (y)k ≤ (1 − ̟k ) kx − yk.
The sequence {Xk , k ∈ N} defines an inhomogeneous Markov chain associated with
the sequence of Markov kernel (Pk )k≥1 on (Rd , B(Rd )) given for all x ∈ Rd and A ∈ Rd
by
Z


1
2
2
Pk (x, A) =
)
dy .
(60)
exp
−
ky
−
h
(x)k
/(2σ
k
k
(2πσk2 )d/2 A
We denote for all n ≥ 1 by Qn the marginal laws of the sequence (Xk )k≥1 and given by
Q n = P1 · · · Pn .

(61)

In this section we compute an upper bound of {kδx Qn − δy Qn kTV , n ∈ 1} which does
not depend on the dimension d.
Let x, y ∈ Rd . We consider for all k ≥ 1 the following coupling (X1 , Y1 ) between
Pk (x, ·) and Pk (y, ·). Define the function E and e from Rd × Rd to Rd by
(
Ek (x, y)/ kEk (x, y)k if Ek (x, y) 6= 0
(62)
Ek (x, y) = hk (y) − hk (x) , ek (x, y) =
0
otherwise ,
For all x, y, z ∈ Rd , x 6= y, define


Fk (x, y, z) = hk (y) + Id −2ek (x, y)ek (x, y)T z
ϕσ2 (kEk (x, y)k − hek (x, y), zi)
,
αk (x, y, z) = k
ϕ σ2 (hek (x, y), zi)

(63)
(64)

k

where ϕ σ2 is the probability density of a zero-mean gaussian variable with variance σk2 .
k
Let Z1 be a standard d-dimensional Gaussian random variable. Set X1 = hk (x) + σk Z1
and
(
hk (y) + σk Z1
if Ek (x, y) = 0
Y1 =
B1 X1 + (1 − B1 ) Fk (x, y, Z1 ) if Ek (x, y) 6= 0 ,
where B1 is a Bernoulli random variable independent of Z1 with success probability
pk (x, y, z) = 1 ∧ αk (x, y, z) .
29

The construction above defines for all (x, y) ∈ Rd × Rd the Markov kernel Kk on (Rd ×
Rd , B(Rd ) ⊗ B(Rd )) given for all (x, y) ∈ Rd × Rd and A ∈ B(Rd ) ⊗ B(Rd ) by
Z
1D (hk (x), hk (y))
−kτk (x̃,x)k2 /(2σk2 )
Kk ((x, y), A) =
1
dx̃
(65)
(x̃,
x̃)e
A
(2πσk2 )d/2
Rd
Z
2
1Dc (hk (x), hk (y))
2
1A (x̃, x̃)pk (x, y, τk (x̃, x))e−kτk (x̃,x)k /(2σk ) dx̃
+
2
d/2
(2πσk )
Rd

Z
−kτk (x̃,x)k2 /(2σk2 )
+
1A (x̃, Fk (x, y, τk (x̃, x))) {1 − pk (x, y, τk (x̃, x))} e
dx̃ ,
Rd


where for all x̃ ∈ Rd , τk (x̃, x) = x̃ − hk (x) and D = (x̃, ỹ) ∈ Rd × Rd x̃ = ỹ . It is
shown in [4, Section 3.3] that for all x, y ∈ Rd and k ≥ 1, Kk ((x, y), ·) is a transference
plan of Pk (x, ·) and Pk (y, ·). Furthermore, we have for all x, y ∈ Rd and k ≥ 1


kEk (x, y)k
Kk ((x, y), D) = 2Φ −
.
(66)
2σk

eµ and E
e µ denote the
For all initial distribution µ0 on (Rd × Rd , B(Rd ) ⊗ B(Rd )), P
0
0
probability and the expectation respectively, associated with the sequence of Markov
kernels (Kk )k≥1 defined in (65) and µ0 on the canonical space ((Rd × Rd )N , (B(Rd ) ⊗
B(Rd ))⊗N ), {(Xi , Yi ), i ∈ N} denotes the canonical process and {Fei , i ∈ N} the corresponding filtration. Then if (X0 , Y0 ) = (x, y) ∈ Rd × Rd , for all k ≥ 1 (Xk , Yk ) is
coupling of δx Qk and δy Qk . Using Lindvall inequality, bounding kδx Qn − δy Qn kTV for
e(x,y) (Xn 6= Yn ). It is the content of the proof of the main
n ≥ 1 amounts to evaluate P
result of this section.
Theorem 23. Assume AR1. Then for all x, y ∈ Rd and n ≥ 1,
(
!)
kx
−
yk
kδx Qn − δy Qn kTV ≤ 1Dc ((x, y)) 1 − 2Φ −
,
1/2
2Ξn
where (Ξi )i≥1 is defined for all k ≥ 1 by Ξk =

Pk

Qi

2
i=1 σi {

We preface the proof by a technical Lemma.

j=1 (1

− ̟j )−2 }.

Lemma 24. For all ς, a > 0 and t ∈ R+ , the following identity holds
Z



ϕ 2 (t − y)
ϕ ς 2 (y) 1 − 1 ∧ ς
ϕ ς 2 (y)
R




|2y − t|
1 − 2Φ −
dy
2a



t
= 1 − 2Φ − 2
2(ς + a2 )1/2

Proof. The proof is postponed to [12, Section 4.2].

30



.

Proof of Theorem 23. Since for all k ≥ 1, (Xk , Yk ) is a coupling of δx Qk and δy Qk ,
e (x,y) (Xk 6= Yk ).
kδx Qk − δy Qk kTV ≤ P
P 2
Q
Define for all k1 , k2 ∈ N∗ , k1 ≤ k2 , Ξk1 ,k2 = ki=k
σ 2 { ij=k1 (1 − ̟j )−2 }. Let n ≥ 1.
1 i
We show by backward induction that for all k ∈ {0, · · · , n − 1},
e(x,y) (Xn 6= Yn )
P

e (x,y)
≤E

"

"

(

1Dc (Xk , Yk ) 1 − 2Φ −

kXk − Yk k

2 (Ξk+1,n )1/2

)##

, (67)

Note that the inequality for k = 0 will conclude the proof.
Since Xn 6= Yn implies that Xn−1 6= Yn−1 , the Markov property and (66) imply
e(x,y) (Xn 6= Yn ) = E
e (x,y)
P

h

i

e (X ,Y ) [1Dc (X1 , Y1 )]
1Dc (Xn−1 , Yn−1 )E
n−1
n−1

e (x,y)
≤E







kEn−1 (Xn−1 , Yn−1 )k
1Dc (Xn−1 , Yn−1 ) 1 − 2Φ −
2σn

Using AR1 and (62), kEn (Xn−1 , Yn−1 )k ≤ (1 − ̟n ) kXn−1 − Yn−1 k, showing (67) holds
for k = n − 1.
Assume that (67) holds for k ∈ {1, . . . , n − 1}. On {Xk 6= Yk }, we have
kXk − Yk k = − kEk (Xk−1 , Yk−1 )k + 2σk ek (Xk−1 , Yk−1 )T Zk ,
which implies


 kX − Y k 
1Dc (Xk , Yk ) 1 − 2Φ − k 1/2 k 

2Ξk+1,n 



 2σ e (X , Y )T Z − kE (X , Y )k 
k k
k−1
k−1
k
k
k−1
k−1
 .
= 1Dc (Xk , Yk ) 1 − 2Φ −
1/2


2Ξ


k+1,n

Since Zk is independent of Fek−1 , σk ek (Xk−1 , Yk−1 )T Zk is a real Gaussian random variable
with zero mean and variance σk2 , therefore by Lemma 24, we get
 



e Fek−1 1Dc (Xk , Yk ) 1 − 2Φ − kXk − Yk k 
E
(x,y)
1/2

2Ξk+1,n 
"




(

kEk (Xk−1 , Yk−1 )k
≤ 1Dc (Xk−1 , Yk−1 ) 1 − 2Φ −
1/2
2 σk2 + Ξk+1,n

)#

.

Using by AR1 that kEk (Xk−1 , Yk−1 )k ≤ (1 − ̟k ) kXk−1 − Yk−1 k concludes the induction.

31

Acknowledgements
The authors would like to thank Arnak Dalalyan for helpful discussions. The work
of A.D. and E.M. is supported by the Agence Nationale de la Recherche, under grant
ANR-14-CE23-0012 (COSMOS).

References
[1] J. H. Albert and S. Chib. Bayesian analysis of binary and polychotomous response
data. Journal of the American Statistical Association, 88(422):669–679, 1993.
[2] A. N. Borodin and P. Salminen. Handbook of Brownian motion—facts and formulae.
Probability and its Applications. Birkhäuser Verlag, Basel, second edition, 2002.
[3] S. Boucheron, G. Lugosi, and P. Massart. Concentration inequalities. Oxford University Press, Oxford, 2013. A nonasymptotic theory of independence, With a
foreword by Michel Ledoux.
[4] R. Bubley, M. Dyer, and M. Jerrum. An elementary analysis of a procedure for
sampling points in a convex body. Random Structures Algorithms, 12(3):213–235,
1998.
[5] Faes C., Ormerod J. T., and Wand M. P. Variational bayesian inference for parametric and nonparametric regression with missing data. Journal of the American
Statistical Association, 106(495):959–971, 2011.
[6] M. F. Chen and S. F. Li. Coupling methods for multidimensional diffusion processes.
Ann. Probab., 17(1):151–177, 1989.
[7] H. M. Choi and J. P. Hobert. The Polya-Gamma Gibbs sampler for Bayesian logistic
regression is uniformly ergodic. Electron. J. Statist., 7:2054–2064, 2013.
[8] N. Chopin and Ridgway J. Leave pima indians alone: binary regression as a benchmark for bayesian computation. Technical Report 0864, arXiv, une 2015.
[9] A. S. Dalalyan. Further and stronger analogy between sampling and optimization:
Langevin monte carlo and gradient descent. In Proceedings of the 30th Annual
Conference on Learning Theory.
[10] A. S. Dalalyan. Theoretical guarantees for approximate sampling from smooth and
log-concave densities. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), pages n/a–n/a, 2016.
[11] A. Durmus and É. Moulines. Non-asymptotic convergence analysis for the unadjusted langevin algorithm. Accepted for publication in Ann. Appl. Probab.
[12] A. Durmus and É. Moulines. Supplement to “high-dimensional bayesian inference
via the unadjusted langevin algorithm”, 2015. https://hal.inria.fr/hal-01176084/.
32

[13] A. Eberle. Quantitative contraction rates for Markov chains on continuous state
spaces. In preparation.
[14] D. L Ermak. A computer simulation of charged particles in solution. i. technique
and equilibrium properties. The Journal of Chemical Physics, 62(10):4189–4196,
1975.
[15] S. Frühwirth-Schnatter and R. Frühwirth. Data augmentation and MCMC for
binary and multinomial logit models statistical modelling and regression structures.
In Thomas Kneib and Gerhard Tutz, editors, Statistical Modelling and Regression
Structures, chapter 7, pages 111–132. Physica-Verlag HD, Heidelberg, 2010.
[16] R. B. Gramacy and N. G. Polson. Simulation-based regularized logistic regression.
Bayesian Anal., 7(3):567–590, 09 2012.
[17] U. Grenander. Tutorial in pattern theory. Division of Applied Mathematics, Brown
University, Providence, 1983.
[18] U. Grenander and M. I. Miller. Representations of knowledge in complex systems.
J. Roy. Statist. Soc. Ser. B, 56(4):549–603, 1994. With discussion and a reply by
the authors.
[19] T. E. Hanson, A J. Branscum, and W. O. Johnson. Informative g-priors for logistic
regression. Bayesian Anal., 9(3):597–611, 2014.
[20] C. C. Holmes and L. Held. Bayesian auxiliary variable models for binary and
multinomial regression. Bayesian Anal., 1(1):145–168, 03 2006.
[21] I. Karatzas and S.E. Shreve. Brownian Motion and Stochastic Calculus. Graduate
Texts in Mathematics. Springer New York, 1991.
[22] B. Klartag. A central limit theorem for convex sets. Invent. Math., 168(1):91–131,
2007.
[23] D. Lamberton and G. Pagès. Recursive computation of the invariant distribution
of a diffusion. Bernoulli, 8(3):367–405, 2002.
[24] D. Lamberton and G. Pagès. Recursive computation of the invariant distribution
of a diffusion: the case of a weakly mean reverting drift. Stoch. Dyn., 3(4):435–451,
2003.
[25] V. Lemaire. Estimation de la mesure invariante d’un processus de diffusion. PhD
thesis, Université Paris-Est, 2005.
[26] T. Lindvall and L. C. G. Rogers. Coupling of multidimensional diffusions by reflection. Ann. Probab., 14(3):860–872, 1986.

33

[27] J. C. Mattingly, A. M. Stuart, and D. J. Higham. Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise. Stochastic Process.
Appl., 101(2):185–232, 2002.
[28] R. M. Neal. Bayesian learning via stochastic dynamics. In Advances in Neural Information Processing Systems 5, [NIPS Conference], pages 475–482, San Francisco,
CA, USA, 1993. Morgan Kaufmann Publishers Inc.
[29] Y. Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Applied Optimization. Springer, 2004.
[30] G. Parisi. Correlation functions and computer simulations. Nuclear Physics B,
180:378–384, 1981.
[31] N. G. Polson, J. G. Scott, and J. Windle. Bayesian inference for logistic models using
Polya-Gamma latent variables. Journal of the American Statistical Association,
108(504):1339–1349, 2013.
[32] G. O. Roberts and R. L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341–363, 1996.
[33] P. J. Rossky, J. D. Doll, and H. L. Friedman. Brownian dynamics as smart Monte
Carlo simulation. The Journal of Chemical Physics, 69(10):4628–4633, 1978.
[34] D. Sabanés Bové and L. Held. Hyper-g priors for generalized linear models. Bayesian
Anal., 6(3):387–410, 2011.
[35] D. Talay and L. Tubaro. Expansion of the global error for numerical schemes solving
stochastic differential equations. Stochastic Anal. Appl., 8(4):483–509 (1991), 1990.
[36] C. Villani. Optimal transport : old and new. Grundlehren der mathematischen
Wissenschaften. Springer, Berlin, 2009.
[37] M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient langevin
dynamics. In Proceedings of the 28th International Conference on Machine Learning
(ICML-11), pages 681–688, 2011.
[38] J. Windle, N. G. Polson, and J. G. Scott. Bayeslogit: Bayesian logistic regression, 2013. http://cran.r-project.org/web/packages/BayesLogit/index.html R package version 0.2.

34

This figure "var_tv_0_v2.png" is available in "png" format from:
http://arxiv.org/ps/1605.01559v3

This figure "var_tv_1_2_v2.png" is available in "png" format from:
http://arxiv.org/ps/1605.01559v3

This figure "var_tv_1_4_v2.png" is available in "png" format from:
http://arxiv.org/ps/1605.01559v3

This figure "var_tv_3_4_v2.png" is available in "png" format from:
http://arxiv.org/ps/1605.01559v3

