Kinetic walks for sampling

TL;DR

This paper introduces and analyzes discrete-time kinetic walk algorithms, including a Zig-Zag sampler and hybrid jump/diffusion methods, for efficient sampling of complex distributions, extending continuous kinetic models to discrete settings.

Contribution

It defines and studies a discrete-space Zig-Zag sampler and hybrid kinetic samplers, bridging continuous models and discrete algorithms for sampling in high-dimensional spaces.

Findings

Discrete Zig-Zag sampler effectively samples from target distributions.

Hybrid jump/diffusion samplers handle multi-scale potentials efficiently.

The methods extend continuous kinetic processes to discrete frameworks.

Abstract

The persistent walk is a classical model in kinetic theory, which has also been studied as a toy model for MCMC questions. Its continuous limit, the telegraph process, has recently been extended to various velocity jump processes (Bouncy Particle Sampler, Zig-Zag process, etc.) in order to sample general target distributions on $\mathbb R^d$. This paper studies, from a sampling point of view, general kinetic walks that are natural discrete-time (and possibly discrete-space) counterparts of these continuous-space processes. The main contributions of the paper are the definition and study of a discrete-space Zig-Zag sampler and the definition and time-discretisation of hybrid jump/diffusion kinetic samplers for multi-scale potentials on $\mathbb R^d$.

Equations

$$(X_{k+1},V_{k+1})\ =\ \left\{\begin{array}[]{ll}(X_{k}+V_{k},V_{k})&\text{with probability }1-\alpha\\ (X_{k},-V_{k})&\text{else,}\end{array}\right.$$

$$(X_{k+1},V_{k+1})\ =\ \left\{\begin{array}[]{ll}(X_{k}+V_{k},V_{k})&\text{with probability }\min\left(\frac{\pi(X_{k}+V_{k})}{\pi(X_{k})},1\right)\\ (X_{k},-V_{k})&\text{else,}\end{array}\right.$$

$$(Y_{k+1},W_{k+1})\ =\ \left\{\begin{array}[]{ll}(Y_{k}+W_{k},-W_{k})&\text{with probability }\min\left(\frac{\pi(Y_{k}+W_{k})}{\pi(Y_{k})},1\right)\\ (Y_{k},W_{k})&\text{else.}\end{array}\right.$$

$$\frac{1}{n}\sum _{k=0}^{n-1}f(X_k,V_k)\underset{n\to +\infty }{⟶}\mu (f).$$

$$F(x)=\frac{1}{2}g(x)+\text{mathbbm}{1}_{x⩾1}\sum _{i=1}^{x-1}g(i)+\text{mathbbm}{1}_{x⩽-1}\sum _{i=x+1}^{-1}g(i),$$

$$\frac{1}{\sqrt{n}}\sum _{k=0}^{n-1}f(X_k,V_k)\stackrel{law}{\underset{n\to +\infty }{⟶}}\mathcal{N}(0,{\sigma }_f^2),$$

$$\mathbb{P}\left(T_1⩾k\right)=\prod _{j=0}^{k-1}{e}^{U(j)-U(j+1)}=e^{U(0)-U(k)},$$

$$A_n=\sum _{k=S_n}^{S_{n+1}-1}f(X_k,V_k).$$

$$\mathbb{E}\left(\left|\sum _{k=0}^{2T_1+1}f(X_k,V_k)\right|\right)$$

$$\mathbb{E}\left(A_0\right)=\sum _{k\in \mathbb{Z}}{e}^{U(0)-U(k)}\left(f(k,1)+f(k,-1)\right)=\lambda \mu (f)$$

$$\frac{1}{n}\sum _{j=0}^{K(n)}{A}_j⩽\frac{1}{n}\sum _{k=0}^{n-1}f(X_k,V_k)⩽\frac{1}{n}\sum _{j=0}^{K(n)+1}{A}_j.$$

$$\frac{K(n)}{n}\times \frac{S_{K(n)}}{K(n)}⩽1⩽\frac{K(n)+1}{n}\times \frac{S_{K(n)+1}}{K(n)+1}.$$

$$\frac{K(n)}{n}\times \frac{1}{K(n)}\sum _{j=0}^{K(n)}{A}_j⩽\frac{1}{n}\sum _{k=0}^{n-1}f(X_k,V_k)⩽\frac{K(n)+1}{n}\times \frac{1}{K(n)+1}\sum _{j=0}^{K(n)+1}{A}_j,$$

$$\sqrt{\frac{\lambda }{n}}\sum _{j=0}^{K(n)}{A}_j\stackrel{law}{\underset{n\to \infty }{⟶}}\mathcal{N}(0,{\sigma }_A^2),$$

$$\frac{1}{\sqrt{n}}\left(\sum _{k=0}^{n-1}f(X_k,V_k)-\sum _{k=0}^{K(n)}{A}_k\right)\stackrel{\mathbb{P}}{\underset{n\to \infty }{⟶}}0,$$

$$\frac{1}{\sqrt{n}}\sum _{k=0}^{n-1}f(X_k,V_k)\stackrel{law}{\underset{n\to \infty }{⟶}}\mathcal{N}(0,{\sigma }_A^2/\lambda ).$$

$$A_0^{\prime }=\sum _{k=1}^{2T_1}f(X_k,V_k),A_0^{\prime \prime }=\sum _{k=2T_1+2}^{S_1-1}f(X_k,V_k),$$

$$\mathbb{E}\left((A_0^{\prime }{)}^2\right)$$

$$\frac{1}{\lambda }\mathbb{E}\left(A_0^2\right)$$

$${\left(\frac{1}{\sqrt{n}}\sum _{k=0}^{\lfloor nt\rfloor }\left(f(X_k,V_k)-\mu (f)\right)\right)}_{t⩾0}\stackrel{law}{\underset{n\to +\infty }{⟶}}({\sigma }_f{B}_t{)}_{t⩾0}$$

$${\tau }_{\epsilon }=\inf \{n\in \mathbb{N},X_n^{\epsilon }\notin ]]a,b[[\}.$$

$$\mathbb{E}\left({\tau }_{\epsilon }\right)$$

$$\mathbb{P}(X_{{\tau }_{\epsilon }}^{\epsilon }=a)\underset{\epsilon \to 0}{⟶}\frac{1}{2}\left(1+\text{mathbbm}{1}_{U(a)⩽U(b)}-\text{mathbbm}{1}_{U(b)⩽U(a)}\right).$$

$$2(\beta -\alpha +1)\left(G_{\epsilon }-1\right)\stackrel{sto}{⩽}{\tau }_{\epsilon }\stackrel{sto}{⩽}2(b-a+1)G_{\epsilon }$$

$$e^{-U(b)/\epsilon }+e^{-U(a)/\epsilon }-e^{-\left(U(b)+U(a)\right)/\epsilon }.$$

$$p:=\mathbb{P}\left(\tau <S_1\right)$$

$$\mathbb{P}\left(X_{\tau }=b\right)=\mathbb{P}\left(X_{\tau }=b|\tau <S_1\right)=\frac{1}{p}{e}^{-U(b)/\epsilon },$$

$$\mathbb{E}(\tau )$$

$$\left|\mathbb{E}(\tau -S_{K-1})\right|⩽2(b-a)=e^{E_1/\epsilon }\underset{\epsilon \to 0}{\mathcal{O}}\left(e^{-E_2/\epsilon }\right),$$

$$\mathbb{P}\left(S_1=2(\beta -\alpha +1)\right)$$

Full text

1 Introduction

The classical persistent walk on $\mathbb{Z}$ is the Markov chain $(X_{k},V_{k})_{k\in\mathbb{N}}$ on $\mathbb{Z}\times\{-1,1\}$ with transitions

[TABLE]

for some $\alpha\in[0,1]$. It describes the constant-speed motion of a self-propelled particle, $X_{k}$ denoting the position of the particle and $V_{k}$ its velocity. Since the time between two changes of the velocity follows a geometric distribution with parameter $\alpha$, $(\alpha X_{\lfloor\alpha t\rfloor},V_{\lfloor\alpha t\rfloor})_{t\geqslant 0}$ naturally converges as $\alpha$ vanishes to the so-called telegraph process, for which the flips of the velocity are governed by a Poisson process [25]. From the seminal work of Goldstein [19], these two processes, and various extensions, have been studied in details, in particular from the point of view of statiscial physics and kinetic theory, or for other modelling motivations in physics, finance or biology (see for instance [26, 50, 10, 22, 1, 45, 21] and references within).

Meanwhile, the search for efficient Markov Chain Monte Carlo (MCMC) methods led to the development of so-called rejection-free or lifted chains (see e.g. [28, 12, 2] and references within). In this context, the persistent walk has been a toy model to understand the efficiency of these algorithms, especially when compared to the reversible simple walk [12, 13, 38]. For instance, correctly scaled, the persistent walk shows a ballistic behaviour, which means its expected distance to its initial position after $K$ steps is of order $K$, while the simple walk shows a diffusive behaviour, moving to a distance $\sqrt{K}$ after $K$ steps. Since the efficiency of the MCMC schemes is related to the speed at which the space is explored, this is an argument in favour of non-reversible kinetic processes. The model being simple, it is even possible to determine the optimal $\alpha$ (in the sense that it gives the maximal rate of convergence toward equilibrium on the periodic torus $\mathbb{Z}/(N\mathbb{Z})$), which turns to be of order $1/N$ [38]. This is consistent, as $N$ goes to infinity, with the ballistic scaling that yields the telegraph process (by contrast, if $\alpha$ is constant with $N$ and if time is accelerated by $N^{2}$, the persistent walk converges to the Brownian motion).

Of course, both the persistent walk and the telegraph process sample the uniform measure in dimension one, which is not of practical interest. These last years, the telegraph process has been extended to several continuous-space processes, such as the Zig-Zag sampler [4, 3, 6, 5] or the Bouncy Particle Sampler [43, 39, 14, 8], which are velocity jump processes designed to target any given distribution in any dimension. Many variants like randomized bounces [49, 36] are currently being developped and we refer to the review [48] for more details, considerations and references on this vivid topic.

The present paper is concerned with similar extensions, but conducted at the level of the persistent walk rather than of the continuous kinetic process. Or, from another viewpoint, we are interested in persistent walks, but through the prism of MCMC sampling rather than kinetic theory. The motivations are the following: first, the discrete chain yields some insights on their continuous-time limits (for instance, we will see that the Zig-Zag process can be seen as the continuous limit of a Gibbs algorithm). Second, used in an MCMC scheme on $\mathbb{Z}^{d}$, a persistent walk shares, as will be detailed in this work, the following advantages with its continuous counterparts: thinning, factorization and ballistic behaviour. Finally, although continuous-time velocity jump processes can sometimes be sampled exactly thanks to thinning methods, it is not necessarily the case for mixed diffusion/jump kinetic samplers (see Section 5), in which case the corresponding chain obtained through an integration scheme (say, Euler scheme) is a discrete-time kinetic walk.

The rest of the paper is organized as follows. We start in Section 2 with the definition of an analogous on $\mathbb{Z}$ of the Zig-Zag process on $\mathbb{R}$ (or, equivalently, of the persistent walk but in a general potential landscape). The simplicity of the chain allows an elementary study of its ergodicity, of its metastable behaviour at small temperature via an Eyring-Kramers formula and of its convergence toward the continuous Zig-Zag process on $\mathbb{R}$ under proper scaling. Section 3 is a general and informal discussion about kinetic walks on $\mathbb{R}^{d}$, their simulation, invariant measures and continuous-time scaling limits. Finally, the last two sections present two particular applications, which are the main contributions of this work: the discrete Zig-Zag walk in a general potential in Section 4 and hybrid jump/diffusion kinetic samplers with a numerical integrator in Section 5. Although related by their motivation (the understanding of sampling with kinetic walks), the four sections are sufficiently independent to be read separately one from the other. The definition of the Zig-Zag walk associated to a general potential and all the results on this topic in Section 2 and Section 4 are new. The specific numerical scheme introduced in Section 5 is also new, although straightforwardly obtained from the well-known general method of Strang splittings.

Notations. If $x,y\in\mathbb{R}^{d}$, we denote $x\cdot y$ their scalar product and $|x|=\sqrt{x\cdot x}$. The Dirac mass at $x$ is denoted $\delta_{x}$ and $\mathbbm{1}_{A}$ is 1 if $A$ and 0 else. For $r\in\mathbb{R}$, $(r)_{+}=\max(r,0)$. The set of $\mathcal{C}^{k}$ functions on $\mathbb{R}^{d}$ with compactly supported supported is denoted $\mathcal{C}^{k}_{c}(\mathbb{R}^{d})$. The Gaussian distribution on $\mathbb{R}^{d}$ with mean $m$ and variance $\Sigma^{2}$ is denoted $\mathcal{N}(m,\Sigma^{2})$. We denote respectively $\mathcal{P}(E)$, $\mathcal{M}(E)$ and $\mathcal{M}_{b}(E)$ the sets of probability measures, measurable functions and bounded measurable functions on a measurable space $E$, and for $\mu\in\mathcal{P}(E)$ and $f\in L^{1}(\mu)$ we write $\mu f=\mu(f)=\int f\text{d}\mu$. When $(X_{t}^{\varepsilon})_{t\geqslant 0}$ for $\varepsilon>0$ and $(Y_{t})_{t\geqslant 0}$ are cádlág processes on $\mathbb{R}^{d}$, we write $(X_{t}^{\varepsilon})_{t\geqslant 0}\overset{law}{\underset{\varepsilon\rightarrow 0}{\longrightarrow}}(Y_{t})_{t\geqslant 0}$ for the convergence in law in the Skorohod topology. Recall a sequence $(x_{n})_{n\in\mathbb{N}}$ of càdlàg functions from $\mathbb{R}_{+}$ to $\mathbb{R}^{d}$ is said to converge to $x$ if, on all finite time interval, it converges uniformly up to a uniformly small change of time, i.e. if there exists a sequence $(\gamma_{n})_{n\in\mathbb{N}}$ with $\gamma_{n}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ increasing so that $\sup_{t\in[0,T]}(|x_{n}(\gamma_{n}(t))-x(t)|+|\gamma_{n}(t)-t|)\rightarrow 0$ as $n\rightarrow\infty$ for all $T>0$.

2 The Zig-Zag walk on $\mathbb{Z}$

Let $U:\mathbb{Z}\rightarrow\mathbb{R}$ be such that $\mathcal{Z}=\sum_{x\in\mathbb{Z}}\exp(-U(x))<+\infty$, $\pi(x)=\exp(-U(x))/\mathcal{Z}$ be the associated Gibbs distribution and $\mu(x,v)=\pi(x)/2$ for $v\in\{-1,1\}$ and $x\in\mathbb{Z}$. We consider the Markov chain $(X_{k},V_{k})_{k\in\mathbb{N}}$ on $\mathbb{Z}\times\{-1,1\}$ with transitions

[TABLE]

which we call the Zig-Zag walk on $\mathbb{Z}$. This transition can be seen as the composition of two Markov transitions. Indeed, consider on $\mathbb{Z}\times\{-1,1\}$ the Markov kernel given by $x,v\mapsto\delta_{(x+v,-v)}$. Since $(y,w)=(x+v,-v)$ implies that $(x,v)=(y+w,-w)$, this kernel is symmetric. If a Metropolis accept/reject step with target measure $\mu$ is added, the transition of the resulting chain is simply

[TABLE]

By construction of the Metropolis-Hastings algorithm, this transition leaves $\mu$ invariant. Now if we compose this transition with the deterministic transition $(Y_{k+1},W_{k+1})=(Y_{k},-W_{k})$, which obviously leaves $\mu$ invariant, we obtain the initial chain. We have thus obtained that $\mu$ is invariant for the Zig-Zag walk. Note however that, although both intermediate transition kernels are reversible with respect to $\mu$, their composition is not. Indeed, $\mathbb{P}((X_{2},V_{2})=(X_{0},V_{0}))=0$ for all initial condition.

The chain is clearly irreducible, and it is periodic. Indeed, $(-1)^{X_{k}}V_{k}=(-1)^{k+X_{0}}V_{0}$ for all $k\in\mathbb{N}$, so that, if $X_{k}$ is even with $V_{k}=1$ or $X_{k}$ is odd with $V_{k}=-1$, then $X_{k+1}$ is odd with $V_{k+1}=1$ or $X_{k}$ is even with $V_{k+1}=-1$. In particular, the period is even. If $U$ admits a strict local minimum $x_{0}$ then there is a path of length 2 with strictly positive probability from $(x_{0},1)$ to itself (which is $(x_{0},1)\rightarrow(x_{0},-1)\rightarrow(x_{0},1)$), so that the period is exactly 2, but this may not be the case in general. For instance, with $U(k)=|\lfloor k/2\rfloor|$, the reader can check that the period is 4.

In the following, we prove an ergodic Law of Large Number and a Central Limit Theorem (CLT) in Theorem 1 (in the unimodal case; other cases are treated in any dimension in Section 4.4.4), an Eyring-Kramers formula in Theorem 2 and the convergence toward the continuous Zig-Zag process in Theorem 3.

2.1 Asymptotic results

Although the chain is already quite simple, let us focus for now on the case where $\pi$ is unimodal. In that case, and similarly to the continuous-time case [3], ergodicity can be established through elementary considerations on renewal chains.

Theorem 1**.**

Suppose that $U$ is decreasing on $\rrbracket-\infty,0\rrbracket$ and increasing on $\llbracket 0,+\infty\llbracket$, and let $f\in L^{1}(\mu)$. Then, for all initial conditions, almost surely,

[TABLE]

Moreover, denoting $g(x)=f(x,1)+f(x,-1)$ and

[TABLE]

suppose that $M_{f}:=\mathbb{E}_{\pi}(g(X)F(X))<\infty$ and that $\mu(f)=0$. Then

[TABLE]

with some explicit variance $\sigma_{f}^{2}\leqslant 3M_{f}$.

Proof.

Consider first the case where $(X_{0},V_{0})=(0,1)$ and denote $T_{1}=\inf\{n\in\mathbb{N}\ :\ X_{n+1}=X_{n}\}$, $T_{2}=\inf\{n\in\mathbb{N}\ :\ X_{2T_{1}+n+2}=X_{2T_{1}+1+n}\}$. The monotonicities of $U$ implies that almost surely $X_{n}$ increases for $n\in\llbracket 0,T_{1}\rrbracket$, decreases for $n\in\llbracket T_{1}+1,2T_{1}+1+T_{2}\rrbracket$ with $X_{2T_{1}+1}=0$, and finally, denoting $S_{1}=2(T_{1}+T_{2}+1)$, increases for $n\in\llbracket 2T_{1}+2+T_{2},S_{1}\rrbracket$ with $(X_{S_{1}},V_{S_{1}})=(0,1)$ (cf. Fig. 1). Remark that

[TABLE]

so that $T_{1}<\infty$ almost surely (since we assumed that $\mathcal{Z}<\infty$, $U$ necessarily goes to $\infty$ at $\infty$). The same goes for $T_{2}$, hence for $S_{1}$. By the strong Markov property, $(X_{n},V_{n})_{n\geqslant S_{1}}$ has the same law as $(X_{n},V_{n})_{n\in\mathbb{N}}$ and is independent from $(X_{n},V_{n})_{n\in\llbracket 0,S_{1}-1\rrbracket}$. Denote $S_{0}=0$ and, for all $n\in\mathbb{N}$, $S_{n+1}=\inf\{k>S_{n}\ :\ (X_{k},V_{k})=(0,1)\}$ and, given a function $f\in L^{1}(\mu)$,

[TABLE]

The $A_{n}$’s are i.i.d. and

[TABLE]

The sum for $k\in\llbracket 2T_{1}+2,S_{1}-1\rrbracket$ is treated the same way, so that $\mathbb{E}|A_{0}|<\infty$ and

[TABLE]

with $\lambda=2e^{U(0)}\mathcal{Z}$. The proof then follows from classical renewal arguments, which we recall for completeness. Considering the case $f=1$, the law of large numbers implies that $S_{n}/n$ converges almost surely toward $\lambda$ as $n$ goes to infinity. For $n\in\mathbb{N}$ set $K(n)=\sup\{k\in\mathbb{N}\ :\ S_{k}\leqslant n\}$. If $f$ is positive then for all $n\in\mathbb{N}$,

[TABLE]

Applied with $f=1$, this reads

[TABLE]

Since $K(n)$ almost surely goes to infinity with $n$, we get that $K(n)/n$ almost surely converges to $1/\lambda$. Applied again with a general positive $f$, now,

[TABLE]

and letting $n$ go to infinity concludes. If $f$ is not positive, the same conclusion follows from the decomposition $f=(f)_{+}-(-f)_{+}$.

Now, consider the case of any general initial condition $(X_{0},V_{0})=(x,v)$, and let $R=\inf\{n\in\mathbb{N}\ :\ (X_{n},V_{n})=(0,1)\}$. By similar arguments as above, $R<\infty$ almost surely so that $n^{-1}\sum_{k<R}f(X_{k},V_{k})$ almost surely goes to zero, while by the Markov property, $n^{-1}\sum_{R\leqslant k\leqslant n}f(X_{k},V_{k})$ converges toward $\mu(f)$, which concludes.

The proof of the CLT is similar, and we refer to [3, Lemma 4] to get that, if $\mu(f)=0$,

[TABLE]

provided that $\sigma_{A}^{2}:=\mathbb{E}\left(A_{0}^{2}\right)<\infty$. Now, even if $(X_{0},V_{0})\neq(0,1)$, as before,

[TABLE]

hence, provided that $\sigma_{A}^{2}<\infty$,

[TABLE]

Decompose $A_{0}=f(0,1)+A_{0}^{\prime}+f(0,-1)+A_{0}^{\prime\prime}$ with

[TABLE]

and remark that by the Markov property, $A_{0}^{\prime\prime}$ is independent from $A_{0}$. Compute

[TABLE]

The case of $A_{0}^{\prime\prime}$ is similar and, using that $1/\lambda=\pi(0)/2$, we get

[TABLE]

which concludes. In fact $\sigma_{f}^{2}=\sigma_{A}^{2}/\lambda$ can be computed since $\mathbb{E}(A_{0}^{\prime})=\lambda\sum_{x\in\mathbb{N}_{*}}g(x)\pi(x)$, and similarly for $A_{0}^{\prime\prime}$. ∎

For $N\in\mathbb{N}_{*}$ and $t_{1}<\dots<t_{N}$, considering $K(t_{i})=\sup\{k\in\mathbb{N}\ :\ S_{k}\leqslant\lfloor nt\rfloor\}$ and decomposing $\sum_{n=1}^{K(t_{N})-1}A_{n}=\sum_{i=1}^{N}\sum_{n=K(t_{i})}^{K(t_{i+1})-1}A_{n}$, the previous elementary proof is easily extended to obtain a functional CLT, namely the convergence

[TABLE]

where $(B_{t})_{t\geqslant 0}$ is a one-dimensional Brownian motion. See also Section 4.4.4.

2.2 Metastability

We now consider the question of escape times from local minima at low temperature, as in [39] for the Zig-Zag process on $\mathbb{R}$. Recall that a random variable $G$ on $\mathbb{R}$ is said to be stochastically larger than a random variable $F$ on $\mathbb{R}$ if $\mathbb{P}\left(G<t\right)\leqslant\mathbb{P}\left(F<t\right)$ for all $t\in\mathbb{R}$. In that case we write $F\overset{sto}{\leqslant}G$.

Theorem 2**.**

Let $U:\mathbb{Z}\rightarrow\mathbb{R}$ and, for all $\varepsilon>0$, let $(X^{\varepsilon}_{k},Y^{\varepsilon}_{k})_{k\in\mathbb{N}}$ be the persistent walk on $\mathbb{Z}$ associated to $U/\varepsilon$ and with initial condition $(0,1)$. Suppose that $U(0)=0$ and that $U$ is decreasing on $\rrbracket-\infty,0\rrbracket$ and increasing on $\llbracket 0,+\infty\llbracket$. Let $a<\alpha\leqslant 0\leqslant\beta<b$ be such that $\llbracket\alpha,\beta\rrbracket=\{k\in\mathbb{Z},\ U(k)=0\}$, and let

[TABLE]

Then

[TABLE]

with $E_{1}=\min\left(U(a),U(b)\right)$, $E_{2}=\min\left(U(\alpha-1),U(\beta+1)\right)$ and $E_{3}=|U(a)-U(b)|$. Moreover, $\tau_{\varepsilon}/\mathbb{E}(\tau_{\varepsilon})$ converges in law as $\varepsilon$ vanishes to an exponential random variable with parameter 1, and

[TABLE]

Finally, for all $\varepsilon>0$,

[TABLE]

where $G_{\varepsilon}$ is a geometric random variable with parameter

[TABLE]

Proof.

The proof is similar to the one appearing in [39]. To alleviate notations, we only write $X_{k},V_{k}$ and $\tau$ for $X_{k}^{\varepsilon},V_{k}^{\varepsilon}$ and $\tau_{\varepsilon}$. Like in the previous proof, set $S_{0}=0$ and, by induction, $S_{n+1}=\inf\{k>S_{n},\ (X_{k},V_{k})=(0,1)\}$. For all $n\in\mathbb{N}$, let $\tilde{S}_{n}=\inf\{k>S_{n},\ (X_{k},V_{k})=(0,-1)\}$, and let $K=\inf\{n\in\mathbb{N},\ S_{n}>\tau\}$. Keep Figure 1 in mind. By the strong Markov property, $K$ follows a geometric distribution with parameter

[TABLE]

and

[TABLE]

which indeed converges as $\varepsilon$ vanishes to 0 if $U(a)<U(b)$, 1 if $U(a)>U(b)$ and $1/2$ if $U(a)=U(b)$. Decomposing $\tau\ =\ \tau-S_{K-1}+\sum_{i=1}^{K-1}(S_{i}-S_{i-1})$, remark that almost surely $2(\beta-\alpha+1)\leqslant S_{i}-S_{i-1}\leqslant 2(b-a+1)$ for all $i<K$ and $\tau-S_{K-1}<2(b-a+1)$, which proves the last claim of the theorem. Besides, again by the strong Markov property, conditionally to $K$, $(S_{i}-S_{i-1})_{i\in\llbracket 1,K-1\rrbracket}$ are i.i.d. random variables independent from $K$, so that

[TABLE]

Since $0\leqslant\tau-S_{K-1}\leqslant 2(b-a)$ almost surely,

[TABLE]

where we used that $E_{1}\geqslant E_{2}$. For $\varepsilon$ small, the most likely trajectory of the process between times [math] and $S_{1}$ is the following: starting from $(0,1)$, it deterministically goes to $(\beta,1)$, then jumps to $(\beta,-1)$ with high probability, then deterministically goes to $(\alpha,-1)$ and jumps to $(\alpha,1)$ with high probability before going back to $(0,1)$ deterministically. More precisely, $S_{1}\geqslant 2(\beta-\alpha+1)$ almost surely, and

[TABLE]

On the other hand, conditionally to $\tau>S_{1}$, almost surely, $S_{1}\leqslant 2(b-a)$, so that

[TABLE]

Thus, we get that

[TABLE]

where we used again that $E_{1}\geqslant E_{2}$. Using in (2) this estimate together with (3) and the fact that $\mathbb{E}(K)=1/p$ concludes the proof of the Eyring-Kramers formula (1).

Finally, $K$ being an exponential random variable whose parameter vanishes with $\varepsilon$, $pK$ converges in law toward an exponential random variable with parameter 1. By the Markov inequality, for any $\delta>0$,

[TABLE]

Since $\mathbb{E}\left(S_{1}\ |\ \tau>S_{1}\right)$ and $p\mathbb{E}(\tau)$ both converges toward $2(\beta-\alpha+1)$ as $\varepsilon$ vanishes,

[TABLE]

From Slutsky’s theorem, $\sum_{i=1}^{K-1}(S_{i}-S_{i-1})/\mathbb{E}(\tau)$ converges in law to an exponential random variable with parameter 1 as $\varepsilon\rightarrow 0$. Finally, $|(\tau-S_{K-1})/\mathbb{E}(\tau)|\leqslant 2(b-a)/\mathbb{E}(\tau)$ almost surely goes to zero, which concludes.

∎

The simulated annealing chain obtained by considering a non-constant temperature $(\varepsilon_{k})_{k\in\mathbb{N}}$ and a potential $U$ with possibly several local minima could also be studied by similar arguments as in [39, Theorem 3.1] to get a necessary and sufficient condition on the cooling schedule for convergence in probability toward the global minima of $U$.

2.3 Continuous scaling limit

The (continuous-time) Zig-Zag process on $\mathbb{R}$ associated to a potential $H\in\mathcal{C}^{1}(\mathbb{R})$ (also known as the (integrated) telegraph or run-and-tumble process) is the Markov process on $\mathbb{R}\times\{-1,1\}$ with generator

[TABLE]

In other words, it is a piecewise deterministic Markov process that, starting from an initial condition $(y,w)$, follows the flow $(Y_{t},W_{t})=(y+tw,w)$ up to a random time $T$ with distribution $\mathbb{P}(T>t)=\exp(-\int_{0}^{t}(wH^{\prime}(y+sw))_{+}\text{d}s)$, at which point $(Y_{T},W_{T})=(y+Tw,-w)$, after which it follows again the deterministic flow up to a new random jump time, etc.

Theorem 3**.**

For $H\in\mathcal{C}^{2}(\mathbb{R})$ that goes to infinity at infinity, for all $\varepsilon>0$, define $U_{\varepsilon}:\mathbb{Z}\mapsto\mathbb{R}$ by $U_{\varepsilon}(k)=H(\varepsilon k)$ for all $k\in\mathbb{Z}$. Let $(X^{\varepsilon}_{k},V^{\varepsilon}_{k})_{k\in\mathbb{N}}$ be the persistent walk on $\mathbb{Z}$ associated to $U_{\varepsilon}$ and with some initial condition $(x_{0}^{\varepsilon},v_{0})$. Suppose that $\varepsilon x_{0}^{\varepsilon}$ converges to some $x_{0}^{*}\in\mathbb{R}$ as $\varepsilon$ vanishes. Then,

[TABLE]

where $(Y_{t},W_{t})_{t\geqslant 0}$ is a Zig-Zag process on $\mathbb{R}$ associated to $H$ and with $(Y_{0},W_{0})=(x_{0}^{*},v_{0})$.

Proof.

Denote $T^{\varepsilon}_{1}=\varepsilon\inf\{n\in\mathbb{N}\ :\ X_{n+1}^{\varepsilon}=X_{n}^{\varepsilon}\}$. Its cumulative function is

[TABLE]

From

[TABLE]

we get that $T_{1}^{\varepsilon}$ converges in law as $\varepsilon$ vanishes to a random variable $T_{1}^{0}$ with cumulative function

[TABLE]

Remark that

[TABLE]

so that $T_{1}^{0}$ is almost surely finite, and similarly for $T_{1}^{\varepsilon}$ for all $\varepsilon>0$. In particular,

[TABLE]

Let $(A_{j})_{j\in\mathbb{N}}$ be an i.i.d. sequence of random variables uniformly distributed over $[0,1]$. For all $\varepsilon\geqslant 0$, set $(Z_{0}^{\varepsilon},R_{0}^{\varepsilon},S_{0}^{\varepsilon})=(\varepsilon x_{0}^{\varepsilon},v_{0},0)$ (with, in the case where $\varepsilon=0$, $Z_{0}^{0}=x_{0}^{*}$). Suppose by induction that, for some $n\in\mathbb{N}$, $(Z_{n}^{\varepsilon},R_{n}^{\varepsilon},S_{n}^{\varepsilon})$ has been defined for all $\varepsilon\geqslant 0$ and is independent from $(A_{j})_{j\geqslant n}$. Then, for all $\varepsilon\geqslant 0$, set

[TABLE]

and $R_{n+1}^{\varepsilon}=-R_{n}^{\varepsilon}$. Remark that, for all $t\geqslant 0$, $(x,v)\mapsto F_{x,v}^{\varepsilon}(t)$ is continuous, uniformly in $\varepsilon$. As a consequence, for all $N\in\mathbb{N}$, by the previous result, almost surely,

[TABLE]

By construction, for all $\varepsilon>0$,

[TABLE]

with $T_{0}^{\varepsilon}=0$ and, by induction, $T^{\varepsilon}_{k+1}=\varepsilon\inf\{k>T^{\varepsilon}_{k}\ :\ X_{n+1}^{\varepsilon}=X_{n}^{\varepsilon}\}$, and similarly

[TABLE]

with $J_{0}=0$ and by induction $J_{k+1}=\inf\{t>J_{k},\ W_{t}=-W_{J_{k}}\}$. At this point, we have thus proved that the skeleton chain of the persistent walk (namely the persistent walk observed at its jump times, and those jump times) converges in law toward the skeleton chain of the Zig-Zag process (namely the process observed at its jump times, and those jump times). The convergence of the full chain is then a consequence from the fact that the latter is a deterministic function of its skeleton chain, as we detail now.

Note that $(S_{n}^{0})_{n\in\mathbb{N}}$ has the same distribution as $(J_{n})_{n\geqslant 0}$. Moreover, for any $t\geqslant 0$ and for all $s\in[0,t]$, $|Y_{s}-x_{0}^{*}|\leqslant t$, so that the jump rate of the Zig-Zag process is bounded for times $s\in[0,t]$ by $\omega(t)=\sup_{x\in[x_{0}^{*}-t,x_{0}^{*}+t]}|H^{\prime}(x)|$, which is finite. In particular $\sup\{n\in\mathbb{N},\ J_{n}<t\}$ the number of jumps of the Zig-Zag process on $[0,t]$ is stochastically smaller than a Poisson process with rate $\omega(t)$, hence is almost surely finite. As a consequence, almost surely $S_{n}^{0}\rightarrow+\infty$ as $n\rightarrow\infty$, and similarly for $S_{n}^{\varepsilon}$ for all $\varepsilon\geqslant 0$.

Now the continuous-time processes are obtained by interpolating the skeleton chains. For all $\varepsilon\geqslant 0$ and $n\in\mathbb{N}$, set $\tilde{R}_{t}^{\varepsilon}=R_{n}^{\varepsilon}$ for all $t\in[S_{n}^{\varepsilon},S_{n+1}^{\varepsilon}[$. For all $t\geqslant 0$, set $\tilde{Z}_{t}^{0}=x_{0}^{*}+\int_{0}^{t}\tilde{R}_{s}^{0}\text{d}s$ and remark that, by construction, $\tilde{Z}_{T_{n}}^{0}=Z_{n}^{0}$ for all $n\in\mathbb{N}$. Finally, for all $\varepsilon>0$, all $n\in\mathbb{N}$ and all $k\in\llbracket T_{n}/\varepsilon,T_{n+1}/\varepsilon\rrbracket$, set

[TABLE]

and $\tilde{Z}_{t}^{\varepsilon}=\tilde{Z}_{k\varepsilon}^{\varepsilon}$ for all $t\in[k\varepsilon,(k+1)\varepsilon[$. This construction ensures that, for all $\varepsilon>0$,

[TABLE]

For all $\varepsilon>0$, consider the increasing continuous change of time $\gamma_{\varepsilon}:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ given by: for all $n\in\mathbb{N}$, $\gamma_{\varepsilon}(T_{n}^{\varepsilon})=T_{n}^{0}$ and $\gamma_{\varepsilon}$ is linear on $[T_{n}^{\varepsilon},T_{n+1}^{\varepsilon}]$. In particular, $\tilde{R}_{\gamma_{\varepsilon}(t)}^{\varepsilon}=R_{t}^{0}$ for all $t\geqslant 0$ (they start at the same value and both change sign at all times $(T_{n}^{0})_{n\geqslant 1}$). Moreover, from the convergence of the skeleton chains, $\sup_{t\in[0,T_{n}^{0}]}\left(|\gamma_{\varepsilon}(t)-t|+|\tilde{Z}_{\gamma_{\varepsilon}(t)}^{\varepsilon}-Z_{t}^{0}|\right)$ almost surely goes to [math] as $\varepsilon$ vanishes. Together with the fact $T_{n}^{0}$ almost surely goes to $+\infty$, this means that for all fixed $T>0$, $\sup_{t\in[0,T]}\left(|\gamma_{\varepsilon}(t)-t|+|\tilde{Z}_{\gamma_{\varepsilon}(t)}^{\varepsilon}-Z_{t}^{0}|\right)$ almost surely goes to [math], which concludes.

∎

3 Kinetic walks on $\mathbb{R}^{d}$

In the following, we will be interested in a class of Markov chains $(X_{n},V_{n})_{n\in\mathbb{N}}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ for $d\in\mathbb{N}_{*}$, with transitions given by

[TABLE]

for some $\delta>0$ and some kernel $p:(x,v)\in\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto p(x,v;\cdot)\in\mathcal{P}(\mathbb{R}^{d})$. We call such a chain the kinetic walk on $\mathbb{R}^{d}$ associated to $p$ with timestep $\delta$. Up to a rescaling of the velocities, we can always consider that $\delta=1$.

This definition is close to – but distinct from – the definition of second-order Markov chains on $\mathbb{R}^{d}$ (sometimes also called correlated random walks like in [21]). Indeed, $(X_{n},X_{n-1})_{n\in\mathbb{N}}$ is a Markov chain if and only if $(X_{n},X_{n}-X_{n-1})_{n\in\mathbb{N}}$ is, with a simple way to express the transition of one of these chains from the transition of the other. Denoting $V_{n}=(X_{n}-X_{n-1})/\delta$ would yield $X_{1}=X_{0}+\delta V_{0}$. On the contrary, consider the chain defined in Section 2, which satisfies (5) with $\delta=1$. For this chain, $(X_{n},X_{n-1})_{n\in\mathbb{N}}$ is not Markovian: if $X_{n}=X_{n-1}$ is at a strict local minimum of the potential $U$, then that only means that the velocity $V_{n}$ has changed between times $n-1$ and $n$, but it could be from $1$ to $-1$ or the converse (which we could know by looking farther in the past trajectory, for instance with the fact that $(-1)^{X_{k}}V_{k}=(-1)^{k+X_{0}}V_{0}$ for all $k\in\mathbb{N}$), and this affects the law of $X_{n+1}$. Our present definition is only motivated by the fact it gives a simple and unified framework for the cases studied in Sections 2, 4 and 5. Second-order Markov or related chains (like the discrete-time bounce sampler of [49]) may be studied with the same arguments (especially concerning their continuous-time scaling limits). We use the term kinetic rather than persistent in order to keep the latter for cases where the velocity is typically constant for some times, and this is not always the case for the different kinetic walks we will be interested in.

Note that discrete-space walks can be seen as particular cases of walks on $\mathbb{R}^{d}$ as follows. Let $(X_{n},V_{n})_{n\in\mathbb{N}}$ be a kinetic walk on $\mathbb{Z}^{d}\times\mathbb{Z}^{d}$ with transitions given by (5) with $\delta>0$ and $p:\mathbb{Z}^{d}\times\mathbb{Z}^{d}\rightarrow\mathcal{P}(\mathbb{Z}^{d})$ and let $\eta,\kappa>0$. Consider the chain $(\tilde{X}_{n},\tilde{V}_{n})_{n\in\mathbb{N}}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with transitions given by

[TABLE]

if $v_{1}\in\kappa\mathbb{Z}^{d}$ and zero else. In particular, whatever the initial condition, $\tilde{X}_{n}\in\eta\mathbb{Z}^{d}$ and $\tilde{V}_{n}\in\kappa\mathbb{Z}^{d}$ for all $n\geqslant 1$. If $(\tilde{X}_{0},\tilde{V}_{0})=(\eta X_{0},\kappa V_{0})$, then $(\eta X_{n},\kappa V_{n})_{n\in\mathbb{N}}$ and $(\tilde{X}_{n},\tilde{V}_{n})_{n\in\mathbb{N}}$ have the same law. For this reason, in the rest of this section, only kinetic walks on $\mathbb{R}^{d}$ will be considered. See Section 4 for an example of kinetic walk on $\mathbb{Z}^{d}$.

This section is more concerned with a general and informal discussion than with rigorous results, the latter possibly requiring technical details that can be checked on explicit examples (see Section 4 in particular). In particular the results of this section (Proposition 4 and Theorem 5) are not new results.

3.1 First examples

Example 1. Let $U\in\mathcal{C}^{1}(\mathbb{R}^{d})$. Then the Hamiltonian dynamics $\partial_{t}(x_{t},v_{t})=(v_{t},-\nabla U(x_{t}))$ can be discretized as

[TABLE]

for some time-step $\delta>0$. This is a slight modification of the classical velocity Verlet integrator. It is a second-order scheme and, contrary to the basic Euler scheme, it is symplectic, like the Hamiltonian dynamics. From KAM theory and backward error analysis, it can be shown to conserve up to a high precision an approximate Hamiltonian, which ensures long-time stability, see [23] and in particular [23, Theorem 5.1] for long-time energy conservation.

Example 2. The Langevin diffusion (or sometimes underdamped Langevin diffusion)

[TABLE]

where $\gamma>0$ and $(B_{t})_{t\geqslant 0}$ is a standard Brownian motion on $\mathbb{R}^{d}$, can be approximated by similar second-order schemes (see [29, 7], references within and Section 5 for more details on this topic). For instance, the Ricci-Ciccotti scheme [44] reads

[TABLE]

where $(G_{n})_{n\in\mathbb{N}}$ is an i.i.d. sequence with law $\mathcal{N}(0,I_{d})$.

3.2 Sampling by thinning

The continous-time thinning and superposition method for sampling inhomogeneous Poisson processes, hence piecewise-deterministic Markov processes, is detailed e.g. in [32]. See Section 5.2 for an example of application. This section is concerned with its discrete analogous, which is essentially a rejection method applied to Bernoulli random variables (see also [34, 37] on similar topics).

Suppose that the transition $p$ can be decomposed as

[TABLE]

where, from a numerical point of view, computing $q$ and sampling according to $p_{1}$ is expensive, and sampling according to $p_{2}$ is not (for instance, $p_{2}(x,v;\cdot)=\delta_{v}(\cdot)$ for persistent chains). Suppose moreover that $q(x,v)\leqslant\tilde{q}(x,v)$ where $\tilde{q}$ is cheaper to compute thant $q$. Then, for $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$, a random variable $V\sim p(x,v;\cdot)$ can be sampled as follows. Draw two independent random variables $U_{1}$ and $U_{2}$ uniformly distributed over $[0,1]$. If $U_{1}\leqslant\tilde{q}(x,v)$ and $U_{2}\leqslant q(x,v)/\tilde{q}(x,v)$, draw $V$ according to $p_{1}(x,v:\cdot)$ else draw $V$ according to $p_{2}(x,v;\cdot)$. That way, obviously, $V\sim p(x,v;\cdot)$. The trick is that if $U_{1}>\tilde{q}(x,v)$ then we already know that $V$ has to be drawn according to $p_{2}(x,v;\cdot)$ and in that case we don’t even have to compute $q(x,v)$. The smaller is $\tilde{q}$, the higher is the computational gain.

We can go a bit further in two cases for which the first step $K$ where $q$ has to be computed, i.e. where $U_{1}<\tilde{q}(X_{K},V_{K})$, can be computed more efficiently than with Bernoulli variables at each step.

• •

If $\tilde{q}(x,v)=\tilde{q}\in(0,1)$ is constant. In that case, $K+1$ follows a geometric law with parameter $\tilde{q}$ and can be sampled through the inverse transformation method i.e. by setting $K=\lfloor\ln U/\ln\tilde{q}\rfloor$ with $U$ uniformly distributed over $[0,1]$ (which is particularly more efficient than with Bernoulli variables when $\tilde{q}$ is small).

• •

If $p_{2}(x,v;\cdot)=\delta_{f(x,v)}(\cdot)$ is deterministic, for some $f:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ (typically, for a persistent walk, $f(x,v)=v$). In that case, $K$ follows the distribution

[TABLE]

where $\varphi_{0}(x,v)=(x,v)$, $\varphi_{1}(x,v)=(x+\delta(v+f(x,v))/2,f(x,v))$ and $\varphi_{k+1}=\varphi_{k}\circ\varphi_{1}$ for all $k\in\mathbb{N}$. In particular cases, depending on $f$ and $\tilde{q}$, this distribution may again be sampled through the inverse transformation method.

In both cases, the algorithm is thus the following: draw $K$ as above and an independent $U_{2}$ uniformly distributed over $[0,1]$. Sample $(X_{n},V_{n})_{n\in\llbracket 0,K\rrbracket}$ as a kinetic chain associated to the transition $p_{2}$. If $U_{2}\leqslant q(X_{K},V_{K})/\tilde{q}(X_{K},V_{K})$, draw $V_{K+1}$ according to $p_{1}(X_{K},V_{K};\cdot)$, else according to $p_{2}(X_{K},V_{K};\cdot)$, and in both cases set $X_{K+1}=X_{K}+\delta(V_{K}+V_{K+1})/2$. Then, draw a new $K^{\prime}$ in a similar way as $K$, etc.

In the general case, of course each of the kernels $p_{1}$ and $p_{2}$ may also be decomposed in a similar way as $p$, so that we end up for some $N\geqslant 1$ with a decomposition $p=\sum_{i=1}^{N}p_{n}q_{n}$ where for each $n\in\llbracket 1,N\rrbracket$, $p_{n}$ is a transition kernel and $q_{n}\in[0,1]$ with some $\sum_{i=1}^{N}q_{n}=1$. Similarly, if $\tilde{q}(x,y)\leqslant\hat{q}(x,y)$ with $\hat{q}(x,v)$ cheaper to compute than $\tilde{q}(x,v)$, then we can sample a Bernoulli variable with parameter $q(x,v)$ as the product of three Bernoulli variables with respective parameters $\hat{q}(x,v)$, $\tilde{q}(x,v)/\hat{q}(x,v)$ and $q(x,v)/\tilde{q}(x,v)$. In other words, in the general decomposition we can decompose each weight $q_{n}$ as a product $\prod_{k=1}^{r_{n}}q_{n,k}$ for some $r_{n}\geqslant 1$ and $q_{n,k}\in[0,1]$ for each $k\in\llbracket 1,r_{n}\rrbracket$. At the end of the day we get a representation of the form

[TABLE]

that we can use to sample according to $p$ in such a way that the average cost of computation is minimized. See Sections 4 and 5 for examples and related questions, in particular the link with factorization for Metropolis acceptance probabilities in Section 4.3. See also [42] for an application in molecular dynamics.

3.3 Invariant measure

For MCMC applications, usual continuous-time kinetic Markov processes are designed to sample according to a given probability measure of the form $\mu(\text{d}x,\text{d}v)=\pi(\text{d}x)\otimes\nu(\text{d}v)$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$. The target is the position marginal $\pi$ and the velocity marginal $\nu$ can be chosen by the user, usual choices being Gaussian or uniform (over the sphere or a discrete set of velocities) distributions. By definition, $\mu$ is invariant for the kinetic walk on $\mathbb{R}^{d}$ associated to some kernel $p$ and time-step $\delta$ if

[TABLE]

for all $f\in\mathcal{M}_{b}(\mathbb{R}^{2d})$. Nevertheless, there are cases (such as Examples 1 and 2 of Section 3.1) where this condition is only approximately satisfied, typically $\mu$ is invariant only for the continuous-time process that is approximated by the discrete walk. In that case, for a fixed $\delta>0$, the kinetic walk, with timestep $\delta$ and associated to some transition $p_{\delta}$, typically admits some invariant measure $\mu_{\delta}$, which is also an invariant measure for the continuous-time Markov chain with generator

[TABLE]

In the typical case where $L_{\delta}$ converges, as $\delta$ vanishes, to some $L$ for which $\mu$ is invariant, it is possible to obtain explicit error bounds between $\mu_{\delta}$ and $\mu$. We recall here a general argument based on Stein’s method [35, Section 6.2].

Proposition 4**.**

Let $L$ and $L_{\delta}$ be two Markov generators on some Polish space $E$ and let $\mu,\mu_{\delta}\in\mathcal{P}(E)$ be invariant measures of, respectively, $L$ and $L_{\delta}$. Denote $(P_{t})_{t\geqslant 0}$ the semi-group associated to $L$. Suppose that there exist $C,\rho,h_{\delta}>0$ and two norms $N_{1}$ and $N_{2}$ on a subspace $\tilde{\mathcal{M}}$ of $M_{b}(E)$ such that the following holds: 1) $\tilde{\mathcal{M}}$ is dense in $(M_{b}(\mathbb{R}^{2d}),\|\cdot\|_{\infty})$ and is contained by the domains of $L$ and $L_{\delta}$; 2) for all $f\in\tilde{\mathcal{M}}$ and all $t\geqslant 0$, $N_{1}(P_{t}f-\mu f)\leqslant Ce^{-\rho t}N_{2}(f)/(1\vee\sqrt{t})$; 3) for all $f\in\tilde{\mathcal{M}}$, $\|(L-L_{\delta})f\|_{\infty}\leqslant h_{\delta}N_{1}(f)$. Then, denoting $N_{2}(\mu-\mu_{\delta})=\sup\{\mu(f)-\mu_{\delta}(f),\ N_{2}(f)\leqslant 1\}$, it holds:

[TABLE]

Proof.

Let $f\in\tilde{\mathcal{M}}$ and $g=\int_{0}^{\infty}P_{t}(f-\mu f)\text{d}t$, which is well defined and satisfies $N_{1}(g)\leqslant C(2/3+1/\rho)N_{2}(f)$ and solves the Poisson equation $Lg=\mu f-f$. Using that $\int L_{\delta}g\mu_{\delta}=0$ by invariance of $\mu_{\delta}$, we get that

[TABLE]

∎

Remark that, in Proposition 4, $C$ and $\rho$ only depend on the limit process $L$. For particular processes, their existence usually follows from regularization and ergodicity results for $L$ (with $N_{1}$ and $N_{2}$, typically, $L^{2}(\mu)$, $H^{1}(\mu)$ or $V$-norms associated to some Lyapunov function), see e.g. Section 4.4.4, [18, 35] or, for velocity jump processes, [40, Section 3]. Then, if we are given a family of generators $L_{\delta}$ for all $\delta>0$ such that $\|(L-L_{\delta})f\|_{\infty}\leqslant h_{\delta}N_{1}(f)$ with $h_{\delta}\rightarrow 0$ for $\delta\rightarrow 0$, we get a quantitative estimate on the convergence $\mu_{\delta}\rightarrow\mu$ (remark that, though the uniqueness of the invariant measure is ensured for $L$ from the geometric ergodicity assumption, we haven’t assumed the uniqueness of the invariant measure for $L_{\delta}$).

In fact, reiterating this argument, following the Talay-Tubaro method [46], an expansion of the bias $\mu_{\delta}f-\mu f$ in term of powers of $\delta$ can be computed and a Romberg extrapolation (or related methods) can be used to kill the first order terms, see [46, Section 2.3] and [31, Section 3.3.4].

Proposition 4 can be shown to apply in Example 2 of Section 3.1, i.e. the Langevin dynamics, see for instance [30]. On the contrary, it does not apply to the Hamiltonian dynamics of Example 1. Indeed, in that case, the limit process admits many invariant measures (because of energy conservation) and thus $C$ and $\rho$ cannot exist.

With Proposition 4, we have seen that the convergence of the invariant measures is related to the convergence of the generators. Now, the latter is classically related to the convergence of the processes and we address this question in the next section. See Theorem 5 below for some considerations on the convergence of the generators in the specific case of kinetic walks.

3.4 Scaling limits

In this section, we consider for all $\varepsilon\in(0,1]$ a kinetic walk $(X_{n}^{\varepsilon},V_{n}^{\varepsilon})_{n\in\mathbb{N}}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with timestep $1$, kernel $p_{\varepsilon}$ and initial condition $(x_{0}^{\varepsilon},v_{0}^{\varepsilon})$. We are interested in the possible convergence of this chain, possibly rescaled as $\varepsilon$ vanishes, toward a continuous-time process. The regime for which $(X_{n}^{\varepsilon})_{n\in\mathbb{N}}$ converges toward an elliptic diffusion has been abundantly studied for second-order chains, see [20, Section 5] and references therein, and for this reason we will mostly focus on the cases where the full system $(X_{n}^{\varepsilon},V_{n}^{\varepsilon})_{n\in\mathbb{N}}$ converges toward a continuous-time kinetic process $(Y_{t},W_{t})_{t\geqslant 0}$, where kinetic means that $Y_{t}=Y_{0}+\int_{0}^{t}W_{s}\text{d}s$.

To alleviate notations, unless otherwise specified, we drop the $\varepsilon$ superscript in all the rest of the section and simply write $(X_{n},V_{n})$. We start with an informal discussion on the scaling in the simple case where the dynamics are homogeneous with respect to the space variable $x$. Since, in order to expect a continuous-time limit, $X_{n}$ should be nearly constant over a large number of steps as $\varepsilon$ goes to zero, this homogeneous case should be expected to describe the short time dynamics of the general case.

3.4.1 The space homogeneous case

Consider the case where $p_{\varepsilon}(x,v;\cdot)=h_{\varepsilon}(v;\cdot)$ for some transition kernel $h_{\varepsilon}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. In that case, $(V_{n})_{n\in\mathbb{N}}$ is a Markov chain by itself. Since

[TABLE]

then, provided that, say, the variance of the last term is bounded uniformly in $n$, the situation is essentially the same as the case of correlated random walks. In term of time/space scaling, different cases may be distinguished concerning the limit of $(X_{n})_{n\in\mathbb{N}}$:

• •

If $(V_{k})_{k\in\mathbb{N}}$ is in fact an i.i.d. sequence – namely if $v\mapsto h_{\varepsilon}(v,\cdot)$ is constant – then, up to a vanishing term, $(X_{n})_{n\in\mathbb{N}}$ is a simple random walk. If it admits a continuous-time limit, then the latter is necessarily a Levy process, and conversely any Levy process $(L_{t})_{t\geqslant 0}$ may be obtained as a scaling limit of such a walk, even if we restrict the question to walks on $\mathbb{Z}^{d}$: indeed, considering $V_{k}=\lfloor(L_{(k+1)\varepsilon}-L_{k\varepsilon})/\varepsilon\rfloor$, then $(\varepsilon X_{\lfloor t/\varepsilon\rfloor})_{t\geqslant 0}\rightarrow(L_{t})_{t\geqslant 0}$.

Let us consider two particular cases. First, suppose that there exist $\eta:(0,1]\rightarrow(0,+\infty)$ such that

[TABLE]

(with possibly $\mu=0$ or $\Sigma=0$). Then, from Donsker’s Theorem, provided that $\eta(\varepsilon)x_{0}^{\varepsilon}$ converges to some $x_{0}^{*}\in\mathbb{R}^{d}$, we get the drifted Brownian motion

[TABLE]

where $(B_{t})_{t\geqslant 0}$ is a standard Brownian motion. Second, suppose that there exist $\lambda>0$, $\mu\in\mathbb{R}^{d}$, $\nu\in\mathcal{P}(\mathbb{R}^{d})$, $\eta:(0,1]\rightarrow(0,+\infty)$ and $\mu:(0,1]\rightarrow\mathbb{R}^{d}$ such that

[TABLE]

Then the cardinality $N_{t}^{\varepsilon}$ of $\{n\in\mathbb{N}\ :\varepsilon n\leqslant t,\ V_{k}\neq\mu(\varepsilon)\}$ converges as $\varepsilon$ vanishes to a Poisson process $N_{t}$ with intensity $\lambda$. As a consequence, considering an i.i.d. sequence $(W_{k})_{k\in\mathbb{N}}$ with law $\nu$ independent from $(N_{t}^{\varepsilon})_{t\geqslant 0}$ for all $\varepsilon$ and from $(N_{t})_{t\geqslant 0}$, provided that $\eta(\varepsilon)x_{0}^{\varepsilon}$ converges to some $x_{0}^{*}\in\mathbb{R}^{d}$, we get the drifted compound Poisson process

[TABLE]

In the more general situation where $(V_{k})_{k\in\mathbb{N}}$ is not necessarily an i.i.d. sequence, if there exists a $\kappa:(0,1]\rightarrow(0,+\infty)$ such that $\left(\kappa(\varepsilon)V_{\lfloor t/\varepsilon\rfloor}\right)_{t\geqslant 0}$ converges in distribution toward a continuous-time Markov process $(W_{t})_{t\geqslant 0}$ then, denoting $\eta(\varepsilon)=\varepsilon\kappa(\varepsilon)$, (6) reads

[TABLE]

Integration being continuous with respect to Skorohod convergence, provided that $\eta(\varepsilon)x_{0}^{\varepsilon}\rightarrow x_{0}^{*}$ as $\varepsilon$ vanishes, this yields

[TABLE]

Remark that, of course, the scaling factor for the space variable is fixed by the scaling factors of the time and velocity variables.

• •

For instance, if $(V_{k})_{k\in\mathbb{N}}$ is a random walk on $\mathbb{Z}^{d}$, then as seen before it can converge toward a drifted Brownian motion, in which case the scaling limit of $(X_{n},V_{n})_{n\in\mathbb{N}}$ is the Langevin diffusion, i.e. the solution of the SDE

[TABLE]

where $(B_{t})_{t\geqslant 0}$ is a standard Brownian motion on $\mathbb{R}^{d}$.

• •

Alternatively, if there exist $\lambda>0$, $\nu\in\mathcal{P}(\mathbb{R}^{d})$ and $\kappa:(0,1]\rightarrow(0,+\infty)$ such that

[TABLE]

then $(\kappa(\varepsilon)V_{\lfloor t/\varepsilon\rfloor})_{t\geqslant 0}$ converges as $\varepsilon$ vanishes toward $(Y_{N_{t}})_{t\geqslant 0}$, where $(Y_{k})_{k\in\mathbb{N}}$ is an i.i.d. sequence with law $\nu$ and $(N_{t})_{t\geqslant 0}$ is a Poisson process with intensity $\lambda$, independent from $(Y_{k})_{k\in\mathbb{N}}$. In that case, $(X_{t},V_{t})_{t\geqslant 0}$ is the velocity jump process associated to the linear Boltzmann (or BGK) equation [9].

These different examples highlighted three cases: if there is no inertia, the velocity tends to mix fast and $(X_{n})_{n\in\mathbb{N}}$ is Markovian. If there is some inertia in the sense that the velocity tends to stay aligned from one step to the other but with possible small fluctuations, the limit is a kinetic diffusion. If the velocity is rigorously constant for large times, the chain converges toward a velocity jump process. Of course this is a non-exhaustive list: in general, drift, diffusion, Poisson or $\alpha$-stable jumps can all be present in the limit, either kinetic or not (see in particular Section 5). But with these three regimes we cover the cases of the processes classically used in MCMC sampling.

3.4.2 The case of kinetic limits

We keep the notations of the beginning of Section 3.4. In particular, for all $\varepsilon\in(0,1]$, $(X_{n}^{\varepsilon},V_{n}^{\varepsilon})_{n\in\mathbb{N}}$ is a kinetic walk on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ with transition $p_{\varepsilon}$ and time-step $1$. Let $\kappa:(0,1]\rightarrow(0,+\infty)$ and, for all $\varepsilon\in(0,1]$, set $\eta(\varepsilon)=\varepsilon\kappa(\varepsilon)$. For any fixed $x\in\mathbb{R}^{d}$, let $(\tilde{V}_{k}^{\varepsilon,x})_{k\in\mathbb{N}}$ be the Markov chain on $\mathbb{R}^{d}$ with transitions

[TABLE]

and let $N_{t}^{\varepsilon}$ be a Poisson process with intensity $1/\varepsilon$. Denote $\tilde{L}_{\varepsilon}$ the operator defined on $\mathcal{M}_{b}(\mathbb{R}^{2d})$ by

[TABLE]

which is the infinitesimal generator of the Feller process $(\tilde{X}_{t}^{x},\tilde{V}_{N_{t}^{\varepsilon}}^{\varepsilon,x})_{t\geqslant 0}$, where we simply set $\tilde{X}_{t}^{x}=x$ for all $t\geqslant 0$ (hence the link with the space homogeneous case). A direct corollary of [27, Theorem 17.28] is the following:

Theorem 5**.**

Suppose that there exists a Feller generator $\tilde{L}$ on $\mathbb{R}^{2d}$ with domain containing $\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$ and such that $\|\tilde{L}_{\varepsilon}f-\tilde{L}f\|_{\infty}\rightarrow 0$ as $\varepsilon\rightarrow 0$ for all $f\in\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$. Define the operator $L$ by

[TABLE]

Suppose that $L$ is the infinitesimal generator of a Feller process $(Y_{t},W_{t})_{t\geqslant 0}$ and that $\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$ is a core of $L$. Suppose that for all compact set $\mathcal{K}\subset\mathbb{R}^{d}$,

[TABLE]

Finally, suppose that $(\eta(\varepsilon)X_{0}^{\varepsilon},\kappa(\varepsilon)V_{0}^{\varepsilon})$ converges in law toward $(Y_{0},W_{0})$ as $\varepsilon$ vanishes. Then

[TABLE]

See some applications with Proposition 7, Section 5.2 or [42].

Proof.

Denoting $(Y_{t}^{\varepsilon},W_{t}^{\varepsilon})=(\eta(\varepsilon)X_{N_{t}^{\varepsilon}}^{\varepsilon},\kappa(\varepsilon)V_{N_{t}^{\varepsilon}}^{\varepsilon})$, the generator $L_{\varepsilon}$ of $(Y_{t}^{\varepsilon},W_{t}^{\varepsilon})_{t\geqslant 0}$ is defined on $\mathcal{M}_{b}(\mathbb{R}^{2d})$ by

[TABLE]

Note that $L_{\varepsilon}f\in\mathcal{M}_{b}(R^{2d})$ for all $f\in\mathcal{M}_{b}(R^{2d})$, so that if $(X_{0}^{\varepsilon},V_{0}^{\varepsilon})=(x,v)$,

[TABLE]

with a negligible term uniform in $(x,v)$ (which will be the case of all the negligible terms in the rest of the proof). This means that $\mathcal{M}_{b}(\mathbb{R}^{2d})$, hence $\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$, is included in the strong domain of $L_{\varepsilon}$ for all $\varepsilon\in(0,1]$. From [27, Theorem 17.28] and the assumption that $\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$ is a core for $L$, it only remains to check that $\|L_{\varepsilon}f-Lf\|_{\infty}$ vanishes with $\varepsilon$ for all $f\in\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$. Now, indeed, for $f\in\mathcal{C}^{2}_{c}(\mathbb{R}^{2d})$,

[TABLE]

Considering a ball $\mathcal{B}\subset\mathbb{R}^{d}$ of some radius $R$ such that the support of $f$ is included in $\mathcal{B}^{2}$, we bound

[TABLE]

As a consequence, for some $C_{f}>0$, for all $(x,v)\in\mathbb{R}^{2d}$,

[TABLE]

Condition (7) concludes. ∎

4 The discrete Zig-Zag walk

This section is devoted to the definition and study of a discrete-space analogous of the Zig-Zag process on $\mathbb{R}^{d}$.

4.1 Definition

Let $d\in\mathbb{N}_{*}$, $U:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ be such that $\mathcal{Z}=\sum_{x\in\mathbb{Z}^{d}}\exp(-U(x))<+\infty$, $\pi(x)=\exp(-U(x))/\mathcal{Z}$ be the associated Gibbs distribution and $\mu(x,v)=\pi(x)/{2^{d}}$ for $v\in\{-1,1\}^{d}$ and $x\in\mathbb{Z}^{d}$. For $i\in\llbracket 1,d\rrbracket$, denote $e_{i}$ the $i^{th}$ vector of the canonical basis of $\mathbb{R}^{d}$ and let

[TABLE]

That way, $p:\mathbb{Z}^{d}\times\{-1,1\}^{d}\rightarrow\mathcal{P}(\{-1,1\}^{d})$. We call Zig-Zag walk on $\mathbb{Z}^{d}$ the kinetic walk $(X_{n},V_{n})_{n\in\mathbb{N}}$ associated to this kernel $p$ with timestep $\delta=1$, i.e. the Markov chain on $\mathbb{Z}^{d}\times\{-1,1\}^{d}$ whose transitions are given by (5). Remark that, for $d=1$, we retrieve the chain studied in Section 2.

A random variable $V\sim p(x,v;\cdot)$ in $\{-1,1\}^{d}$ can be sampled as follows. Set $Y_{0}=x$, and suppose by induction that $Y_{k-1}\in\mathbb{Z}^{d}$ has been defined for some $k\in\llbracket 1,d\rrbracket$. Set $V_{k}=v_{k}$ with probability $q_{k}(Y_{k-1},v_{k})$ and $V_{k}=-v_{k}$ else, and in either case set $Y_{k}=Y_{k-1}+(v_{k}+V_{k})/2$. Then $V$ is distributed according to $p(x,v;\cdot)$, and $X:=Y_{d}=x+(v+V)/2$. In other words, this is a Gibbs algorithm based on the Zig-Zag walk on $\mathbb{Z}$: one step of the Zig-Zag walk in $\mathbb{Z}^{d}$ is the result of $d$ successive one-dimensional Zig-Zag steps on each coordinate, the others being fixed.

If we want the coordinates to play a symmetric role in the transition, for $\sigma$ a permutation of $\llbracket 1,d\rrbracket$ we can define $p_{\sigma}(x,v;\cdot)$ to be the law of $W_{\sigma^{-1}}$ when $W\sim p(x_{\sigma},v_{\sigma};\cdot)$, where $u_{\sigma}$ for $u\in\mathbb{Z}^{d}$ and $\sigma\in\mathfrak{S}_{d}$ denotes $(u_{\sigma(1)},\dots,u_{\sigma(d)})$. This accounts to use the order given by $\sigma$ to update the coordinates. Then

[TABLE]

corresponds to a transition where the order is sampled at random at each step of the Zig-Zag walk. There is no particular practical interest to consider $p_{sym}$ rather than $p$, and moreover any result on $p$ can straightforwardly be adapted to $p_{\sigma}$ by renumbering of the coordinates, and then to $p_{sym}$.

4.2 Equilibrium and scaling limit

Proposition 6**.**

The probability distribution $\mu$ is invariant for the Zig-Zag walk on $\mathbb{Z}^{d}$.

Proof.

As proven in Section 2, for all fixed $(x_{j},v_{j})_{j\in\llbracket 2,d\rrbracket}$, the transition on $\mathbb{Z}\times\{-1,1\}$ defined by

[TABLE]

admits the conditional law $(y,w)\mapsto\pi(y,x_{2},\dots,x_{d})/2$ as an invariant measure. As a consequence, the transition of the Markov chain $(\tilde{X}_{n},\tilde{V}_{n})$ on $\mathbb{Z}^{d}\times\{-1,1\}^{d}$ with $(\tilde{X}_{n,1},\tilde{V}_{n,1})=(Y_{n},W_{n})$ and $(\tilde{X}_{n,j},\tilde{V}_{n,j})=(\tilde{X}_{0,j},\tilde{V}_{0,j})$ for $j\neq 1$ also fixes $\mu$. Since the transition of the Zig-Zag walk is the composition of $d$ such transitions, it fixes $\mu$. ∎

Remark that, if the target law is of a tensor form $\pi(x)=\prod_{i=1}^{d}\pi_{i}(x_{i})$, then the coordinates of a Zig-Zag walk are just $d$ independent one-dimensional Zig-Zag walks (which is similar to the continuous-time process).

Recall that the continuous-time Zig-Zag process on $\mathbb{R}^{d}$ associated to a potential $H$ is the Markov process on $\mathbb{R}^{d}\times\{-1,1\}^{d}$ with generator

[TABLE]

where we denote by $v_{-i}$ the vector of $\{-1,1\}^{d}$ obtained from $v$ by multiplying its $i^{th}$ coordinate by $-1$. The following is the extension of Theorem 3 in larger dimension.

Theorem 7**.**

For $H\in\mathcal{C}^{2}(\mathbb{R}^{d})$ that goes to infinity at infinity, for all $\varepsilon>0$, define $U_{\varepsilon}:\mathbb{Z}^{d}\mapsto\mathbb{R}$ by $U_{\varepsilon}(x)=H(\varepsilon x)$ for all $x\in\mathbb{Z}^{d}$. Let $(X^{\varepsilon}_{k},V^{\varepsilon}_{k})_{k\in\mathbb{N}}$ be the Zig-Zag walk on $\mathbb{Z}^{d}$ associated to $U_{\varepsilon}$ and with some initial condition $(x_{0}^{\varepsilon},v_{0})$. Suppose that $\varepsilon x_{0}^{\varepsilon}$ converges to some $x_{0}^{*}\in\mathbb{R}^{d}$ as $\varepsilon$ vanishes. Then

[TABLE]

where $(Y_{t},W_{t})_{t\geqslant 0}$ is a Zig-Zag process on $\mathbb{R}^{d}$ associated to $H$ and with $(Y_{0},W_{0})=(x_{0}^{*},v_{0})$.

Proof.

Let us show that Theorem 5, or rather directly [27, Theorem 17.28], applies. First, following [14], we can see that the continuous-time process can be smoothly and compactly approximated (in the sense of [15, Definition 20]) by replacing its continuous jump rates by $\mathcal{C}^{\infty}$ jump rates (see the case of the BPS in [15, Proposition 23] for details). From [15, Theorem 21], this proves that $\mathcal{C}^{1}_{c}(\mathbb{R}^{d}\times\{-1,1\}^{d})$ is a core for the strong generator $L$ of the Zig-Zag process. Denote $L_{\varepsilon}$ the generator of $(\varepsilon X^{\varepsilon}_{N_{t}},V^{\varepsilon}_{N_{t}})_{k\in\mathbb{N}}$ where $(N_{t})_{t\geqslant 0}$ is a Poisson process with intensity $1/\varepsilon$. Then all bounded measurable functions $f$ are in the domain of $L_{\varepsilon}$ and

[TABLE]

where $p_{\varepsilon}$ is given by (4.1) with $U=U_{\varepsilon}$. Next, for all $i\in\llbracket 1,d\rrbracket$,

[TABLE]

where the negligible term is uniform over all compact set of $\mathbb{R}^{d}\times\{-1,1\}^{d}$ since $H$ is $\mathcal{C}^{2}$ (and this will be the case for all the negligible terms below). More generally,

[TABLE]

from which

[TABLE]

On the other hand, if $f\in\mathcal{C}^{2}_{c}(\mathbb{R}^{d}\times\{-1,1\}^{d})$,

[TABLE]

so that

[TABLE]

and Theorem [27, Theorem 17.28] concludes. ∎

Note that the space/time scaling in Theorem 7 is ballistic. It means in particular that, in $n=\lfloor 1/\varepsilon\rfloor$ steps, the Zig-Zag walk with potential $U_{\varepsilon}$ is at distance of order $n$ (and not $\sqrt{n}$ as in the diffusive case) from its starting point.

4.3 Thinning and factorization

In order to sample the Zig-Zag walk, a priori, at each time step, $U(x)$ has to be computed for $d+1$ values of $x$. However, thanks to the thinning method recalled in Section 3.2, this computational cost may drop if simple bounds are known on the increments of $U$.

Moreover, the factorization principle used for the continuous-time Zig-Zag process in [4] to do subsampling is still available here. Suppose that we can decompose $U(x+v)-U(x)=\sum_{j=1}^{M}f_{j}(x,v)$ for all $(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}$ for some $M\in\mathbb{N}_{*}$ and $f_{j}:\mathbb{Z}^{d}\times\{-1,1\}^{d}\rightarrow\mathbb{R}$. This is for instance the case if $U=\sum_{j=1}^{M}U_{j}$, in which case we can take $f_{j}(x,v)=U_{j}(x+v)-U_{j}(x)$, but in general the $f_{j}$’s are not required to be discrete gradients. Consider the Zig-Zag walk as defined above except that the probability $q_{i}(x,v_{i})$ is replaced by

[TABLE]

Proposition 8**.**

This Zig-Zag walk with $\tilde{q}_{i}$ still admits $\mu$ as an invariant measure.

Proof.

As in the proof of Proposition 6, we just have to prove the result for $d=1$. Following Section 2, this stems from the same result applied to classical Metropolis-Hastings algorithms. Indeed, let $q$ be a symmetric Markov kernel on a space $E$ and let $\alpha:E\times E\rightarrow[0,1]$ be the acceptance probability of a Metropolis-Hastings chain with proposal $q$, namely a chain with transition kernel $p(x,y)=q(x,y)\alpha(x,y)$ for $y\neq x$. Then this chain is reversible with respect to a probability $\mu$ if and only if

[TABLE]

In particular, if $\mu(x)=C\prod_{i=1}^{M}g_{i}(x)$ for some positive $g_{i}$’s and a normalization constant $C$ then, setting $\alpha_{i}(x,y)=\min\left(1,g_{i}(x)/g_{i}(y)\right)$ for all $i\in\llbracket 1,M\rrbracket$ and $x,y\in E$ ensures that

[TABLE]

Taking the product over $i\in\llbracket 1,M\rrbracket$ and noting that $\min(1,a/b)=\exp\left(-(\ln(a/b))_{+}\right)$ for $a,b>0$ concludes. ∎

The bad side of factorization is that it increases the number of rejections (for the Metropolis-Hastings algorithm, hence of collisions for the Zig-Zag process). On the other hand, if $U(x+v)-U(x)$ can be decomposed in a part that is cheap to compute and a part that may be expensive to compute but is small and for which an efficient bound is available for thinning, then it may give a significant computational gain. This is particularly well-adapted for multi-scale potentials, as we can see in Section 5 on a similar problem.

4.4 Irreducibility, Ergodicity, CLT

4.4.1 Irreducibility

Irreducibility is a delicate question for the continuous Zig-Zag process, see [6]. Here, for the discrete Zig-Zag walk, we will only tackle the restrictive case where, following the definitions of [6], all velocities are asymptotically flippable111The nice proof of irreducibility of [6] under weaker conditions on $U$ may possibly be partially adapted for the discrete Zig-Zag walk. Nevertheless, note that the smoothness condition on $U$ has no discrete counterpart. In any case, the potential $U(x)=\|x\|_{\infty}$ would still be a counter-example for which the conclusion of Proposition 9 would not hold. and thus the proof is similar to the Gaussian case with dominant diagonal of [6, Corollary 1]. Anyway we are interested in the exponentially fast convergence toward equilibrium under the assumption of Proposition 10 below, which is even stronger. Moreover, with this restriction, we can focus on the specificities of the discrete realm.

Indeed, let us call $\sigma(x,v):=\left((-1)^{x_{i}}v_{i}\right)_{i\in\llbracket 1,d\rrbracket}\in\{-1,1\}^{d}$ the signature of $(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}$. If $(X_{n},V_{n})_{n\in\mathbb{N}}$ is a Zig-Zag walk on $\mathbb{Z}^{d}$ then, like in the one-dimensional case, $\sigma(X_{n},V_{n})=-\sigma(X_{n-1},V_{n-1})$. In particular, denoting $\mathcal{A}_{s}=\{(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}\ :\ \sigma(x,v)=s\}$ then $\mathcal{A}_{s}\cup\mathcal{A}_{-s}$ is fixed by the Zig-Zag walk (see Figure 2). Therefore, in dimension larger than 1, the Zig-Zag walk is not irreducible on $\mathbb{Z}^{d}\times\{-1,1\}^{d}$.

Proposition 9**.**

Suppose that there exist $R>0$ such that $U(x+v_{i}e_{i})>U(x)$ for all $i\in\llbracket 1,d\rrbracket$, $x\in\mathbb{Z}^{d}$, $v\in\{-1,1\}^{d}$ with $x_{i}v_{i}>R$. Then for all $s\in\{-1,1\}^{d}$ the Zig-Zag walk on $\mathbb{Z}^{d}$ associated to $U$ is irreducible on $\mathcal{A}_{s}\cup\mathcal{A}_{-s}$.

Proof.

Let $s\in\{-1,1\}^{d}$ be fixed, and let $(x,v)\in\mathcal{A}_{s}\cup\mathcal{A}_{-s}$. Remark that $(x,w)\in\mathcal{A}_{s}\cup\mathcal{A}_{-s}$ if and only if $w\in\{-v,v\}$. We say that we can reach $(y,w)$ from $(x,v)$ if there is a path from $(x,v)$ to $(y,w)$ that has a non-negative probability for the Zig-Zag walk. Starting from $(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}$, we can reach all the points $(x+nv,v)$ with $n\in\mathbb{N}$. For $n$ large enough, $v_{i}(x_{i}+nv_{i})>R$ for all $i\in\llbracket 1,d\rrbracket$ so that each coordinate has a non-negative probability to flip its velocity in the next step. As a consequence, from $(x+nv,v)$ with such a $n$, $(x+nv+(v+w)/2,w)$ can be reached for all $w\in\{-1,1\}^{d}$. In particular, if $w=-v$, since $(x,-v)$ can be reached from $(x+nv,-v)$, we see by transitivity that $(x,-v)$ can be reached from $(x,v)$.

Second, let us show that for all $j\in\llbracket 1,d\rrbracket$ and all $a\in\{-1,1\}$, $(x+ae_{j},-v+2v_{j}e_{j})$ can be reached from $(x,v)$. Remark that this will conclude the proof: indeed, repeating this, then for all $x^{\prime}\in\mathbb{Z}^{d}$, there will exist $w^{\prime}\in\{-1,1\}^{d}$ such that $(x^{\prime},w^{\prime})$ (and then $(x^{\prime},-w^{\prime})$ by the previous result) can be reached from $(x,v)$. Since $\mathcal{A}_{s}\cup\mathcal{A}_{-s}$ is fixed by the Zig-Zag walk, if $(x^{\prime},z)\in\mathcal{A}_{s}\cup\mathcal{A}_{-s}$ then necessarily $z\in\{w^{\prime},-w^{\prime}\}$, and thus we will have obtained that all points of $\mathcal{A}_{s}\cup\mathcal{A}_{-s}$ will be reachable from $(x,v)$.

Hence, fix $j\in\llbracket 1,d\rrbracket$, $a\in\{-1,1\}$ and set $(x^{\prime},v^{\prime})=(x+ae_{j},-v+2v_{j}e_{j})$. Since $(x,-v)$ can be reached from $(x,v)$ we can suppose that $v_{j}=-a$. Let $n_{1}$ and $n_{2}$ be large enough so that, for all $i\in\llbracket 1,d\rrbracket$,

[TABLE]

Consider the following path: from $(x,v)$, go to $(x+n_{1}v,v)$, flip the $j^{th}$ velocity, which gives $(x+(n_{1}+1)v-v_{j}e_{j},-v^{\prime})$, go straight to $(x+(n_{1}+1)v-v_{j}e_{j}-n_{2}v^{\prime},-v^{\prime})$, flip all the velocities but the $j^{th}$, which gives $(x+(n_{1}+1)v-2v_{j}e_{j}-n_{2}v^{\prime},-v)$, go straight to $(x+v-2v_{j}e_{j}-n_{2}v^{\prime},-v)$, flip the $j^{th}$ velocity, which gives $(x-v_{j}e_{j}-n_{2}v^{\prime},v^{\prime})$ and go straight to $(x^{\prime},v^{\prime})$ (see Figure 3). It is clear that, in view of the conditions on $n_{1}$ and $n_{2}$ the two first flips have a non-negative probability. For the third one, remark that the $j^{th}$ coordinate of $x+v-2v_{j}e_{j}-n_{2}v^{\prime}$ is $x_{j}^{\prime}-n_{2}v_{j}^{\prime}$ and that $v_{j}^{\prime}=v_{j}$. Hence, the condition that $-v_{i}^{\prime}(x_{i}^{\prime}-n_{2}v_{i}^{\prime})>R$ ensures that the third flip, hence the whole path, has a non-negative probability, which concludes.

∎

4.4.2 A Lyapunov function

For some fixed $a,b>0$ and for $x\in\mathbb{Z}^{d}$, $v\in\{-1,1\}^{d}$ and $i\in\llbracket 1,d\rrbracket$, denote

[TABLE]

and $\mathcal{V}(x,v)=\sum_{i=1}^{d}\mathcal{V}_{i}(x,v)$. Consider $(X_{n},V_{n})_{n\in\mathbb{N}}$ the Zig-Zag walk on $\mathbb{Z}^{d}$ associated to $U:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ and denote $Q$ its transition operator, namely

[TABLE]

Proposition 10**.**

Suppose that there exist $R,h>0$ such that for all $i\in\llbracket 1,d\rrbracket$, $x\in\mathbb{Z}^{d}$ and $v\in\{-1,1\}^{d}$,

[TABLE]

Then, for all choice of $a,b>0$ and for all $x\in\mathbb{Z}^{d}$ and $v\in\{-1,1\}^{d}$,

[TABLE]

Proof.

For all $x\in\mathbb{Z}^{d}$, $v\in\{-1,1\}^{d}$,

[TABLE]

If $x_{1}v_{1}>R$, $q_{1}(x,v_{1})\leqslant e^{-h}$, so that

[TABLE]

If $x_{1}v_{1}<-R$ then $(x_{1}+v_{1})v_{1}\leqslant 0$ and $q_{1}(x,v_{1})=1$, and thus

[TABLE]

If $|x_{1}v_{1}|\leqslant R$ then

[TABLE]

Since $Q$ is the result of $d$ consecutive and identical one-dimensional transitions, we get the result by summing over $i\in\llbracket 1,d\rrbracket$. ∎

Taking $a=h/2$ and $b$ arbitrarily large we get that $Q\mathcal{V}\leqslant\gamma\mathcal{V}+C$ with $\gamma<1$ (in fact $\gamma$ arbitrarily close to $e^{-h/2}$), which means that $\mathcal{V}$ is a Lyapunov function for $Q$.

Remark that similar computations in the case of the continuous-time Zig-Zag process on $\mathbb{R}^{d}$ show that if $(v_{i}\partial_{x_{i}}U(x))_{+}\geqslant h\mathbbm{1}_{\{x_{i}v_{i}>0\}}$ for all $x\in\mathbb{R}^{d}$ such that $|x_{i}|>R$ then

[TABLE]

where $\varphi(s)$ is some smooth approximation of $sign(s)$, is a Lyapunov function for the continuous-time process. Note that this condition on $U$ is not covered by [6, Condition 3] since the latter constrains $|\nabla U(x)|$ to go to infinity at infinity, excluding Laplace-tail distributions. Hence, our computations extends the scope of Lemma 2 (hence Theorem 2) of [6]. Note that condition (9) holds when $U(x)=\sum_{i=1}^{d}|x_{i}|$ but not when $U(x)=|x|$. Of course, condition (9) roughly means that the different coordinates are more or less independent at infinity and thus it is not surprising that we recover, in a discrete-space case, the one-dimensional computations of [17, Proposition 2.8].

In the following, under condition (9) we fix $a=h/2$ and $b$ such that $e^{-b}=e^{-h/4}-e^{-h/2}$, in which case (10) implies

[TABLE]

This classically yields uniform in time exponential moment bounds on the Zig-Zag walk, since for all $n\in\mathbb{N}$,

[TABLE]

and thus, if $(X_{0},V_{0})=(x,v)$, we bound

[TABLE]

4.4.3 Long-time convergence

For $x\in\mathbb{Z}^{d}$ and $s\in\{-1,1\}^{d}$, let $v(x,s)$ be the (unique) vector of $\{-1,1\}^{d}$ such that $\sigma(x,v(x,s))=s$. Then $\mathcal{A}_{s}=\{(x,v(x,s))\,,x\in\mathbb{Z}^{d}\}$ and $\mathcal{A}_{-s}=\{(x,-v(x,s))\,,x\in\mathbb{Z}^{d}\}$. Consider on $\mathcal{A}_{s}$ the probability measure

[TABLE]

For a given $\mathcal{W}:\mathbb{Z}^{d}\times\{-1,1\}^{d}\rightarrow[1,+\infty)$, we endow $\mathcal{P}_{\mathcal{W}}:=\{\mu\in\mathcal{P}(\mathbb{Z}\times\{-1,1\}^{d})\ :\mu(\mathcal{W})<+\infty\}$ with the norm

[TABLE]

which makes it complete.

Theorem 11**.**

Suppose that $U$ admits a strict local minimum at some $x_{*}\in\mathbb{Z}^{d}$ and that there exist $R,h>0$ such that (9) holds for all $i\in\llbracket 1,d\rrbracket$, $x\in\mathbb{Z}^{d}$ and $v\in\{-1,1\}^{d}$. Set $\mathcal{W}(x,v)=\sum_{i=1}^{d}e^{h|x_{i}|/2}$. Then, there exist $C>0$ and $\rho\in(0,1)$ such that for all $x\in\mathbb{Z}^{d}$, $v\in\{-1,1\}^{d}$ and $n\in\mathbb{N}_{*}$,

[TABLE]

Proof.

Let $(x,v),(x^{\prime},v^{\prime})\in\mathcal{A}_{s}$ for some $s\in\{-1,1\}^{d}$. Proposition 9 gives a path from $(x,v)$ to $(x^{\prime},v^{\prime})$ whose transitions are non-negative under $Q$. Since $Q(\mathcal{A}_{s})\subset Q(\mathcal{A}_{-s})$, the length of such a path is necessarily even, from which $Q^{2}$ is irreducible on $\mathcal{A}_{s}$. The path $(x_{*},v(x_{*},s))\rightarrow(x_{*},-v(x_{*},s))\rightarrow(x_{*},v(x_{*},s))$ having a non-negative probability under $Q$, $Q^{2}$ is aperiodic on $\mathcal{A}_{s}$. As a consequence, for all $(x,v),(x^{\prime},v^{\prime})\in\mathcal{A}_{s}$, there exist $n_{0}$ such that $Q^{2n}\left((x,v),(x^{\prime},v^{\prime})\right)>0$ for all $n\geqslant n_{0}$. From Proposition 6, $\mu_{s}$ is invariant for $Q^{2}$.

Consider $\mathcal{V}$ as defined in Section 4.4.2 with $a=h/2$ and $b$ large enough, so that

[TABLE]

for some $\gamma\in(0,1)$, $C>0$. The set $\{\mathcal{V}\leqslant 4C\}$ is finite. Fix any point $x^{\prime}\in\mathbb{Z}^{d}$ (say, $x^{\prime}=x_{*}$) and let

[TABLE]

Then

[TABLE]

From [24, Theorem 1.2] applied to $Q^{2m}$, the Foster-Lyapunov condition (11) (applied with $n=m$) and the Doeblin condition (12) imply the existence of $\rho\in(0,1)$ and $C^{\prime}>0$ such that for all $x\in\mathbb{Z}^{d}$, $v\in\{-1,1\}^{d}$ and $k\in\mathbb{N}_{*}$,

[TABLE]

In fact, even without the Doeblin condition, following the proof of [24, Theorem 1.2], we also get that the Lyapunov condition given by Proposition 10 alone implies the following: there exist $\tilde{C}>0$ such that, for all probability measures $\nu,\nu^{\prime}\in\mathcal{P}_{\mathcal{V}}$,

[TABLE]

and in particular

[TABLE]

Then for all $n\in\mathbb{N}$, considering the Euclidian division $n=km+r$ we get that

[TABLE]

The equivalence betwenn $\mathcal{V}$ and $\mathcal{W}$, hence between $\|\cdot\|_{\mathcal{V}}$ and $\|\cdot\|_{\mathcal{W}}$, concludes. ∎

4.4.4 Asymptotic theorems

Theorem 12**.**

Suppose that $U$ admits a strict local minimum and that there exist $R,h>0$ such that (9) holds for all $i\in\llbracket 1,d\rrbracket$, $x\in\mathbb{Z}^{d}$ and $v\in\{-1,1\}^{d}$. Let $f:\mathbb{Z}^{d}\rightarrow\mathbb{R}$ be such that $\|f/\mathcal{W}\|_{\infty}<\infty$, where $\mathcal{W}$ is defined in Theorem 11 (here and below we identify $f$ with the function $(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}\mapsto f(x)$). Consider the Zig-Zag walk on $\mathbb{Z}^{d}\times\{-1,1\}^{d}$ associated to $U$ with some initial condition $(x,v)$. Then, almost surely,

[TABLE]

If, moreover, $\|f/\sqrt{\mathcal{W}}\|_{\infty}<\infty$, then

[TABLE]

for some $\sigma_{f}\geqslant 0$, where $(B_{t})_{t\geqslant 0}$ is a one-dimensional Brownian motion.

Proof.

For the first part of the Theorem, simply decompose

[TABLE]

From Proposition 11 and the law of large numbers for $\mathcal{W}$-regular ergodic Markov chains, these terms almost surely converge respectively to $\mu_{\sigma(x,v)}(f)/2$ and $\mu_{-\sigma(x,v)}(f)/2$ (since $(X_{1},V_{1})\in\mathcal{A}_{-\sigma(x,v)}$ almost surely), which are both equal to $\pi(f)/2$.

For the second part, note that by the Jensen inequality and Proposition 10,

[TABLE]

for some $\gamma\in(0,1)$ and $C>0$. Hence, $\sqrt{\mathcal{V}}$ is still a Lyapunov function for $Q$ and the results established with $\mathcal{V}$ and $\mathcal{W}$ in the previous section also hold with $\sqrt{\mathcal{V}}$ and $\sqrt{\mathcal{W}}$. Let $f:\mathbb{Z}^{d}\times\{-1,1\}^{d}\rightarrow\mathbb{R}$ with $\|f/\sqrt{\mathcal{V}}\|_{\infty}<\infty$. From Theorem 11 (applied with $\sqrt{\mathcal{W}}$), for all $(x,v)\in\mathbb{Z}^{d}\times\{-1,1\}^{d}$ and all $n\in\mathbb{N}$,

[TABLE]

and, using (13),

[TABLE]

Here we used that $\mu_{s}Q=\mu_{-s}$, which can be obtained from

[TABLE]

where the limit holds in $\mathcal{P}_{\mathcal{W}}$ thanks to Proposition 11 together with the fact that the support of $\mu_{s}Q$ is included in $\mathcal{A}_{-s}$. We have thus obtained that

[TABLE]

is well-defined and satisfies $\|g/\sqrt{\mathcal{W}}\|_{\infty}\leqslant C\|f/\sqrt{\mathcal{W}}\|_{\infty}$ for some $C>0$ independent from $f$. Now suppose that in fact $f$ is a function of space alone, i.e. $f(x,v)=f(x)$. In that case $\mu_{s}(f)=\mu(f)$ for all $s\in\{-1,1\}^{d}$ and thus

[TABLE]

in other words $g$ is the Poisson solution associated to $Q$ and $f$. Since $\sqrt{\mathcal{W}}\in L^{2}(\mu)$, so does $g$, and [33, Theorem 3.1] concludes. ∎

5 Numerical scheme for hybrid kinetic samplers

5.1 The continuous-time processes

In this section we consider a class of kinetic processes for MCMC that can have a jump, drift and/or diffusion component at the same time in their dynamics.They are defiend as follows.

Let $U\in\mathcal{C}^{\infty}(\mathbb{T}^{d})$ (where $\mathbb{T}=\mathbb{R}/\mathbb{Z}$ is the $1$-periodic torus) and denote by $\mu$ the Gibbs measure on $E=\mathbb{T}^{d}\times\mathbb{R}^{d}$ associated to the Hamiltonian $H(x,v)=U(x)+|v|^{2}/2$, namely the probability law on $E$ with density proportional to $\exp(-H)$. Suppose that $\nabla U(x)=\sum_{i=0}^{N}F_{i}(x)$ where $N\in\mathbb{N}$ and $F_{i}\in\mathcal{C}^{\infty}(\mathbb{T}^{d},\mathbb{R}^{d})$ for all $i\in\llbracket 0,N\rrbracket$. Consider the operator $L$ defined for all $f\in\mathcal{C}^{2}(E)$ by

[TABLE]

where

[TABLE]

and where $\gamma,\lambda\geqslant 0$, $\nu_{d}$ is the standard $d$-dimensional Gaussian distribution and

[TABLE]

is the orthogonal reflection of $v$ with respect to $F_{i}(x)$. We call respectively $A_{1}$ the transport operator, $A_{2}$ the drift one, $A_{3,i}$ the $i^{th}$ bounce one, $A_{4}$ the Ornstein-Uhlenbeck (or friction/dissipation) one and $A_{5}$ the refreshment one.

As particular cases, many usual kinetic processes used in MCMC algorithms can be recovered:

• •

$F_{0}=\nabla_{x}U$ and $\lambda=\gamma=0$ corresponds to the Hamiltonian dynamics.

• •

$F_{0}=\nabla_{x}U$, $\gamma>0$ and $\lambda=0$ to the Langevin diffusion.

• •

$F_{0}=\nabla_{x}U$, $\gamma=0$ and $\lambda>0$ to the Hybrid Monte Carlo (HMC) algorithm.

• •

$F_{0}=0$, $N=1$, $F_{1}=\nabla_{x}U$, $\gamma=0$ and $\lambda>0$ to the Bouncy Particle sampler.

• •

$F_{0}=0$, $N=d$, $F_{i}=\nabla_{x_{i}}Ue_{i}$ (recall $e_{i}$ denotes the $i^{th}$ vector of the canonical basis of $\mathbb{R}^{d}$), $\gamma=\lambda=0$ to the Zig-Zag process.

We could also consider other kinds of jump mechanisms, like the randomized bounces of [49, 36], or different kinds of relaxation operators in the velocity operator rather than $A_{4}$ and $A_{5}$, for instance refreshment of velocities coordinate by coordinate, or partial refreshments for which, at exponential random times with parameter $\lambda$, the velocities $v$ jumps to $(1-\alpha)v+\sqrt{\alpha}G$ where $G\sim\nu_{d}$ (varying $\alpha\in(0,1]$ and $\lambda$ interpolates between $A_{5}$ and $A_{4}$, the latter being the limit $\alpha\rightarrow 0$ and $\lambda\rightarrow+\infty$ with $\lambda\alpha=1$). However, this would just make the notations heavier and the presentation more confused, without adding any particularly new idea with respect to the discussion to come, and thus we stick to (14).

For MCMC purposes, the law of the process associated to a generator $L$ of the form above should converge in large times toward the target measure $\mu$. This is established in the next two results. However, this is not the main motivation of this section, which is the study of the numerical sampling of such a process, and thus we only give sketches of proof and references for these results.

Proposition 13**.**

The set of compactly supported smooth functions $\mathcal{C}^{\infty}_{c}(E)$ is a core for $L$ and $\mu$ is invariant for $L$.

Proof.

The proof that $\mathcal{C}^{\infty}_{c}(E)$ is a core for $L$ is similar to the case of the Bouncy Particle Sampler in [15, Theorem 21 and Proposition 23]. From that, the invariance of $\mu$ straightforwardly follows from the fact that $\mu(Lf)=0$ for all compactly supported $\mathcal{C}^{2}$ functions, which is easily checked through integration by part (see also [39, Section 1.4]). ∎

Similarly to Section 4.4.3, we denote $\mathcal{W}(x,v)=1+|v|^{2}$, $\mathcal{P}_{\mathcal{W}}=\{\nu\in\mathcal{P}(E),\ \nu(\mathcal{W})<\infty\}$, $\|\nu_{1}-\nu_{2}\|_{\mathcal{W}}=\sup_{|f|\leqslant\mathcal{W}}|\nu_{1}f-\nu_{2}f|$.

Proposition 14**.**

Suppose that either $\gamma>0$ or $\lambda>0$. Then there exist $\kappa,C>0$ such that the semi-group $(P_{t})_{t\geqslant 0}$ associated to $L$ satisfies, for all $t\geqslant 0$ and $\nu\in\mathcal{P}_{\mathcal{W}}$,

[TABLE]

Proof.

Similarly to Theorem 11, this is classically obtained by combining a Doeblin and a Foster-Lyapunov conditions [24].

When $\gamma>0$ or $\lambda>0$, the Doeblin condition is obtained through controllability arguments (a notable fact is that under quite general conditions on $U$ the Zig-Zag process is irreducible even if $\lambda=\gamma=0$, see [6]). Remark that the bounce operators play no role here since if, uniformly in $(x,v)$ in a compact of $E$, $\delta_{(x,v)}e^{t_{0}(A_{1}+A_{2}+\gamma A_{4}+\lambda A_{5})}$ is bounded below at some time $t_{0}>0$ by $c_{0}$ times the uniform measure on some compact set with some $c_{0}>0$, then $\delta_{(x,v)}e^{t_{0}L}$ is bounded below by $c_{0}e^{-t_{0}\sum_{i=1}^{N}\|F_{i}\|_{\infty}}$ times this uniform measure, where we have used that the total bounce rate is bounded above by $\sum_{i=1}^{N}\|F_{i}\|_{\infty}$ so that there is a probability at least $e^{-t_{0}\sum_{i=1}^{N}\|F_{i}\|_{\infty}}$ that no bounce occurs on the time interval $[0,t_{0}]$.

The Lyapunov condition stems from the fact that for all $(x,v)\in E$,

[TABLE]

and thus $L\mathcal{W}\leqslant-\min(2\gamma,\lambda)/2\mathcal{W}+C$ for some $C>0$.

∎

As announced, from now and in the rest of Section 5 we focus on the question of the the numerical sampling of the process with generator $L$.

When $F_{0}=0$, and $\gamma=0$, the process is a piecewise deterministic velocity jump process. Between two random jumps, the process simply follows the flow $(x,v)\mapsto(x+tv,v)$, so that the jump time $T_{i}$ associated to the vector field $F_{i}$ follows the law

[TABLE]

Provided that for all $x,v\in\mathbb{R}^{d}$ and $s>0$, $\left(v\cdot F_{i}(x+sv)\right)_{+}\leqslant\varphi_{x,v}(s)$ for some function $\varphi_{x,v}$ such that $\int_{0}^{t}\varphi_{x,v}(s)\text{d}s$ can be computed, the continuous-time thinning algorithm allows for an exact simulation of the jump times, hence of the process [32]. The absence of discretization bias on the invariant measure is an argument in favour of these kinetic processes. This is no longer the case when $F_{0}\neq 0$, except in very particular cases (e.g. harmonic oscillators) where the ODE $\partial_{t}(x,v)=(v,-F_{0}(x))$ can be explicitly solved. In most cases, this ODE is solved numerically and thus the simulation for the stochastic process is not exact. In order to conserve the invariant measure and suppress the bias, a Metropolis step can be added [7], but this slows down the motion of the process, increases the variance and may thus be counter-productive, as observed when comparing Metropolis Adjusted Langevin algorithm and Unadjusted Langevin Algorithm [16].

So why mix a deterministic drift and jump mechanisms if this prohibits exact simulation? The motivation is given by the possible numerical gain given by thinning, as presented in Section 3.2. Indeed, as said before, exact simulation requires a bound on the jump rate, and a poor (i.e. large) bound leads to many jump proposals per time unit and a very low efficiency. In particular, jump mechanisms are not adapted for fastly-varying potentials, like potentials used in molecular dynamics with a singularity at zero (Lennard-Jones, Coulomb…). On the other hand, jump mechanisms are very efficient with Lipschitz potentials. So, a mixed drift/jump part is interesting as soon as the potential exhibits different scales, like fast-varying but numerically cheap parts together with Lipschitz but numerically intensive parts. This is similar to the idea of multi-time-steps algorithms [47, 11] with somehow random adaptive time-steps, except that now we are simply going to discretize (with a unique time-step, no subtlety here) a continuous-time process which is ergodic with respect to the target law and thus there shouldn’t be any resonance problem as exhibited by multi-time-steps algorithms. Besides, for the applications in molecular dynamics that have motivated this question [42], due to stability issues raised by very fast oscillations in some parts of the system, the time-step is anyway constrained to be very small, as compared to a high variance that comes from the problem of exploring a complex, high-dimensional, multi-scale, metastable landscape. So, exact simulation is not necessarily our objective.

5.2 A Strang splitting scheme

A Strang splitting scheme to compute the evolution given by a generator $L=L_{1}+L_{2}$ is based on the fact that, formally,

[TABLE]

Hence, if we can simulate exactly a process with generator $L_{1}$ and $L_{2}$, or more generally if we have second-order approximations of those, we get a second-order scheme for $L_{1}+L_{2}$. Using twice this fact,

[TABLE]

In particular, from the considerations developed in Section 3.3 (and Theorem 5), the equilibrium of the corresponding Markov chain is close (in some senses) to the Gibbs measure $\mu$ at order $\delta^{2}$, where $\delta$ is the timestep.

For instance, for the Langevin diffusion, which corresponds to $L$ given by (14) with $\gamma>0$, $\lambda=0$, $N=0$ and $F_{0}=\nabla_{x}U$, many splitting schemes can be considered, see e.g. the discussions in [29, 7]. A very precise study, both theoretical and empirical, of all the possible schemes for mixed jump/diffusion processes is beyond the scope of the present paper, and we will only consider one particular choice. Consider the splitting

[TABLE]

Each of the three evolutions corresponding to $e^{tL_{i}}$, $i=1,2,3$ can be sampled exactly (remark that, in particular cases, the velocity jump process corresponding to $e^{t(L_{1}+L_{3})}$ could also be sampled exactly). As a consequence, for a given time-step $\delta>0$, we consider the kinetic walk $(X_{n},V_{n})_{n\in\mathbb{N}}$ whose transition is defined as follows:

Set $\widetilde{X}_{n}=X_{n}+\delta V_{n}/2$. 2. 2.

If $\gamma>0$, set

[TABLE]

with $G$ a standard Gaussian random variable. If $\gamma=0$, set $\widetilde{V}_{n}=V_{n}-\delta F_{0}(\widetilde{X}_{n})/2$. 3. 3.

Set $\widehat{V}_{n}=W_{\delta}$ where $(Y_{t},W_{t})_{t\in[0,\delta]}$ is a Markov chain with generator $L_{3}$ and initial condition $(Y_{0},W_{0})=(\widetilde{X}_{n},\widetilde{V}_{n})$ (so that $Y_{t}=\widetilde{X}_{n}$ for all $t\in[0,\delta]$). 4. 4.

If $\gamma>0$, set

[TABLE]

with $G^{\prime}$ a standard Gaussian random variable. If $\gamma=0$, set $V_{n+1}=\widehat{V}_{n}-\delta F_{0}(\widetilde{X}_{n})/2$. 5. 5.

Set $X_{n+1}=\tilde{X}_{n}+\delta V_{n+1}/2=X_{n}+\delta(V_{n}+V_{n+1})/2$.

(Of course, implicitly, the variables $G$, $G^{\prime}$ and $(W_{t})_{t\in[0,\delta]}$ are independent one from the other and from the past trajectory). Then, denoting $Q_{\delta}$ the transition operator of $(X_{n},V_{n})$ defined by $Q_{\delta}f(x,v)=\mathbb{E}\left(f(X_{1},V_{1})\ |\ (X_{0},V_{0})=(x,v)\right)$ for all measurable bounded function $f$, one has

[TABLE]

To make the link with the discussion of Section 3, for all measurable bounded $f$,

[TABLE]

where $p_{\delta}$ is defined by the fact that, in the algorithm above, for all $n\in\mathbb{N}$ the law of $V_{n+1}$ is $p_{\delta}(X_{n},V_{n};\cdot)$.

Proposition 15**.**

Let $(Y_{t},W_{t})_{t\geqslant 0}$ be a Markov process associated to $L$ and, for all $\delta>0$, let $(X_{n}^{\delta},V_{n}^{\delta})_{n\geqslant 0}$ be a Markov chain with transition operator $Q_{\delta}$ and initial condition $(X_{0}^{\delta},V_{0}^{\delta})=(Y_{0},W_{0})$. Then

[TABLE]

Proof.

We start by a smoothing/truncation step similar to the study in [15, Theorem 21]. For all $\varepsilon>0$ small enough and all $i\in\llbracket 1,N\rrbracket$, we consider the jump rate on $E$ given by

[TABLE]

and the associated regularized bounce operator

[TABLE]

We also consider $\nu_{d,\varepsilon}$ the standard Gaussian law conditioned on $\{|v|\leqslant 1/\varepsilon\}$ and the corresponding truncated refreshment operator

[TABLE]

Then, $L_{\varepsilon}:=A_{1}+A_{2}+\sum_{i=1}^{N}A_{3,i,\varepsilon}+\gamma A_{4}+\lambda A_{5,\varepsilon}$ is such that

[TABLE]

for some $C>0$. The two processes $(Y_{t},W_{t})_{t\geqslant 0}$ and $(Y_{t}^{\varepsilon},W_{t}^{\varepsilon})_{t\geqslant 0}$ with respective generators $L$ and $L_{\varepsilon}$ can be defined simultaneously following the synchronous coupling detailed in [15, Section 6], in such a way that $(Y_{t},W_{t})=(Y_{t}^{\varepsilon},W_{t}^{\varepsilon})$ up to a random time that is stochastically bounded above by an exponential variable with parameter $C\varepsilon$. In other words, for all $t\geqslant 0$,

[TABLE]

In particular $(Y_{t}^{\varepsilon},W_{t}^{\varepsilon})_{t\geqslant 0}\overset{law}{\underset{\varepsilon\rightarrow 0}{\longrightarrow}}(Y_{t},W_{t})_{t\geqslant 0}$. The same coupling argument holds for a chain $(X_{n}^{\delta,\varepsilon},Y_{n}^{\delta,\varepsilon})_{n\in\mathbb{N}}$ with transition operator $Q_{\delta,\varepsilon}$ defined like $Q_{\delta}$ but after smoothing/truncation. It is thus sufficient to prove Proposition 15 in the case where $A_{5}$ and $(A_{3,i})_{\in\llbracket 1,N\rrbracket}$ are replaced by the smooth and truncated operators.

Denote $\mathcal{B}=\mathrm{vect}\left(\{A_{1},A_{2},A_{4},A_{5,\varepsilon}\}\cup\{A_{3,i,\varepsilon}\ :\ i\in\llbracket 1,N\rrbracket\}\right)$. If $f\in\mathcal{C}_{c}^{\infty}(E)$ then $Bf\in\mathcal{C}_{c}^{\infty}(E)$ for all $B\in\mathcal{B}$ (this is false without the smoothing/truncation step, which is the reason why it has been added). Since $\mathcal{C}_{c}^{\infty}(E)$ is in the strong domain of all operators of $\mathcal{B}$, we can use for $B_{1},B_{2}\in\mathcal{B}$ the decomposition

[TABLE]

and the fact $\|e^{tB}f\|_{\infty}\leqslant\|f\|_{\infty}$ for all $B\in\mathcal{B}$ to get that

[TABLE]

for all $f\in\mathcal{C}_{c}^{\infty}(E)$. Using this repeatedly, we get that

[TABLE]

for all $f\in\mathcal{C}_{c}^{\infty}(E)$, and [27, Theorem 17.28] concludes. ∎

5.3 Computational complexity

Let us informally discuss the efficiency of the algorithm introduced in the previous section that defines the transition associated with the operator $Q_{\delta}$. It should be compared with the transition of a numercial scheme for classical processes, like Langevin or Hamiltonian dynamics. The numerical cost mainly comes from the computation of the forces. For the Langevin or Hamiltonian dynamics, at each time step, $\nabla U$ is computed once. Similarly, for a transition associated with $Q_{\delta}$, $F_{0}$ is evaluated once per step, at $x=\widetilde{X}_{n}$. What about $(F_{i})_{i\in\llbracket 1,N\rrbracket}$?

First, consider a naive construction of the third step of the algorithm, namely the construction of $W_{\delta}$ when $(Y_{t},W_{t})_{t\geqslant 0}$ is a continuous-time Markov chain with generator $L_{3}$ and initial condition $(x,w)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$. For all $t\geqslant 0$, $Y_{t}=x$. Start by sampling a Poisson process with intensity $\lambda$ on $[0,\delta]$, consider $T_{0}$ its last jump in $[0,\delta]$ (with $T_{0}=0$ if there is no jump, in particular if $\lambda=0$). If $T_{0}=0$, set $W_{T_{0}}=w$, else draw $W_{T_{0}}$ according to the standard Gaussian distribution. Suppose by induction that $T_{n}$ and $W_{T_{n}}$ have been defined for some $n\in\mathbb{N}$. Let $(E_{n,i})_{i\in\llbracket 1,N\rrbracket}$ be i.i.d. exponential random variable and let

[TABLE]

Let $i_{n}\in\llbracket 1,N\rrbracket$ be such that $T_{n,i_{n}}=T_{n+1}$ ($i_{n}$ is almost surely unique), set $W_{t}=W_{T_{n}}$ for all $t\in[T_{n},T_{n+1})$ and $W_{T_{n+1}}=R(x,W_{T_{n}})$. Remark that the norm of $W_{t}$ is conserved at a jump time, so that the jump rate $\sum_{i=1}^{N}(W_{t}\cdot F_{i}(x))_{+}$ is bounded and $\max\{n\in\mathbb{N}\ :\ T_{n}<\delta\}$ is almost surely finite, and thus $W_{\delta}$ is defined after a finite number of jumps.

This is a correct construction, but it relies on the computation of $F_{i}(x)$ for all $i\in\llbracket 1,N\rrbracket$. In other words, with this naive construction, there is no numercial gain: $\nabla U(x)=\sum_{i=0}^{N}F_{i}(x)$ is computed at each time step. As a consequence, the algorithm is only useful if a relevant thinning procedure such as discussed in Section 4.3 is available, to avoid the systematic computation of $(F_{i}(x))_{i\in\llbracket 1,N\rrbracket}$ at each step. Such a thinning relies on suitable bounds on the vector fields $(F_{i})_{i\in\mathbb{N}}$, that depends on the form of $U$ and of the way $\nabla U$ is splitted. For this reason, in the rest of this section, we won’t address this question in a general case, but rather focus on a particular example. This will illustrate the fact that there exist cases where splitting forces between drift and jump mechanisms can decrease the total computational cost of the simulation.

5.3.1 A motivating example

To fix ideas, consider a system of $M$ particles in the torus $(a\mathbb{T})^{3}$ interacting through truncated a Lennard-Jones potential, in other words, given some parameters $U_{0},r,a>0$, the total energy of the configuration $x\in((a\mathbb{T})^{3})^{M}$ is

[TABLE]

where $W(h)=[(r/h)^{12}-(r/h)^{6}]\chi(h/a)$ with $\chi$ a $\mathcal{C}^{2}$ positive function with values in $[0,1]$ such that $\chi(s)=1$ for $s\leqslant 1/2$ and $\chi(s)=0$ for $s\geqslant 1$ (in particular, a particle doesn’t interact with its periodic image, or with several copies of the same other particle). Remark that strictly speaking this problem doesn’t enter the framework considered above where, for the sake of simplicity, $U$ was supposed smooth, but this is typically the kind of singular potentials met in molecular dynamics simulations.

For $i\in\llbracket 1,M\rrbracket$, denote $J_{i}$ the $3M\times 3$ matrix with zeros everywhere except $J_{i}(3(i-1)+1,1)=J_{i}(3(i-1)+2,2)=J_{i}(3(i-1)+3,3)=1$, so that $\nabla U=\sum_{i=1}^{M}J_{i}\nabla_{x_{i}}U$. For some $R\ll a$, we decompose $\nabla U=F_{0}+\sum_{i=1}^{M}\sum_{j\neq i}J_{i}G_{i,j}$ with

[TABLE]

Then $F_{0}$ gathers the (singular) short-range forces, and the $G_{i,j}$’s the (bounded) long-range ones. Set

[TABLE]

In particular the number of jump mechanisms is $N=M(M-1)$. Remark that the jump rate associated with the vector field $F_{i,j}$ can be bounded as $(v\cdot J_{i}G_{i,j}(x))_{+}\leqslant|v_{i}|U_{0}\sup\{|W^{\prime}(s)|,\ s\geqslant R/2\}:=|v_{i}|C_{R}$. As $R$ increases (keeping $R\ll a$), $C_{R}$ decays very fast toward zero (more precisely, as $R^{-5}$).

Conversely, computing $F_{0}(x)$ involves computing, for each $i\in\llbracket 1,N\rrbracket$, a sum over all the particles $j$ such that $\{|x_{i}-x_{j}|\leqslant R\}$. The numerical cost of this computation increases with $R$, but is small with respect to the computation of $\nabla U(x)$ as long as $R\ll a$.

5.3.2 Thinning

For the Lennard-Jones system with the decomposition $\nabla U=F_{0}+\sum_{i=1}^{M}\sum_{j\neq i}F_{i,j}$ introduced above, the third step of the algorithm of Section 5.2 can be achieved as follows:

• •

Set $W=\widetilde{V}_{n}$.

• •

For all $i\in\llbracket 1,M\rrbracket$,

• •

Draw $K_{i}$ according to a Poisson law with parameter $|W_{i}|C_{R}M\delta$.

• •

For all $k\in\llbracket 1,K_{i}\rrbracket$,

• •

Draw $j$ uniformly over $\llbracket 1,M\rrbracket$ and $U_{i,k}$ uniformly over $[0,1]$.

• •

If $j\neq i$ and $U_{i,k}\leqslant\left(W_{i}\cdot G_{i,j}(\widetilde{X}_{n})\right)_{+}/|W_{i}|C_{R}$, do

[TABLE]

else do nothing.

• •

end for all $k$

• •

end for all $i$.

• •

Set $\widehat{V}_{n}=W$.

We call this the thinned algorithm. Remark that $|W_{i}|$ is unchanged by the reflection $R_{i,j}$, which allows to draw $K_{i}$, the number of jump proposals for the $i^{th}$ particle, at the beginning of the loop (but in practice this is not important). More crucially, note also that each particle $i\in\llbracket 1,M\rrbracket$ can be treated in parallel, and there is no need for time synchronization: indeed, two particles $i$ and $j$ only interact through $G_{i,j}(\widetilde{X}_{n})$, and the positions $\widetilde{X}_{n}$ are fixed at this step, so that the velocity jumps of each particle doesn’t affect the law of the other.

5.3.3 Numerical efficiency

Computing $\nabla U$ has a numerical cost of order $\mathcal{O}(M^{2})$, and it is computed $T/\delta$ times if we sample a trajectory of a classical process (Langevin dynamics, etc.) in a time interval $[0,T]$ with a usual integrator with time-step $\delta$. Let us compare this with the cost of computing gradients in the thinned algorithm for the hybrid process.

First, denote $\mathcal{N}(R)$ the average number of particles that are at distance less than $R$ from a given particle. By using a Verlet list of neighbors, computing $F_{0}$ has an average cost of $\mathcal{O}\left(M\mathcal{N}(R)\right)$.

Second, denote $D_{T}$ the number of times that $G_{i,j}$ has been evaluated for some $i,j\in\llbracket 1,M\rrbracket$ in a trajectory of length $T$ with time-step $\delta$. Then

[TABLE]

where, conditionally to the $\widetilde{V}_{n,i}$’s the $S_{n,i}$’s are independent Poisson random variables with parameter $|\widetilde{V}_{n,i}|C_{R}M\delta$. By the law of large numbers for ergodic Markov chains,

[TABLE]

where $H$ is the average of $|w_{i}|$ with respect to the equilibrium distribution of the chain. When $\delta$ is small, the velocity distribution at equilibrium is close to a standard Gaussian one of dimension 3, so that $H\leqslant\sqrt{3}+\mathcal{O}(\delta)$.

Computing $G_{i,j}(x)$ being of order $\mathcal{O}(1)$, the total cost is $\mathcal{O}(M^{2}C_{R}T+M\mathcal{N}(R)T/\delta)$. As noted previously, as $R$ increases, $C_{R}$ decays but $\mathcal{N}(R)$ increases. In the regime where $\mathcal{N}(R)\ll M$, the gain in the thinned algorithm, with respect to an integrator where the full gradient $\nabla U$ is computed at each step, is the ratio of $M^{2}C_{R}T$ and $M^{2}T/\delta$, namely $C_{R}\delta$. In other words, the long-range forces are only evaluated at an average time-step $1/C_{R}$. Depending on the parameters, this can be a significant speed-up. Indeed, $\delta$ is typically constrained to be small in order to handle the singular behaviour of short-range forces while, for intermediate values of $R$, $C_{R}$ can already be small. We refer to [42] for a practical case where $\delta C_{R}\simeq 10^{-6}$.

5.3.4 Sampling rate: the mean-field regime

The numerical efficiency discussed above is encouraging, but one should be careful when comparing two Markov processes for sampling. Indeed, suppose we are given two kinetic processes $(X_{t},V_{t})_{t\in[0,T]}$ and $(\tilde{X}_{t},\tilde{V}_{t})_{t\in[0,T]}$, both $\mu$-ergodic, such that simulating a trajectory of length $T$ for the first one is $10$ times faster than for the seconde one. This is useless if the first process explores the space $100$ times slower than the second one so that, to achieve a similar quality of sampling, much longer trajectories are required.

Such a situation may be feared for the hybrid drift/jump process introduced above in the Lennard-Jones model. Indeed, under the effect of many collisions, the process may show as $M\rightarrow\infty$ a diffusive behaviour, namely the jumps may average and give a velocity close to zero. In that case, the exploration of the space (hence the convergence toward equilibrium) would be very slow. Let us (informally) check that it is not the case here, at least in the mean-field regime: in the following, we suppose that $U_{0}=1/M$ (in which case the constant $C_{R}$ is of order $1/M$).

Consider the Langevin dynamics with generator $L=A_{1}+A_{2}+A_{4}$ with $F_{0}=\nabla U$. Suppose that the initial conditions $(x_{i},v_{i})_{i\in\llbracket 1,M\rrbracket}$ are i.i.d. with the $x_{i}$’s distributed according to some law $\nu$. Conditionally to $(x_{1},v_{1})$, by the law of large numbers, the force felt by the first particle at time [math] satisfies

[TABLE]

Classical propagation of chaos results show that, as $M\rightarrow+\infty$, the particles behaves approximately as $M$ independent processes on $(a\mathbb{T})^{3}\times\mathbb{R}^{3}$ with generator

[TABLE]

where $\nu_{t}$ is the law of the process at time $t$, solution of the non-linear equation $\partial_{t}(\nu_{t}f)=\nu_{t}L_{t}f$ for all nice $f$.

Similarly, for the hybrid jump/diffusion process obtained from the splitting of the forces introduced in Section 5.3.1, conditionally to $(x_{1},v_{1})$, by the law of large numbers, denoting $g$ the function such that $G_{i,j}(x)=g(x_{i}-x_{j})/M$ for $i,j\in\llbracket 1,M\rrbracket$, the jump rate at time 0 for the first particle satisfies

[TABLE]

and a similar limit $Q_{\nu}$ is obtained for the jump kernel at time [math]. From the propagation of chaos phenomenon, the system is expected to behave as $M\rightarrow\infty$ as $M$ independent non-linear processes on $(a\mathbb{T})^{3}\times\mathbb{R}^{3}$ with non-homogeneous generators given by

[TABLE]

with a suitable function $\tilde{g}$, and $\nu_{t}$ the law of the process (see [41] and references within). The non-homogeneous Markov process with generator $L_{t}$ is not degenerated in the sense its velocity is not averaged to zero, which means that for large $M$ the system of interacting particles is not expected to have a diffusive behaviour, nor to converge to equilibrium with a rate that vanishes as $M\rightarrow+\infty$.

Of course the dynamics (Langevin versus hybrid) are different and thus they may converge toward equilibrium at different speed. Nevertheless, the ratio of the convergence rate may be expected to be of order independent from $M$, since the mean-field limits as $M\rightarrow\infty$ are obtained with the same scaling (namely it was not necessary to accelerate time). In other words, in order to achieve the same convergence as the Langevin process in a time $T$ (for a cost of order $\mathcal{O}(M^{2}T/\delta)$), we may need to sample the hybrid process up to a time $T^{\prime}>T$ (for a cost of order $\mathcal{O}(M\mathcal{N}(R)T^{\prime}/\delta+M^{2}C_{R}T^{\prime})$) but with a ratio $T^{\prime}/T$ bounded below independently from $M$. For $R$ small enough (independent from $M$) so that $\mathcal{N}(R)<(T/T^{\prime})M$, using that $C_{R}$ is $\mathcal{O}(1/M)$, we see that the hybrid process appears as the most efficient for large $M$.

More generally, if $U_{0}$ is not of order $1/M$, the conclusion may be unclear, but the informal discussion above shows that there exist cases where the parameters of the problem are such that a suitable factorization can lower the numerical cost of the simulation at constant long-time convergence quality. We refer to [42] for a practical case where the diffusion constant, used as an indicator of the sampling rate, is of the same order for the two processes.

5.3.5 An alternative splitting of the forces

For the Lennard-Jones model, consider the splitting of the forces $\nabla U=\sum_{i=0}^{M}F_{i}$ where $F_{0}$ is as above and $F_{i}(x)=J_{i}\nabla_{x_{i}}U=J_{i}\sum_{j\neq i}G_{i,j}$ for $i\in\llbracket 1,M\rrbracket$. Since $(v_{i}\cdot\sum_{j\neq i}G_{i,j})_{+}\leqslant\sum_{j=1}^{M}(v_{i}\cdot G_{i,j})_{+}$, this alternative splitting reduces the total jump rate of the process, by comparison with the previous splitting. Increasing the jump rate is known to increase the diffusive behaviour and the asymptotical variance for the Zig-Zag process, at least in dimension $1$ [3]. Since the present case is similar, we may expect the alternative splitting to yield a process that converges faster toward equilibrium than the initial one. Nevertheless, in the mean-field case, similarly to the previous section, this speed-up should be independent from $M$.

On the other hand, let us estimate the numerical cost (in term of computations of forces) of the alternative process. We can bound $|F_{i}(x)|\leqslant MC_{R}$ (with the same $C_{R}$ as previously), and thus a thinning algorithm yields again an average of $MC_{R}T$ jump proposals per particle in a time period $[0,T]$. But now, at each jump proposal, $F_{i}(x)$ has to be computed, which has a cost $\mathcal{O}(M)$, so the total cost is $\mathcal{O}(M^{3}C_{R}T)$.

This little computation indicates that, in the case of forces that come from pairwise interactions between particles, rather than considering one jump mechanisms per particle (as in the usual Zig-Zag process, and in the alternative splitting of this section), it may be better to split $\nabla_{x_{i}}U$ so that each jump mechanism is only associated to a single interaction $(i,j)$ (as in the initial splitting introduced in Section 5.3.1).

Acknowledgements

The author would like to thank Gabriel Stoltz and Louis Lagardère for fruitfull discussions about Strang splitting schemes, and acknowledges partial financial support from the ERC grant MSMATHS (European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement number 614492) and the ERC grant EMC2, and from the French ANR grant EFI (Entropy, flows, inequalities, ANR-17-CE40-0030).