Limit theorems for the Zig-Zag process

Introduction

Markov Chain Monte Carlo methods remain an essential computational tool in statistics and among other things have made it possible for Bayesian inference techniques to be applied to increasingly complex models. Due to its simplicity and wide applicability, the Metropolis-Hastings (MH) algorithm and its numerous variants remain the most widely used MCMC method for sampling from a general target probability distribution, despite having been introduced over 60 years ago. Given a target distribution $\pi$, the Metropolis-Hastings scheme defines a discrete time Markov chain which will be both ergodic and reversible with respect to $\pi$. The fact that the Markov chain is reversible is a serious limitation. Indeed, it is now well known that non-reversible chains can significantly outperform reversible chains, in terms of rate of convergence to equilibrium , asymptotic variance as well as large deviation functionals . One particular approach to improving performance is to introduce a velocity/momentum variable and construct Markovian dynamics which are able to mixing more rapidly in the augmented state space. Such methods include Hybrid Monte Carlo (HMC) methods, inspired by Hamiltonian dynamics, and numerous generalisations. While the standard construction of HMC is reversible, it is straightforward to alter the scheme such that the resulting process is non-reversible .

is satisfied, so that the Zig-Zag process can be used to approximate expectations with respect to $\pi$. Two one-dimensional examples of the Zig-Zag process are displayed in Figure 1.

While the construction and finite-time behaviour of PDMPs is well understood , their use within the context of sampling has only recently been considered and is mostly unexplored. The first such occurrence of a MCMC scheme based on PDMP appeared in the computational physics literature and in one dimension coincides with the Zig-Zag sampler. This scheme was extended and analysed carefully in , where it was rechristened the Bouncy Particle Sampler. In one dimension, the quantitative long-time behaviour of related PDMP schemes has been analysed in detail, see for example . More recently in , the application of the Zig-Zag sampler to big data settings was investigated. It was found that the Zig Zag sampler lends itself very well to such problems since sub-sampling can be introduced without affecting the stationary distribution, as opposed to standard sub-sampling techniques, such as SGLD which are inherently biased. By introducing appropriate control variates a “super-efficient” sampling scheme for big data problems was produced, in the sense that it is able to generate independent samples from the target distribution at a higher efficiency than directly generating IID samples using the entire data set for each sample.

In this paper we seek to better understand the qualitative performance of the Zig Zag sampler. Focusing on the one-dimensional case, we study the important practical question of whether a central limit theorem (CLT) holds for the Zig-Zag process, i.e. whether for a given observable $f$,

where $\sigma^{2}_{f}$ is the asymptotic variance and where $\Rightarrow$ denotes convergence in distribution. Heuristically, once a CLT is known to hold, we know that the ergodic average in (1) converges at rate $\sigma_{f}/\sqrt{t}$, which is the best convergence to be expected in a Monte Carlo simulation. It is also clear that a smaller value of $\sigma_{f}>0$ implies a faster convergence of the ergodic averages. Without a CLT, convergence may be arbitrarily slow. Starting from the case of a unimodal target distribution and extending to more general cases, we obtain sufficient conditions for (2) to hold. Moreover, we identify conditions under with the CLT can be strengthened to an invariance principle or functional central limit theorem (FCLT) . For the one-dimensional Zig-Zag process we obtain explicit expressions for the asymptotic variance, which we illustrate for various examples.

Given a target distribution $\pi$, there is some freedom in choosing the switching rate $\lambda$ in such a way that $\pi$ is invariant for the Zig-Zag process. This freedom is crucial for the ability of the sub-sampling Zig-Zag scheme of to sample without bias. In Section 4 we study the influence of the particular choice of switching rate on the behaviour of the process. We show that as the switching rate is increased the Zig-Zag sampler will exhibit random walk behaviour. In particular, over an appropriate timescale the Zig-Zag sampler will behave asymptotically, as the excess switching rate tends to infinity, as an overdamped Langevin diffusion which is ergodic with respect to $\pi$.

As the Zig-Zag sampler is based upon a continuous time process, it is not immediately clear how its performance can be compared to existing discrete time sampling schemes. With this aim in mind, we derive approximations for the average switching rate of the process per unit time, and apply this to construct an effective sample size (ESS) for the Zig-Zag sampler which quantifies the number of independent samples generated in terms of the number of evaluations of the gradient of the log density. A suitable definition of effective sample size depends in an essential way on the asymptotic variance of the corresponding CLT, which further illustrates the importance of establishing a CLT from an applied viewpoint. Comparing to IID samples in some cases we observe a remarkable feature: the effective sample size of the Zig-Zag sampler will be larger than that of IID samples, behaviour which is strongly tied to the nonreversibility of the scheme.

We structure the paper as follows. In Section 2 we review the construction of the Zig-Zag sampler in the one dimensional case and explore its basic properties. Section 3 describes conditions for a CLT to hold for the one dimensional Zig-Zag sampler and characterises the asymptotic variance. These results are demonstrated numerically for some standard probability distributions. In Section 4 the diffusive regime is investigated where the switching rate $\lambda$ goes to infinity. Finally, in Section 5 an appropriate measure of effective sample size is introduced for the Zig-Zag sampler, and is used to compare the performance of the Zig-Zag sampler with other sampling techniques for some standard probability distributions. The proofs of most of results may be found in Appendix A. In Appendix B we discuss the simulation of the Zig-Zag process, which provides the necessary background for Section 5.

The Zig-Zag process

Furthermore for some $x_{0}>0$, we have $\lambda(x,\theta)>0$ if $\theta x\geq x_{0}$.

An alternative and convenient way of writing (3) is $\lambda(x,\theta)-\lambda(x,-\theta)=\theta U^{\prime}(x)$ for all $(x,\theta)\in E$. It is easy to check that (3) holds if and only if there exists a continuously differentiable function $U$ and a continuous non-negative function $\gamma$ such that

The switching rates $\lambda$ for which $\gamma\equiv 0$ are called canonical switching rates and the corresponding Zig-Zag process is called the canonical Zig-Zag process.

which will service as the generator of the Markov semigroup of the Zig-Zag process, with dynamics as discussed in the introduction. In the following proposition, the notion of ‘petite sets’ can be found in .

Suppose Assumption 2 holds. Then $(L,\mathcal{D}(L))$ is the extended generator of a piecewise deterministic Markov-Feller process $(Z(t))_{t\geq 0}:=(X(t),\Theta(t))_{t\geq 0}$ in $E$. All compact sets are petite for $(X(t),\Theta(t))$. Finally $\mu$ is the unique invariant probability distribution for $(Z(t))_{t\geq 0}$.

The proof of this result is located in Appendix A.1.

Central Limit Theorems for the Zig-Zag process

First, in Section 3.1, we obtain a CLT for the Zig-Zag process in the simple and intuitive case in which the target distribution is unimodal and the excess switching rate $\gamma=0$. Then we describe a general approach to the CLT in Section 3.2. We then illustrate the theory with several examples in Section 3.3.

If the potential $U(x)$ is continuously differentiable and is monotonically non-decreasing (non-increasing) for $x\geq 0$ ($x\leq 0$) then the canonical switching rates associated with $U$ satisfy $\lambda(x,+1)=0$ for $x\leq 0$, and $\lambda(x,-1)=0$ for $x\geq 0$. In this situation trajectories of the canonical Zig-Zag process will always pass through the origin $x=0$ between switches. This regular behaviour makes it possible to obtain a Central Limit Theorem in a very straightforward way: by inspecting the contributions towards the total variance of trajectory segments between crossings of the origin.

$\lambda(x,\theta)$ are the canonical switching rates defined by $\lambda(x,\theta)=(\theta U^{\prime}(x))^{+}$.

Note that the definition of $\pi$ agrees with the definition of $\pi$ below Assumption 2. Furthermore, the fact that $\exp(-U(x))$ is integrable, combined with the monotonicity assumption, implies that the switching rates $\lambda(x,\theta)$ are positive for $\theta x\geq x_{0}$, for some fixed $x_{0}>0$, so that Assumption 3.1 implies Assumption 2.

Suppose Assumption 3.1 holds. Let $(X(t),\Theta(t))$ denote the Zig-Zag process with switching rates $\lambda(x,\theta)$. Then

See Figure 2 for a graphical illustration of these times.

Now for $i=1,2,\ldots$, define the random variables

Let $N(t):=\sup\{i:T_{i}^{+}\leq t\}$. Then

Note that $(Y_{i})$ are i.i.d., with distribution identical to that of the random variable $Y:=Y^{+}+Y^{-}$, where $Y^{+}$ and $Y^{-}$ are independent random variables defined by

Also by this assumption, $\int_{0}^{T_{0}^{+}}g(X(s))\ ds$ and $\int_{T_{N(t)}^{+}}^{t}g(X(s))\ ds$ are bounded in probability. Furthermore

Combining all terms gives the stated expression for the asymptotic variance.

2 General approach to the Central Limit Theorem

The approach of Section 3.1 is intuitively appealing. However the required assumptions are very restrictive. In this section we will employ a far more general approach to obtaining a CLT. In particular, this approach allows us to include non-unimodal cases, as well as situations in which the excess switching rate $\gamma$ in (4) is non-zero.

First we recall two key results from the literature which will be helpful for our purposes. Recall the definition of a petite set from e.g. .

$(Z(t))_{t\geq 0}$ is a $\varphi$-irreducible continuous time Markov process in a Borel space $E$ with extended generator $L$. For a function $f:E\rightarrow[1,\infty)$, a petite set $C\in\mathcal{B}(E)$, a constant $b<\infty$ and a function $V:E\rightarrow[0,\infty)$, $V\in\mathcal{D}(L)$,

Suppose that Assumption 3.2 is satisfied. Then $(Z(t))_{t\geq 0}$ is positive Harris recurrent with invariant probability distribution $\mu$ and $\mu(f)<\infty$. For some $c_{0}<\infty$ and any $|g|\leq f$, the Poisson equation

admits a solution $\phi$ satisfying the bound $|\phi|\leq c_{0}(V+1)$.

The following general result establishes sufficient conditions for a functional Central Limit Theorem to hold. Part of the results in this section can be obtained simply by verifying the conditions of the following theorem, although in particular work needs to be done to find suitable functions $f$ and $V$ satisfying Assumption 3.2.

In situations where $\mu(V^{2})<\infty$ can not be established, we will have to establish a weaker (non-functional) form of the central limit theorem, which will depend on a CLT for martingales such as [21, Theorem 2.1]. We require the following lemmas, the proofs of which may be found in Appendix A.2.

Suppose Assumption 3.2 is satisfied. Let $g\in\mathcal{M}(E)$ be measurable, satisfy $|g|\leq f$ and $\mu(g)=0$. Suppose $\phi$ is a solution to the Poisson equation (8) for the generator $L$ given by (5) and suppose $\mu(|\phi|)<\infty$. Define the process

Suppose Assumption 3.2 is satisfied for the Zig-Zag process with generator (5) and let $g\in\mathcal{M}(E)$ satisfy $|g|\leq f$ and $\mu(g)=0$. Furthermore suppose $V$ satisfies $\mu(V)<\infty$, or alternatively $\mu(|\phi|)<\infty$ where $\phi$ is the solution of the Poisson equation given by Proposition 3.3. Let $\psi$ be given by

The stated result now follows by combining the obtained limits in (9).

We have now obtained two different expressions for the asymptotic variance, namely (6) and (12). In cases where both Theorem 3.1 and Theorem 3.7 apply these expression of course have the same value. In Appendix A.3 we show the equality of both expressions directly.

We will now introduce some specific assumptions on the switching rates which will suffice to establish a CLT for the Zig-Zag process.

$\inf_{x\geq x_{0}}\lambda(x,+1)>\sup_{x\geq x_{0}}\lambda(x,-1)$, and

$\inf_{x\leq-x_{0}}\lambda(x,-1)>\sup_{x\leq-x_{0}}\lambda(x,+1)$.

In other words, there are constants $M^{-}>m^{-}\geq 0$, $M^{+}>m^{+}\geq 0$, such that

It is established in [3, Theorem 5] that under these conditions the Zig-Zag process is exponentially ergodic.

Suppose Assumption 3.2.1 is satisfied. Let $(Z(t))_{t\geq 0}$ denote the Zig-Zag process with generator (5). Then there exists a unique invariant probability distribution $\mu$ on $E$ for $(Z(t))_{t\geq 0}$. Furthermore there are constants $0<\alpha^{+}\leq M^{+}-m^{+}$ and $0<\alpha^{-}\leq M^{-}-m^{-}$, with $M^{\pm},m^{\pm}$ as above, such that for any function $g\in\mathcal{M}(E)$ satisfying $\mu(g)=0$ and, for $\theta=\pm 1$,

and if $\sigma_{g}^{2}$ as given by (12) satisfies $\sigma_{g}^{2}<\infty$, then

where $B$ denotes a standard Brownian motion and the weak convergence is with respect to the Skorohod topology on $D()$.

Although the constants $\alpha^{\pm}$ are not explicitly specified in the formulation of Theorem 3.9, their construction can be traced in the proof of [3, Theorem 5]. Note that, irrespective of the value of $\alpha^{\pm}$, (13) is satisfied for any sub-exponential function $g$.

Assumption 3.2.1 implies Assumption 2. By Proposition 2.1 it follows that $(Z(t))_{t\geq 0}$ admits a unique invariant probability distribution $\mu$. By tracing the proof of [3, Theorem 5], it follows that there exists a Lyapunov function $V:E\rightarrow[0,\infty)$ such that

for some constants $c^{\pm}>0$ and $\alpha^{\pm}$ as specified in the statement of the theorem, and such that Assumption 3.2 is satisfied with $f:=V$. By the stated assumptions on $g$, possibly after a rescaling by a constant factor, it follows that $|g|\leq f$. By Proposition 3.3, $\mu(f)<\infty$ and there exists a solution $\phi$ for the Poisson equation (8) satisfying $\mu(\phi)=0$ and $|\phi|\leq c_{0}(V+1)$ for some constant $c_{0}>0$. In particular $\mu(|\phi|)<\infty$. The CLT is therefore a result of Theorem 3.7. Under the stronger assumption, $\mu(V^{2})<\infty$ and therefore the FCLT follows by Proposition 3.4.

A sufficient condition for $\sigma_{g}^{2}<\infty$ is that $g\in\mathcal{M}(E)$ and $\lambda:E\rightarrow[0,\infty)$ are of polynomial growth in $x$. Indeed if $g(x,\theta)=O(|x|^{\beta})$ then by Lemma 3.6, for any $\delta>\beta$, $\psi(x)=o(|x|^{\delta})$. Then since $\pi(x)=O(\exp(-(M^{+}-m^{+})x))$ for $x\geq x_{0}$ (and similarly for $x\leq-x_{0}$), it follows that $\psi^{2}(x)\lambda(x,\theta)\pi(x)$ has bounded integral.

2.2 Heavy-tailed distributions

$\lambda:E\rightarrow[0,\infty)$ is continuous. There exist constants $\alpha>0$ and $0\leq\kappa\leq 1$ such that $\lambda(x,+1)\geq\alpha x^{-\kappa}$ for $x>x_{0}$ and $\lambda(x,-1)\geq\alpha(-x)^{-\kappa}$ for $x<-x_{0}$, with $\alpha>2$ in case $\kappa=1$. Furthermore $\lambda(x,-1)=0$ for $x>x_{0}$ and $\lambda(x,+1)=0$ for $x<-x_{0}$.

Suppose Assumption 3.2.2 is satisfied. Let $1\leq\beta<\alpha$ in case $\kappa=1$, and $1\leq\beta<\infty$ in case $\kappa<1$. There exists a norm-like function $V:E\rightarrow[0,\infty)$, and a function $f$ of the form $f(x,\theta)=c|x|^{\beta-1}$ for some $c>0$, and $x_{1}>0$ such that

Let $V$ be given for $x>x_{0}$ by $V(x,+1)=kx^{\beta}$ and $V(x,-1)=\frac{1}{\beta}x^{\beta}$, with

Then for $x>x_{0}$, $LV(x,-1)=-x^{\beta-1}$ and

In the case $\kappa<1$, the negative term will dominate for $x$ sufficiently large. It follows in either case that for a suitable constant $c>0$ and $x_{1}\geq x_{0}$, $LV(x,\pm 1)\leq-cx^{\beta-1}\leq-1$ for all $x\geq x_{1}$. The situation for $x\leq-x_{0}$ is completely analogous, and within $[-x_{0},x_{0}]$, the function $V$ can be continuously and differentiably extended.

In fact for Lemma 3.12 we only require $\alpha>1$ in case $\kappa=1$, because this allows us to choose $\beta\in[1,\alpha)$. However in order to obtain $\mu(V)<\infty$ as required for the proof of the following theorem we need the stronger assumption $\alpha>2$ in case $\kappa=1$.

Suppose Assumption 3.2.2 is satisfied. Let $(Z(t))_{t\geq 0}$ denote the Zig-Zag process with generator (5). Then there exists a unique invariant probability distribution $\mu$ on $E$ for $(Z(t))_{t\geq 0}$. Suppose $g\in\mathcal{M}(E)$ with $\mu(g)=0$ and $g(x,\theta)=O(|x|^{\beta-1})$ where $1\leq\beta<\alpha-1$ in case $\kappa=1$ and $1\leq\beta<\infty$ in case $\kappa<1$. Furthermore suppose $\sigma_{g}^{2}:=4\int_{E}\lambda(x,\theta)\psi^{2}(x)\ d\mu(x,\theta)<\infty$, where $\psi$ is given by (11).

$\kappa=1$, $\alpha>3$ and $1\leq\beta<(\alpha-1)/2$,

where $B$ denotes a standard Brownian motion and the weak convergence is with respect to the Skorohod topology on $D()$.

Assumption 3.2.2 implies Assumption 2 so that by Proposition 3.3 there is a unique invariant probability distribution $\mu$. If $\kappa=1$ in Assumption 3.2.2 then $\frac{d\mu}{dx}(x,\theta)=O(|x_{0}/x|^{\alpha})$. Because $\alpha>2$ we can choose $1\leq\beta<\alpha-1$ in Lemma 3.12, and it follows that the Lyapunov function $V(x,\theta)=O(|x|^{\beta})$ satisfies $\mu(V)<\infty$. If $0\leq\kappa<1$ then $\frac{d\mu}{dx}(x,\theta)=O(\exp(-\alpha/(1-\kappa)|x|^{1-\kappa}))$ and again $\mu(V)<\infty$. The CLT now follows from Theorem 3.7. Under the stronger assumptions, $\mu(V^{2})<\infty$ using the above asymptotic analysis, so that the FCLT follows from Proposition 3.4.

A sufficient condition for $\sigma_{g}^{2}<\infty$ in case $\kappa=1$ is that $\alpha>2$, 1\leq\beta<\min(\alpha-1,\mbox{\frac{1}{2}}\alpha) and $\lambda(x,+1)=O(x^{-1})$. Indeed, in this case there exists a $\delta\in(\beta,\alpha/2)$. Since $\pi(x)=O(|x|^{-\alpha})$ and $\delta<\alpha$ we have that $\pi(x)|x|^{\delta}\rightarrow 0$. Furthermore (10) is satisfied as $g(x)=O(|x|^{\beta-1})$ and $\pi(x)/\pi^{\prime}(x)=O(|x|^{-1})$, so we may deduce from Lemma 3.6 that $\psi(x)=o(|x|^{\delta})$. Hence $\lambda(x)\psi^{2}(x)\pi(x)=o(|x|^{2\delta-1-\alpha})=o(|x|^{-1})$ using that $\delta<\alpha/2$.

2.3 Comparison with Langevin diffusion

Let $A$ denote the generator of the Langevin diffusion with invariant density $\pi(x)=\exp(-U(x))/k$, i.e.

with domain including at least all twice continuously differentiable functions $f$ for which $Af$ is a bounded continuous function.

The proof of this result may be found in Appendix A.4.

In cases where both a CLT holds for the Langevin diffusion and the Zig-Zag process, and the function of interest $g$ does not depend on $\theta$, we can compare the asymptotic variances, given by

where we used (4) to obtain the last equality.

3 Examples

To illustrate the effectiveness of the developed theory we consider several examples. We consider (i) Gaussian distributions, which have light tails and for which the associated Zig-Zag process is exponentially ergodic, and (ii) Student t-distributions, which are heavy tailed so that the associated Zig-Zag process is not exponentially ergodic. For both families of distributions we will consider two types of observables: (a) moments and (b) tail probabilities.

The family of centered one-dimensional Gaussian distributions $\mathcal{N}(0,\nu^{2})$ is described by the potential functions and canonical switching rates

First we consider the asymptotic variance associated with the $k$-th moment for positive integer values of $k$. This corresponds to the mean-zero functional $g(x)=x^{k}-m_{k}$, where

Assumption 3.1 is satisfied for any $k\geq 0$ so that a CLT holds by Theorem 3.1. The asymptotic variance can be computed using (6) to be

The variance of $g$ under $\pi$ is given by

In order to compare the asymptotic variance of the Langevin diffusion, we compute

Expressions for $\psi(x)$ for different values of $k$ are given, along with the computed asymptotic variance for the Zig-Zag process ($\sigma_{g}^{2}$) and Langevin diffusion ($\widetilde{\sigma}_{g}^{2}$), in the following table.

For each of these moments we note that $\sigma^{2}_{g}/\widetilde{\sigma}^{2}_{g}\propto\nu^{-1}$, which suggests that for large variance distributions, the variance of an estimator for $\pi(g)$ using the Zig-Zag process will be considerably lower than that of an estimator generated from a Langevin trajectory.

The result of Theorem 3.1 can be strengthened since by Theorem 3.9 the Functional Central Limit Theorem holds for this entire family of examples.

while the variance of $g$ is given by $\operatorname{Var}_{\pi}(g)=p_{a}(1-p_{a})$.

In Figure 3 we compare the expression (15) with the variance estimated from $10^{5}$ independent simulations of the Zig-Zag process, for different values of $\nu^{2}$.

3.2 Student t-distribution

Consider the family of Student-t distributions with $\nu>0$ degrees of freedom,

and let $\lambda$ denote the canonical switching rates, given by

For integer values of $k$ with $0\leq k<\nu$ we can compute the values of the moments to be

The mean-zero function representing the observable of interest is $g(x)=x^{k}-m_{k}$. Assumption 3.1 is satisfied if $k<(\nu-1)/2$. Moreover we may apply Theorem 3.15 with $\alpha<\nu+1$, $\gamma=1$ and $\beta=k+1$ to see that in the above cases a functional CLT is satisfied under the stated assumption that $k<(\nu-1)/2$.

This may be compared to the Random Walk Metropolis algorithm. In [17, p. 796] it is established that for a finite variance proposal distribution, the range of parameter values for which a CLT holds is $k<\nu/2-1$ which is slightly more restrictive. By tuning the proposal distribution in RWM to have the same decay in the tails, this range can be improved to $k<\nu/2$.

Using (6) we obtain, for the Zig-Zag process,

For $k$ even an also explicit but more cumbersome expression can be obtained from (6).

After evaluating the necessary integrals in (6), we find the asymptotic variance of the Zig-Zag process to be

and, writing ${}_{2}F_{1}$ for the hypergeometric function,

For $\nu=2$, the above expressions simplify to

whereas for other values of $\nu$ the expression for the asymptotic variance can typically not be significantly simplified. See Figure 4 for an experimental verification of these results. We see good agreement with theoretical predictions. Also from Figure 4(b) the rescaled variance of the estimator for $\nu=1$ appears to diverge to infinity as $T\rightarrow\infty$, which suggests that no CLT holds in this case, and thus the condition $\nu>1$ is indeed tight.

Diffusion limit of the Zig-Zag process

In this section we will consider the one dimensional Zig-Zag process with switching rates of the form

for a general non-vanishing space-dependent switching rate $\gamma$. An example arising from applications where $\gamma$ is positive is when Zig-Zag sampling is used in combination with sub-sampling, as discussed in . It is observed in simulations that this gives rise to diffusive behaviour. In this section we show that under an appropriate time change the Zig-Zag process converges weakly to an Itô diffusion, ergodic with respect to $\pi$, with space dependent diffusion coefficient inversely proportional to the switching rate $\gamma$.

We shall focus on behaviour of the Zig-Zag process in the large $\|\gamma\|_{\infty}$ limit. To this end, we shall introduce the rescaling $\gamma^{\epsilon}=\epsilon^{-1}\gamma,$ and denote by $Z^{\epsilon}(t)=(X^{\epsilon}(t),\Theta^{\epsilon}(t))$ the corresponding Zig-Zag process, with generator defined by

where $\lambda^{0}(x,\theta)=\max(0,\theta U^{\prime}(x)).$ Our objective is to prove the following result.

If the process $(\xi(t))_{t\geq 0}$ exists and is non-explosive, then it is ergodic with unique stationary distribution $\pi(x)\propto\exp(-U(x))$.

To prove this result, we will follow an approach similar to that of [13, Theorem 1.5]. The main distinction is that, in [13, Theorem 1.5] the authors introduce a random time-change for the PDMP which produces a limiting SDE with additive noise. On the other hand, the limiting SDE (19) is qualitatively different, in particular it will have multiplicative noise dependent on the switching rate $\gamma$ and moreover is ergodic with respect to the unique stationary disitribution $\pi$. The proof of Theorem 4.1 will be deferred to Section A.5.

We demonstrate the conclusions of Theorem 4.1 using a simple example. Given $U(x)=x^{2}/(2\sigma^{2})$ consider the family of Zig-Zag processes $Z^{\epsilon}(t)=(X^{\epsilon}(t),\Theta^{\epsilon}(t))$ with switching rates

where we choose $\gamma(x)=(1+x^{2})$ for a positive parameter $\epsilon>0$. The resulting process is ergodic, with unique invariant distribution $\pi\sim\mathcal{N}(0,\sigma^{2})$. Applying Theorem 4.1 we know that, in the limit $\epsilon\rightarrow 0$, the time-changed process $X^{\epsilon}(t/\epsilon)$ will converge weakly to an Itô diffusion process $\xi(t)$ given by the unique solution of

It is straightforward to show that $(\xi(t))_{t\geq 0}$ is an ergodic process with unique invariant distribution $\pi$. In Figure 5 we demonstrate this result numerically. Choosing $\sigma^{2}=1$ and for $\epsilon=10,1,0.1$ we plot a histogram of the values of $Z^{\epsilon}(t)$ at values $t/\epsilon=1$,$10$, $20$ and $50$ over $10^{4}$ independent realisations starting from $X^{\epsilon}(0)=2.0$. We compare the result with the corresponding distribution of the diffusion process (21) denoted by the solid line. While for larger values of $\epsilon$ there is a clear discrepancy between $X^{\epsilon}(t)$ and $\xi(t)$, as the speed of the switching rate increases, the Zig-Zag process displays increasing random walk behaviour and shows very good agreement with the diffusion process.

Effective Sample Size for the Zig-Zag process

Provided that a central limit theorem holds, for large $T$, the variance of the estimator $\pi_{T}(f)$ is given to leading order by $T^{-1}\sigma^{2}_{f}$, where $\sigma^{2}_{f}$ is the asymptotic variance for the observable $f$. Suppose we wish to obtain an approximation of $\pi(f)$ within a given error tolerance $\epsilon^{2}$ (in the sense of mean-square error), one can obtain an estimate of the amount of time $T$ that the Zig-Zag process must be simulated, namely

In general, (22) does not reflect the true cost of simulating the Zig-Zag sampler. Indeed, as with all continuous time processes, one can accelerate the mixing of a process simply by introducing a time change $Z^{a}(t)=Z(at)$, for $a>0$. In reality, introducing such a time change will increase the number of switches which occur per unit time, thus increasing the computational effort required to simulate the process up to a given final time $T$.

Assume that $Z(t)$ is simulated using the direct method (see Algorithm 1 in Appendix B). The switching times are determined by a Poisson process with inhomogeneous rate $\int_{0}^{t}\lambda(X(s),\Theta(s))\,ds$. Therefore, the average number of switches occurring in time $[0,T]$ is given by

To quantify the average computational cost of simulating a Zig-Zag sampler we introduce the average switching rate $N_{S}=\lim_{t\rightarrow\infty}t^{-1}N(t)$, which measures the average number of switches occurring per unit time. Since $Z(t)$ is ergodic, then we have that

where we used the explicit formula for $\lambda(x,\theta)$ given in (4). Thus, assuming that $N_{S}$ is finite, after an initial transient period the number of switchings will increase linearly in time with rate $N_{S}$. In terms of computational cost per simulated unit time interval, it is clear that using canonical switching (i.e. $\gamma=0$) is the cheapest option. In this case, the average switching rate will be determined entirely by the target distribution.

For the purpose of comparison with other sampling schemes, it would be ideal to obtain an expression for the variance of the estimator $\frac{1}{T}\int_{0}^{T}f(X_{s})\,ds$ as a function of the number of switches required to simulate the Zig-Zag process up to time $T$. For large $T$ the average number of switches that occurred over $[0,T]$ is approximately $TN_{S}$ where $N_{S}$ is given by (23). Over large time-scales the variance of the estimator $\pi_{T}(f)$ is thus given (for the canonical switching rates, $\gamma=0$), by

where $N(T)$ is the number of switches that occured up to time $T$ and $\psi$ is given by (11).

A useful measure of the effectiveness of a sampling scheme is the effective sample size (ESS), which provides a measurement of the equivalent number of IID draws from $\pi$ which would be required to obtain an estimate for $\pi(f)$ with similar variance. For the Zig-Zag sampler, it is natural to define the ESS as follows

This expression provides a far more natural measure of the effectiveness of the Zig-Zag sampler than e.g. (22). In particular, it is trivial to check that $\mbox{Var}_{\pi}[f]/(\sigma^{2}_{f}N_{S})$ is invariant under time rescaling $t\rightarrow at$, for $a>0$. The use of the number of switches $N(T)$ as a measure of computational cost is also well-justified. One can see from Algorithm 1 that this coincides with the number of evaluations of the gradient of the log target distribution $U(x)$, which in high dimensions, or in the large data regime for Bayesian inference problems (as considered in ) would be the most expensive operation required to compute the next term in the event chain. The ESS is linearly increasing with $N(T)$ by a factor equal to $\operatorname{Var}_{\pi}[f]/(\sigma_{f}^{2}N_{S})$, which determines the efficiency of the Zig-Zag sampler.

Consider the problem of computing moments $x^{k}$ of the Gaussian distribution $\mathcal{N}(0,\nu^{2})$, where $k$ is a natural number. In this case, we can compute the effective switching rate to be $N_{S}=(2\pi\nu^{2})^{-1/2}$, so that using the expression for the asymptotic variance obtain in Example 3.19 we have for $k$ odd

which is independent of $\nu$. A tedious calculation reveals that $ESS>N(T)$, for all such $k$. A similar computation gives, for $k$ even

Evaluating numerically the first few moments using (25) and (26) we obtain

we see that the relation $ESS>N(T)$ appears to hold for general $k$. This demonstrates a non-intuitive phenomenon: the effective sample size of the Zig-Zag process is higher than the number of IID samples. Thus an ergodic average generated from a trajectory of the Zig-Zag process with $N$ switches will tend to have lower variance than a Monte Carlo average of $N$ IID samples of $\pi$. To demonstrate the performance of the Zig-Zag sampler, we generate $10^{5}$ independent realisations of the process ergodic with respect to $\mathcal{N}(0,4)$, and in Figure 6 plot the variance for estimators of the first two moments, as a function of $N$ (the maximum number of switches). We also plot the variance for a MC average generated from IID samples, as well as for a Random Walk Metropolis-Hastings (RWMH) scheme with manually tuned step-size. We see that even after manually tuning the step-size of the RWMH chain, the asymptotic variance of the corresponding estimator is still an order of magnitude higher that that of the IID chain and Zig-Zag sampler. In both cases, the ratio of variances for the Zig-Zag sampler and IID average is constant, independent of $N$, as predicted by (25) and (26).

The fact that the Zig-Zag sampler is able to achieve effective sample sizes which beat IID is a property which is closely tied to the non-reversible nature of the Zig-Zag process. While we have demonstrated this property for the Gaussian case, one should not interpret this as a general result. Indeed, in the following example we repeat the above experiment for the Student t-distribution, and we show that although the Zig-Zag sampler outperforms the corresponding RWMH chain, it will not have ESS higher than that of an IID chain.

Following Example 3.21, we consider once again the problem of the first moment of the Student t-distribution with $\nu$ degrees of freedom. In Figure 7 we plot the variance of estimates for the first moment obtained from the Zig-Zag process using canonical switching rate (37), for $\nu=4$, $6$ and $8$. Each point is generated from $M=10^{5}$ independent realisations of the process. Note that for the observable $f(x)=x$, Assumption 3.1 holds for each value of $\nu$. As in the previous example, we also plot the variance of a Monte-Carlo estimator generated from $M$ IID samples, as well a from a manually tuned RWMH chain.

In this case the effective sample size of the Zig-Zag sampler will not be higher than that of the IID estimator, in general. However, as the degrees of freedom $\nu$ goes to infinity, the target distribution becomes increasingly Gaussian, and for sufficiently large $\nu$, the Zig-Zag sampler will exhibit lower variance than the corresponding IID scheme.

Appendix A

Because $\lambda$ is locally bounded, [7, Assumption 3.1] is satisfied, and a piecewise deterministic Markov process can be constructed as described in . Then, by [7, Theorem 5.5], $L$ is the extended generator. The Feller property is established by tracing the proof of [3, Proposition 4], for which only continuity of $\lambda$ is required. Since $\lambda$ is continuous and because $\lambda(x,\theta)>0$ for $\theta x\geq x_{0}$, we have in fact that, for any $x_{1}>x_{0}$, there exists a $c>0$ such that

The proof that compact sets are petite is now a straightforward adaptation of the proof of [3, Lemma 15], and a Markov process with this property is $\varphi$-irreducible; in particular there exists at most a single invariant measure. The stationarity of $\mu$ is established in [3, Proposition 5].

A.2 Technical results towards the CLT

The following lemma is a continuous time variant of [10, Exercise 2.4.6].

Let $\varepsilon>0$ and $\gamma>0$. Let $\beta=\varepsilon\gamma^{2}/(2\sigma^{2})$. Pick $T>0$ such that for all $t\geq T$, $|N(t)/a(t)-1|>\beta$ with probability at most $\varepsilon/2$. For fixed $t\geq T$, let $\Omega(t)$ denote the event in which $|N(t)/a(t)-1|\leq\beta$. On $\Omega(t)$, $|N(t)-a(t)|\leq\lfloor\beta a(t)\rfloor\leq\beta a(t)$. By Kolmogorov’s maximal inequality,

This establishes that $\frac{1}{\sqrt{a(t)}}\left(\sum_{i=1}^{N}(t)Y_{i}-\sum_{i=1}^{a(t)}Y_{i}\right)$ converges in probability to 0. The stated result now follows from the classical central limit theorem applied to $\frac{1}{\sqrt{a(t)}}\sum_{i=1}^{a(t)}Y_{i}$.

Since $\phi\in\mathcal{D}(L)$ it follows that $M$ is a local martingale. Due to stationarity we have

where we used that $|g|\leq f$ and $\mu(f)<\infty$ by Proposition 3.3. It follows that $M$ is a martingale. We have

where \psi(x)=\mbox{\frac{1}{2}}(\phi(x,+1)-\phi(x,-1)). Using [18, Theorem 26.6 (vii), (viii)] the quadratic variation of $M$ and predictable quadratic variation are given by the stated expressions.

Assume without loss of generality that $\mu(g)=0$. Writing out the relation $L\phi(x,\theta)=-g(x,\theta)$ for $\theta=\pm 1$ and adding the two equations gives

It remains to verify that the constant $c$ vanishes. By Proposition 3.3, we have $|\phi|\leq c_{0}(V+1)$ and hence

By the assumption that $\pi(x)V(x,\pm 1)\rightarrow 0$, it therefore follows that $\pi(x)\psi(x)\rightarrow 0$ as $|x|\rightarrow\infty$. Multiplying (27) by $\pi$, we have that

A.3 Equivalence of expressions for asymptotic variance

A natural question to ask is whether the two expressions for asymptotic variance, given by (6) and (12) are equivalent in cases where both expressions are valid. Suppose for an observable $g$ such that $\pi(g)=0$,

Assuming that (28) and (29) hold, and that the potential $U$ satisfies $U(0)=0$, then we can show that both expressions are equal. Considering the term

where we use (29) to eliminate the contribution due to the upper integration limit. Similarly, we have

for which the second term is zero, by (28). Exchanging the integrals we obtain

Combining (30) and (31) it follows immediately that the expressions for asymptotic variance respectively given by (6) and (12) are equal.

A.4 Proof of Proposition 3.18

Write $P^{s}$ for the Markov semigroup corresponding to the Langevin diffusion, with generator $A$. By [19, Corollary 1.9], a CLT is satisfied if there exists a constant $c>0$ such that

which is satisfied for $c=\|\psi\|_{L^{2}(\pi)}$. In this case, by [19, Corollary 1.9], the asymptotic variance admits the expression

where $\varphi$ satisfies the Poisson equation $A\varphi=-g$. By the Poisson equation for $\varphi$,

By a similar argument as in the proof of Lemma 3.6, using that $\varphi\in\mathcal{D}(A)$ and hence $\varphi^{\prime}\in L^{2}(\pi)$, it follows that $c=0$ and hence $\varphi^{\prime}(x)=-\psi(x)$.

We now prove the converse. To this end, suppose that

where the equality holds due to [5, Lemma 2.3]. For any $t>0$ define

Note that $g_{t}\in\mathcal{D}(A)$ and satisfies

We follow the approach of [5, Theorem 3.3]. Below, let $f^{\prime}$ denote $\frac{d}{dx}f$. Given $s\leq t$,

It follows that the family $(g^{\prime}_{t})_{t>0}$ is Cauchy in $L^{2}(\pi)$, so that it strongly converges to a limit $-\eta\in L^{2}(\pi)$. The weak formulation of (33) is given by

We have $\lim_{t\rightarrow\infty}P^{t}g=\pi(g)=0$, so that by dominated convergence $\langle P^{t}g,v\rangle_{L^{2}(\pi)}\rightarrow 0$ as $t\rightarrow\infty$, and thus taking the $t\rightarrow\infty$ limit in (34) gives

A.5 Proof of Theorem 4.1

In this section we prove Theorem 4.1, following the approach of . To this end, consider the function

is a remainder term which is measurable and independent of $\epsilon$. Defining

it follows (using that $f$ is in the domain of the extended generator, see [7, Theorem 5.5]), that

is a local martingale with respect to the filtration $\mathcal{F}^{\epsilon}_{t}$ generated by $\{Z^{\epsilon}(t):t\in[0,T]\}.$ Similarly, applying the generator to $g(x,\theta):=f^{2}(x,\theta)$, we obtain

where $b(x)$ is as above, $a(x)=\frac{1}{\gamma(x)}$, and $R_{2}(x,\theta)$ can be written as $R_{2}=R^{(1)}_{2}+\epsilon R^{(2)}_{2}+\epsilon^{2}R^{(3)}_{2}$, where the terms

are measurable and independent of $\epsilon$. We thus obtain that

is a local martingale with respect to the filtration $\mathcal{F}^{\epsilon}_{t}$. We now decompose the square local martingale $(M^{\epsilon}(t))^{2}$ into a local martingale term and a remainder. To this end, defining $J^{\epsilon}(t)=\int_{0}^{t}j^{\epsilon}(s)\,ds,$ use integration by parts to obtain

where the terms of order $\epsilon^{2}$ or higher are collected in the remainder term $R_{3}(x,\theta)$. It follows that

is a local martingale with respect to $\mathcal{F}^{\epsilon}_{t}.$ Applying the time change $t\rightarrow t/\epsilon$ we see that

are local martingales with respect to the filtration $\mathcal{F}^{\epsilon}_{t}:=\mathcal{F}_{t/\epsilon}$, $t\geq 0$. We now verify the conditions of [11, Theorem VII.4.1] to derive the diffusive limit. To this end, define

as well as the stopping time $\tau_{R}^{\epsilon}:=\inf\left\{t\geq 0\,:\,\left|X^{\epsilon}(t/\epsilon)\right|\geq R\mbox{ or }\left|X^{\epsilon}(t/\epsilon-)\right|\geq R\right\}$. From our assumptions we have that, for each $R\geq 0:$

and thus converges to $0$ almost surely as $\epsilon\rightarrow 0$. Similarly

almost surely as $\epsilon\rightarrow 0$. Finally, noting that $X^{\epsilon}(t)$ is continuous, we have that

Since the well-posedness of this martingale problem is equivalent to the existence and uniqueness of a weak solution $(\xi(t))_{t\geq 0}$ for (19), the proof is complete.

Appendix B Simulation of the Zig-Zag process

In this section we describe some computational methods for simulating the process $Z(t)=(X(t),\Theta(t))$ and use results from previous sections in analyzing these methods. As with the rest of this paper, we shall focus in particular on the one-dimensional case, referring the reader to for specifics of the general case.

Given the state $(X(T_{0}),\Theta(T_{0}))=(x_{0},\theta_{0})$ at switching time $T_{0}$, the next random switching time is given by $T_{1}=T_{0}+\tau$ where $\tau$ satisfies

In the case where $G(t)=\int_{0}^{t}\lambda(x_{0}+s\theta_{0},\theta_{0})\,ds$ has an explictly computable generalised inverse

then applying an inverse transformation, the random variable $\tau=H(-\log u)$, $u\sim U$ satisfies (35). An algorithm for simulating $Z(t)$ based on this approach is detailed in Algorithm 1.

The computational cost of Algorithm 1 clearly depends on the switching intensity, i.e. a Zig-Zag sampler with a higher switching intensity will require more computational cost to be simulated up to a fixed time $T$. Indeed, while the Zig-Zag sampler does not reject samples like a Metropolis-Hastings scheme, frequent switching will cause the process $Z(t)$ to mix slowly.

A straightforward calculation shows that, given $(x,\theta)\in E$ the generalised inverse of $G(t)=\nu^{-2}\int_{0}^{t}\max(0,x_{0}+s\theta)\,ds$ can be written explicitly as

for $z>0$. In this case, the average switching rate is then given by $N_{S}=(2\pi\nu^{2})^{-1/2}$.

It is also possible to sample from a Student t-distribution with $\nu$ degrees of freedom, i.e.

using the direct Zig-Zag sampling approach. For this distribution, the canonical switching function is given by

Given $(x,\theta)\in E$, the generalised inverse of $G(t)=\int_{0}^{t}\lambda(x+\theta s,\theta)\,ds$ can be written as

The average switching rate is equal to the normalization constant for (36), i.e. $N_{S}=\frac{\Gamma((\nu+1)/2)}{\sqrt{\nu\pi}\Gamma(\nu/2)}$. The resulting process will be ergodic with respect to the target distribution $\pi$, for all $\nu>0$. Conditions under which a central limit theorem holds will be studied in Section 3.3.

B.2 Sampling with Poisson Thinning

In general we will not be able to compute the generalized inverse of $G$ explictly. In many cases however, it is possible to obtain an upper bound $\Lambda(t;x,\theta_{0})$ such that $m(t):=\lambda(x_{0}+\theta_{0}t,\theta_{0})\leq\Lambda(t;x_{0},\theta_{0})$, for all $t\geq 0$, $(x_{0},\theta_{0})\in E$, and where $\Lambda(t)$ has an explicitly computable inverse $\widetilde{H}$. In this case, one can simulate the random switching times using a standard Poisson thinning approach . Using the upper bound $\Lambda(t;x_{0},\theta_{0})$ a candidate switching time $T_{1}=T_{0}+\tau$ is generated, such that

A switch (i.e. $\Theta_{1}=-\Theta_{0})$ will occur at $T_{1}$ with probability $m(t_{1})/\Lambda(t_{1};x_{0},\theta_{0})$. An algorithm for sampling $Z(t)$ based on this approach is detailed in Algorithm 2.

Identifying such a computable upper bound is highly problem specific, however we can highlight two frequently arising scenarios where upper bounds can be easily constructed.

For fixed $(x,\theta)\in E$, the integrated intensity function $G(t)=\int_{0}^{t}\Lambda(s;x,\theta)\,ds$ has generalised inverse

This case arises naturally in various Bayesian inference problems, in particular logistic regression, see [2, Section 6.5].

The number of switches that occur in a given time interval will depend on the intensity function $\Lambda(t;x,\theta)$, and clearly, a poor choice of this upper bound will cause Algorithm 2 to undergo many potential switch events which are rejected. In particular, if the process $Z(t)$ is in stationarity, then the average switching rate will always be higher or equal to that of the direct scheme described in Algorithm 1.

B.3 Computing ergodic averages

While the event chain $(X(T_{k}),\Theta(T_{k}))_{k=0}^{\infty}$ defines a Markov chain, it will not be ergodic with respect to the target distribution $\pi$. To compute an ergodic average for a given observable $f$, the entire continuous time realisation must be used as follows

Since the Zig-Zag process moves linearly between switches, this can be decomposed into a sum of integrals over straight lines. Indeed, for $T=T_{K}$, for some $K$ we have

where $\tau_{k}=T_{k+1}-T_{k}$. In many cases, the integral in (38) can be computed exactly. For example, first and $p^{th}$ moment can be computed ergodically via the expressions

respectively. For more complicated observables it will not be possible to evaluate (38) analytically, and one must resort to some form of quadrature scheme, for example Euler or other higher order methods.

The authors acknowledge the EPSRC for support under grants EP/D002060/1, EP/K014463/1 (Joris Bierkens), EP/J009636/1, EP/L020564/1 (Andrew Duncan) as well as EP/K009788/2 (both authors). Furthermore we acknowledge support of the Lloyds Registry Foundation through the Alan Turing Institute. We are grateful to the referee and associate editor for useful suggestions with regards to the mathematical exposition, which have certainly helped to improve this paper.

References