Piecewise Deterministic Markov Processes for Continuous-Time Monte
Carlo
Paul Fearnhead1,† , Joris Bierkens2 , Murray Pollock2 and Gareth O Roberts2
1

Department of Mathematics and Statistics, Lancaster University
2

Department of Statistics, University of Warwick

arXiv:1611.07873v1 [stat.CO] 23 Nov 2016

†

Correspondence: p.fearnhead@lancaster.ac.uk
November 24, 2016

Abstract
Recently there have been exciting developments in Monte Carlo methods, with the development of new MCMC and
sequential Monte Carlo (SMC) algorithms which are based on continuous-time, rather than discrete-time, Markov
processes. This has led to some fundamentally new Monte Carlo algorithms which can be used to sample from,
say, a posterior distribution. Interestingly, continuous-time algorithms seem particularly well suited to Bayesian
analysis in big-data settings as they need only access a small sub-set of data points at each iteration, and yet are still
guaranteed to target the true posterior distribution. Whilst continuous-time MCMC and SMC methods have been
developed independently we show here that they are related by the fact that both involve simulating a piecewise
deterministic Markov process. Furthermore we show that the methods developed to date are just specific cases of
a potentially much wider class of continuous-time Monte Carlo algorithms. We give an informal introduction to
piecewise deterministic Markov processes, covering the aspects relevant to these new Monte Carlo algorithms, with
a view to making the development of new continuous-time Monte Carlo more accessible. We focus on how and why
sub-sampling ideas can be used with these algorithms, and aim to give insight into how these new algorithms can be
implemented, and what are some of the issues that affect their efficiency.
Keywords: Bayesian Statistics, Big data, Bouncy particle sampler, Continuous-time importance sampling, Control
variates, SCALE, Zig-zag sampler.

1

Introduction

Monte Carlo methods, such as MCMC and SMC, have been central to the application of Bayesian statistics to
real-world problems (Robert and Casella, 2011; McGrayne, 2011). These established Monte Carlo methods are based
upon simulating discrete-time Markov processes. For example MCMC algorithms simulate a discrete-time Markov
chain constructed to have a target distribution of interest, the posterior distribution in Bayesian inference, as its
stationary distribution. Whilst SMC methods involve propagating and re-weighting particles so that a final set of
weighted particles approximate a target distribution. The propagation step here also involves simulating from a
discrete-time Markov chain.
1

In the past few years there have been exciting developments in MCMC and SMC methods based on continuoustime versions of these Monte Carlo methods. For example, continuous-time MCMC algorithms have been proposed
(Peters and de With, 2012; Bouchard-Côté et al., 2015; Bierkens and Roberts, 2015; Bierkens et al., 2016) that involve
simulating a continuous-time Markov process that has been designed to have a target distribution of interest as its
stationary distribution. These continuous-time MCMC algorithms were originally motivated as they are examples of
non-reversible Markov processes. There is substantial evidence that non-reversible MCMC algorithms will be more
efficient than standard MCMC algorithms that are reversible (Neal, 1998; Diaconis et al., 2000; Neal, 2004; Bierkens,
2015), and there is empirical evidence that these continuous-time MCMC algorithms are more efficient than their
discrete-time counterparts (see e.g. Bouchard-Côté et al., 2015). Similarly a continuous-time version of SMC has
also been recently proposed (Fearnhead et al., 2016), which involves propagating particles using a continuous-time
Markov process. The original motivation for this was to be able to target distributions related to infinite-dimensional
stochastic processes, such as diffusions, without resorting to any time-discretisation approximations. However, we
show below that one application of this methodology is to generate weighted-samples from a target distribution of
interest, giving an alternative interpretation of the recently proposed SCALE algorithm of Pollock et al. (2016a).
The purpose of this paper is to show that continuous-time MCMC and continuous-time SMC methods are linked
through the fact that they both are based upon simulating continuous-time processes called piecewise-deterministic
Markov processes. These are processes that evolve deterministically between a countable set of random event times.
The stochasticity in the process is due to the randomness in when these events occur, and possibly random dynamics
at the event times. These processes are natural building blocks of continuous-time Monte Carlo methods, as they
involve a finite amount of computation to simulate as only a finite number of events and transitions are simulated
in any fixed time-interval.
Furthermore we aim to show that the methods that have been developed to date are just specific examples of a much
wider class of potential continuous-time MCMC and SMC methods that are based on piecewise-deterministic Markov
processes. By giving an informal introduction to theory of piecewise-deterministic Markov processes, with emphasis
on aspects most relevant to the development of valid Monte Carlo methods, we hope to make the development of
new continuous-time Monte Carlo methods more accessible and help stimulate work in this area.
One aspect of continuous-time Monte Carlo that is particularly relevant for modern applications of Bayesian statistics
is that they seem well-suited to big-data problems. If we denote our target distribution by π(x) then the dynamics
of these methods depend on the target through ∇ log π(x). Now if π(x) is a posterior distribution, then it will often
be in product-form, where each factor relates to a data-point or set of data-points. This means that ∇ log π(x) is
a sum, and hence is easy to approximate unbiasedly using sub-samples of the data. It turns out we can use these
unbiased estimators within the continuous-time Monte Carlo methods without effecting their validity. That is, the
algorithms will still target π(x). This is in comparison to other discrete-time MCMC methods that use sub-sampling
(Welling and Teh, 2011; Bardenet et al., 2015; Ma et al., 2015; Quiroz et al., 2015), where the approximation in
the sub-sampling means that the algorithms will only target an approximation to π(x). It also compares favourably
with big-data methods that independently run MCMC on batches of data, and then combines the MCMC samples
in some way (Neiswanger et al., 2013; Scott et al., 2016; Srivastava et al., 2015). As the combination steps involved
will also introduce some approximation error.
The outline of the paper is as follows. In the next section we give an informal introduction Piecewise Deterministic
Markov processes. Our aim is to cover key relevant concepts linked to these processes whilst avoiding technical
details. Those interested in a more rigorous introduction should see Davis (1984) and Davis (1993). Sections 3 and 4
then cover continuous-time versions of SMC and MCMC respectively. These have been written so that either section

2

could be read independently of the other. Our aim for each section is to introduce the continuous-time Monte Carlo
algorithm, show how it relates to a piecewise deterministic Markov process, and how we can use the theory for these
processes to see that the Monte Carlo algorithms target the correct distribution. We also cover how these algorithms
can be implemented using sub-sampling ideas, and highlight the importance of low-variance sub-sampling estimators
for obtaining highly efficient samplers for big-data.

2

Piecewise Deterministic Markov Processes

The continuous-time versions of SMC, or sequential importance sampling, and MCMC that we will consider later are
all examples of time-homogeneous piecewise-deterministic Markov processes. We will henceforth call these piecewise
deterministic processes or PDPs.
A PDP is a continuous-time stochastic process. Throughout we will denote the state of a PDP at time t by Zt .
The dynamics of the PDP involves random events, with deterministic dynamics between events and possibly random
transitions at events. These dynamics are thus defined through specifying three quantities:
(i) The deterministic dynamics. We will assume that these are specified through an ordinary differential
equation:
(i)

dzt
= Φi (zt ),
dt

(1)

for i = 1, . . . , d, for some known vector-valued function Φ = (Φ1 (z), . . . , Φd (z)). This will lead to a deterministic
transition function, so that the solution of the differential equation starting from value zt and run for a time
interval of length s will give
zs+t = Ψ(zt , s)
for some function Ψ.
(ii) The event rate. Events will occur singularly at a rate, λ(zt ), that depends only on the current position of the
process. The probability of an event in interval in [t, t + h] given the state at time t, zt , is thus λ(zt )h + o(h).
(iii) The transition distribution at events. At each event the state of the process will change, according to
some transition kernel. For an event at time τ , if zτ − denotes the state immediately prior to the event, then
the state at time τ will be drawn from a distribution with density q(·|zτ − ).
To define a PDP process we will also need to specify its initial condition. We will assume that Z0 is drawn from
some known distribution with density function p0 (·).

2.1

Simulating a PDP

To be able to use a PDP as the basis of an importance sampling or MCMC algorithm, we will need to be able to
simulate from it. A general approach to simulating a PDP is to iterate the following steps:
(S1) Given the current time t and state of the PDP zt , simulate the next event time, τ say.
(S2) Calculate the value of the process immediately before the next event time
zτ − = Ψ(zt , τ − t).

3

(S3) Simulate the new value of the process, immediately after the event, from q(zτ |zτ − ).
The simulation algorithm is initiated with a current time t = 0 and with Z0 drawn from the initial distribution of
the process. To simulate the process for a time interval T these steps can be iterated until the first event time after
T . If we wish to then simulate the value of the process at a time, s, between events we just find the event time, τ ,
immediately prior to s; the value of the process immediately after the event, zτ ; and then set
zs = Ψ(zτ , s − τ ).
If s is a time before the first event, then we would use τ = 0.
Below we will assume that our PDP has been chosen so that Ψ(·, ·) is known analytically and that the proposal
distribution at events, q(·|·) can be easily simulated from. Thus the only challenging step to simulating a PDP will
be simulating the next event in step S1. This involves simulating the next event in a time-inhomogeneous Poisson
process.
The first thing to note is that the event rate in (S1) can be written as a deterministic function of time, as the state
dynamics are deterministic until the next event. If we are currently at time t with state zt , then for any future time
t + s before the next event, the state will be zt+s = Ψ(zt , s). Thus the event rate will be
λ(zt+s ) = λ(Ψ(zt , s)) = λ̃zt (s),
for a suitably defined function λ̃zt (·). We can simulate the time until the next event, s, as the time of the first event
in a Poisson process of rate λ̃zt (u).
If the event rate is a simple function, then we can simulate events directly. If we define Λz (s) =

Rs
0

λ̃z (u)du, then we

can simulate such an event by (i) simulating u, the realisation of an exponential random with rate 1, and (ii) finding
s > 0 the solution of Λz (s) = u (e.g. Cinlar, 2013).
For more complicated rate functions either calculating Λz (s) or solving the equation in step (ii) may not be tractable.
In such cases the most general approach to simulating event times is by thinning, or adaptive thinning (e.g. Lewis
and Shedler, 1979).
If we can upper-bound the event rate, λ̃zt (u) < λ+ (u), then thinning works by simulating possible events of a Poisson
process of rate λ+ (u) and accepting a possible event at time u as an actual event with probability λ̃zt (u)/λ+ (u).
The time of the first accepted event will be the time until the next event for our PDP. This requires the upper bound
λ+ (u) to be such that simulating events from a Poission process of rate λ+ (u) is straightforward – for example λ+ (u)
is constant or linear in u, or piecewise constant or piecewise linear. Obviously the lower the bound λ+ , the more
computationally efficient this approach will be.

2.2

Analysing a PDP

We now give informal introductions to some of the mathematical tools for analysing a PDP. These are introduced as
they are used later to show that the continuous-time Monte Carlo methods we review have appropriate properties. For
example, we introduce the generator of a PDP in the following section, and this is used to show that the continuoustime importance samplers of Section 3 produce properly weighted samples (Liu and Chen, 1998). We then introduce
the Fokker-Planck equation for a PDP, which can be used to showed that the continuous-time MCMC methods of
Section 4 have the correct invariant distribution. Understanding both the generator of a PDP and its Fokker-Planck
equation is important if we wish to develop new versions of these continuous-time Monte Carlo methods.
4

2.2.1

The Generator

The generator of a continuous-time, time-homogeneous, Markov process is an operator that acts on functions of the
state variable. We will denote the generator by A. For suitable functions f (z), the generator is defined by
Af (z) = lim

δt→0

E(f (Zt+δt )|Zt = z) − f (z)
.
δt

The fact that the process is time-homogeneous means that the right-hand side does not depend on t. We can
interpret the generator applied to a function f (z), as giving the derivative of the expectation of f (Zt ) conditional on
the current value of Zt . This quantity is only well-defined for certain functions, and the set of such functions is called
the domain of the generator. As the generator specifies how the expectation of any suitable function f (·) changes
over time, it uniquely defines the dynamics of the underlying continuous-time stochastic process, in a similar way
that knowing the moment generating function of a random variable will uniquely determine its distribution (Ethier
and Kurtz, 2005).
If we are interested in the derivative of the expectation of a function of our PDP at a time t, then we can write this
as
dE(f (Zt ))
E(f (Zt+δt ) − E(f (Zt ))
= lim
= lim Et
δt→0
δt→0
dt
δt



Et+δt|t (f (Zt+δt )|Zt ) − f (Zt )
δt


,

where in the last expression the inner expectation is with respect to Zt+δt given Zt and the outer expectation with
respect to Zt . Assuming we can exchange the outer expectation and the limit we get
dE(Zt )
= Et (Af (Zt )) .
dt

(2)

Thus the derivative of the expectation of our function is the expectation of the generator applied to the function.
Davis (1984) gives the generator for a piecewise deterministic process:
Z
Af (z) = Φ(z) · ∇f (z) + λ(z) q(z0 |z)[f (z0 ) − f (z)]dz0 ,

(3)

for functions f (·) such that t 7→ f (Ψ(z, t)) is absolutely continuous. The form of the generator has a simple
interpretation. The first term on the right-hand side relates to the deterministic dynamics. For deterministic
dynamics the generator is just the time-derivative of f (zt ), which by the product rule is
d

(i)

df (zt ) X ∂f (zt ) ∂zt
=
= Φ(z) · ∇f (z),
(i)
dzt
∂t
i=1 ∂zt
where Φ(z) is define in (1). The second term on the right-hand side is the change in expectation at events. The
probability of an event in time [t, t + h] is λ(zt )h + o(h) and the change in expectation conditional on event occuring,
up to terms of o(h), is given by the integral on the right-hand side.

2.2.2

The Forward Operator and Fokker-Planck Equation

We can define the adjoint of a generator of a continuous-time Markov process, A∗ , such that for suitable functions
g(z) and f (z)
Z

Z
g(z)Af (z)dz =

f (z)A∗ g(z)dz.

For function g(z) and f (z) which are both continuously differentiable and have compact support, it is straightforward
to show that the adjoint of the generator of a PDP (3) is
∗

A g(z) = −

d
X
∂Φi (z)g(z)
i=1

∂z(i)

Z
+

5

g(z0 )λ(z0 )q(z|z0 )dz0 − g(z)λ(z)

(4)

The first term equates to the adjoint of the first term of the generator, and is obtained using a straightforward
integration by parts argument. To obtain this formula for more general functions f and g requires more careful
arguments to justify the change of order of integrations.
Now if we define the density function of our PDP at time t to be pt (z) then from (2) we get that, for suitable function
f , the derivate of the expectation of f (Zt ) with respect to t is
Z
Z
dE(Zt )
= Et (Af (Zt )) = pt (z)Af (z)dz = f (z)A∗ pt (z)dz.
dt
However we can equally write this derivative as
d
dE(Zt )
=
dt
dt

Z

Z
pt (z)f (z)dz =

∂pt (z)
f (z)dz,
∂t

assuming we can interchange differentiation and integration. This gives that
Z
Z
dpt (z)
f (z)dz = f (z)A∗ pt (z)dz.
dt
As this holds for sufficiently many functions f (z) we get
∂pt (z)
= A∗ pt (z).
∂t
This is a partial differential equation for the distribution of the stochastic process, known as the Fokker-Planck or
Forward Kolmogorov equation. If p(z) is an invariant distribution of our PDP than it will satisfy A∗ p(z) = 0, which
gives
−

d
X
∂Φi (z)p(z)
i=1

∂z(i)

Z
+

p(z0 )λ(z0 )q(z|z0 )dz0 − p(z)λ(z) = 0.

The first term here relates to the change in probability mass caused by the deterministic dynamics, the second term
relates to the probability flow into state z and the final term the probability flow out of state z. For an invariant
distribution these will cancel for all states z.

3

Continuous-Time Sequential Importance Sampling

A continuous-time version of sequential importance sampling was first developed to solve the problem of simulating
from a diffusion. In this situation we have a diffusion process, Xt , defined as the solution to an SDE
dXt = b(Xt )dt + σ(Xt )dBt ,
where b(x) is the d-dimensional drift, Bt is d dimensional Brownian motion, and σ(x) is a d by d matrix that defines
the instantaneous variance of the process. We have an initial distribution p0 (x) for the diffusion, and wish to sample
from the distribution of the process at some future time or times. If we denote the density of this process at time t,
by pt (x), then the challenge is to sample from pt (x) for diffusion processes where we cannot write down what pt (x)
is.
The Exact Algorithm (Beskos and Roberts, 2005; Beskos et al., 2006) and its variants (Beskos et al., 2008; Pollock
et al., 2016b) have given a number of algorithms for simulating from such a diffusion process, but only under strong
conditions on the drift and instantaneous variance. For example it is commonly required that σ(x) is constant, and
that the drift can be written as the gradient of some potential function. Almost all uni-variate diffusion processes
can be transformed to satisfy these requirements, but few multivariate diffusion processes can.

6

Whilst we do not know what pt (x) is for any t > 0, we do know that it solves the Fokker-Planck equation for the
diffusion
d
d
d
X
∂pt (x)
∂bi (x)pt (x) 1 X X ∂ 2 Σij (x)pt (x)
=−
+
,
∂t
∂xi
2 i=1 j=1
∂xi xj
i=1

where Σ = σ T σ. This motivates the following question: can we use our knowledge of the Fokker-Planck equation for
the process of interest in order to develop a valid importance sampling algorithm to sample from pt (x)?
The continuous-time importance sampling (CIS) procedure of Fearnhead et al. (2016), which we describe below, will
in fact enable us to do so. We will present it in a slightly more general form, in that we will use CIS to sample from
a distribution pt (x) that is the solution to a partial differential equation
∂pt (x)
= L∗ pt (x)
∂t
for some known operator L∗ and subject to a known initial condition p0 (x). This then allows for sampling from more
general continuous-time Markov processes, where L∗ would be the adjoint of the generator for that process.

3.1

The CIS Algorithm

The idea of CIS is similar in spirit to a standard importance sampler. We will choose a tractable proposal process, for
the problem of sampling from a diffusion this is most naturally chosen to be Brownian motion. This proposal process
must have a know transition density which is simple to sample from. We will simulate paths from this proposal
process up to time t, and then construct an importance sampling weight. The challenge is that we need to calculate
an importance sampling weight without knowing pt (x). The property we want from our importance sampler is that
if we simulate a value and weight at time t, (X, W ), then for suitable functions f (x) we will have that W f (X) will
be an unbiased estimator of the expectation of f (Xt ) for our target process,
Z
E(f (Xt )) = f (x)pt (x)dx.
The original CIS algorithm was derived by taking a limit of a discrete-time sequential importance sampler. Below
we give the general form of the resulting CIS process. The key to making this a valid importance sampler is choosing
the incremental weight (see step 2 below) appropriately. We will show how we can derive the form of the incremental
weight by viewing the CIS process as a PDP, and using the generator of the PDP to calculate expectations with
respect to the CIS process.
Denote the transition density of the proposal process over time interval s as qs (x|y), and assume we have chosen an
event-rate λ̃(s) > 0 for s > 0. The CIS algorithm is of the form:
(0) Set τ = 0, W0 = 1 and simulate X0 from the initial distribution of the target process, p0 (x).
(1) Simulate a new event time τ 0 > τ , with the inter-event time s = τ 0 − τ being drawn from a Poisson process
with rate λ̃(s).
(2) If τ 0 > t then simulate Xt from qt−τ (·|Xτ 0 ), and set Wt = Wτ 0 . Otherwise simulate Xτ 0 from qs (·|Xτ ), update
the weight as


ρ(Xτ 0 , Xτ , s)
,
Wτ 0 = Wτ 1 +
λ̃(s)
and set τ = τ 0 and go to (1).

7

(5)

Step (0) is just an initialisation step. The idea of the algorithm is that we use random event times, simulated in
step (1), at which we evaluate the proposal process. Based on the value of the process at both this event and the
preceeding event we then update the importance sampling weight. As is standard in sequential importance sampling,
the new weight (calculated in step 2) is just the old weight multiplied by an incremental weight. Figure 1 gives an
example of the output of this algorithm.
The key to making this a valid importance sampling algorithm is working out an appropriate form of the incremental
weight. Without loss of generality we have written the incremental weight in the form given inside the square brackets
in (5), as this will simplify the derivation later. To specify the incremental weight, we need to appropriately choose
the function ρ.

3.1.1

CIS as a PDP

It can be seen that the CIS process is just a piecewise deterministic Markov process. The only randomness in the
process is the timing of the events and the transitions at the events. We can formalise this by defining a PDP with
state Zt = (Yt , Wt , St ) where Yt is the value of the CIS process simulated at the most recent event prior to t, Wt
is the importance sampling weight at time t and St is the time since the most recent event. See Figure 1 for an
example.
This PDP has deterministic dynamics given by differential equations
dYt
= 0,
dt

dWt
= 0,
dt

and

dSt
= 1.
dt

Events occur at a rate λ(Zt ) = λ̃(St ). At an event at time τ we simulate Yτ from qSτ − (·|Yτ − ), set Sτ = 0 and
update Wτ as described in (5).
The generator of Zt , given by (3), is
∂h(y, w, s)
Ah(y, w, s) =
+ λ̃(s)
∂s


 

Z  
ρ(y0 , y, s)
0
h y ,w 1 +
, 0 − h(y, w, s) qs (y0 | y) dy0 .
λ̃(s)

(6)

Given the state at time t, Zt = (Yt , Wt , St ), our estimate of the expectation of f (Xt ) under our target process will
be Wt f (X), where X is drawn from qSt (· | Yt ). The requirement that our algorithm is a valid importance sampler
then becomes that




Z

EX,Z (Wt f (X)) := EZ Wt

f (x)qSt (x|Yt )dx

Z
=

f (x)pt (x)dt.

Here the first expectation is with respect to Z, the PDP, and X the simulated value for Xt , and the second expectation
is just with respect to Z.
As this must hold for all appropriate functions f , we need that for almost all x
E (Wt qSt (x|Yt )) = pt (x).
We will denote the left-hand side by p̃t (x). Due to the initialisation of the CIS algorithm we have p̃0 (x) = p0 (x).
Thus we need to choose the incremental weights such that p̃t (x) satisfies the Fokker-Planck equation for the target
process
∂ p̃t (x)
= L∗ p̃t (x).
∂t
3.1.2

(7)

Obtaining the Incremental Weight

We will give an informal outline of how to derive the function ρ. Throughout this argument we will assume that
we can interchange expectation with various operators. See Fearnhead et al. (2016), and the discussion below, for
8

S_t

1.0
S_t

1.0

0.0

0.0

0.5

0.5

X_t

1.5

1.5

X_t

0

1

2

3

4

5

0

1

2

Time

3

4

5

3

4

5

1.0

1.5

●

Y_t

●
●

0.0

0.0

5

0.5

1.0

●
●

0.5

X_t

●

●

0

1

2

3

4

5

0

1

2

Time

1.4
1.3
1.2

W_t

0.8

0.9

1.0

1.1
0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

W_t

1.5

W_t

W_t

4

Y_t

1.5

X_t Simulated

3

0

1

2

3

4

5

0

1

2

Figure 1: Example of the CIS process. Left-hand column shows the original importance sampling process. Conceptually we think of a proposal process, Xt (top plot), simulated in continuous-time. However we need only simulate
this process at event times (denoted by vertical dashed lines) and at times we are interested in. The middle plot
shows these simulated values. The weight process is shown in the bottom plot, and the weight changes only at event
times. The right-hand column shows the corresponding components of our PDP process: St (top); Yt (middle) and
Wt (bottom).

9

conditions under which this is valid.
We need to choose ρ so that (7) holds. The right-hand side of (7) is just
L∗ p̃t (x) = E (L∗ Wt qSt (x|Yt )) .
The left hand side of (7) is the derivative of an expectation, and thus can be written in terms of the generator of the
PDP, A,
∂ p̃t (x)
∂E(Wt qSt (x|Yt ))
=
= E (AWt qSt (x|Yt ))) .
∂t
∂t
So both the left and right hand sides of (7) can be written as an expectation of a function of the PDP. For these
expectations to be equal it is sufficient for the two functions to be equal to each other.
Using the form of the generator for the CIS process (6) we get that for current state z = (y, w, s),


Z  
ρ(y0 , y, s)
∂wqs (x|y)
+ λ̃(s)
w 1+
q0 (x|y0 ) − qs (x|y) qs (y0 | y)dy0
Awqs (x|y) =
∂s
λ̃(s)
Z 


Z
∂qs (x|y)
ρ(y0 , y, s)
0
0
0
0
0
= w
+ wλ̃(s)
1+
q0 (x|y )qs (y | y)dy − qs (x|y) qs (y | y)dy
∂s
λ̃(s)




∂qs (x|y)
ρ(x; y, s)
= w
+ λ̃(s) wqs (x|y) 1 +
− wqs (x|y)
∂s
λ̃(s)
∂qs (x|y)
= w
+ wqs (x|y)ρ(x, y, s).
∂s
The above argument is informal, as for s = 0 (a state visited every time a jump occurs) the function qs (x|y) does not
R
exist. However it acts, informally, like a dirac delta function. So for any function h(y0 ) we have q0 (x|y0 )h(y0 )dy0 =
h(x), and this is used above. We refer the reader to Fearnhead et al. (2016) for formal justification of this step and
others.
Thus for the two functions of the PDP that we are taking expectations of to be equal we need
w

∂qs (x|y)
+ wqs (x|y)ρ(x, y, s) = L∗ wqs (x|y),
∂s

which can be re-arranged to give
ρ(x, y, s) =

(x|y)
L∗ qs (x|y) − ∂qs∂s
.
qs (x|y)

(8)

The incremental weight is then 1 + ρ(x, y, s)/λ̃(s).
The form of this incremental weight is quite intuitive. It is based on the difference between how transition densities
change under the target process and under the proposal process. The optimal proposal process would, obviously, be
the target process. For this choice ∂qs (x|y)/∂s = L∗ qs (x|y), ρ(x, y, s) = 0, and the importance sampling weights
would always be equal to 1. As expected, a proposal process that more closely mimics the target process will have
less variable incremental weights. Furthermore we see a trade-off in the choice of the event rate λ̃(s), as larger values
of this rate will lead to more events and a higher computational cost, but reduce the variance in the incremental
weight. If we were to double the event rate, we would double the expected number of events, but the variance of the
incremental weight at an event would reduce by a factor of 4. We expect that the net effect is an overall reduction
of Monte Carlo variance in the weights but an increase in computational cost.
One issue with the CIS algorithm is that, for certain combinations of target and proposal processes, it can be possible
to get negative weights (similar, for example, to the Russian roulette sampler of Lyne et al., 2015). This can occur if
it is possible for ρ(x, y, s) to be smaller than −λ̃(s). For the applications of CIS to simulating from a general diffusion
process that are described in Fearnhead et al. (2016), ρ(x, y, s) cannot be bounded below, and thus negative weights
10

are possible. Obviously the probability of getting negative weights can be controlled, with larger λ̃(s) values reducing
this probability. This in turn leads to important theoretical and practical issues. The above argument for deriving
the incremental weight required the interchange of expectation and operators. The main condition for justifying this
is that the we need the importance sampling weights to be such that E(|Wt |) is finite for any t > 0. This condition
may not hold due to the possibility of negative weights. In fact, Fearnhead et al. (2016), show that some naive
implementations of CIS will not have importance sampling weights that satisfy this condition. Furthermore, for
problems where the target process is a diffusion, they give sufficient conditions on both the proposal process and the
event rate that will ensure E(|Wt |) is finite. In general, and within the more specific context of diffusion proposals
and targets, Fearnhead et al. (2016) provide conditions under which Wt has p-th moments for p ≥ 1. The paper also
gives practically implementable and intuitive strategies for ensuring that these moments exist.

3.2

Continuous-time Sequential Monte Carlo

To date we have reviewed a continuous-time importance sampling algorithm. This is most naturally viewed as an
extension of sequential importance sampling (Liu and Chen, 1998) to continuous-time. However it is then possible
to extend this to a continuous-time version of SMC. All we need to do is to jointly simulate multiple CIS processes,
and then introduce resampling steps. The simplest implementation of this is to choose a set of resampling times, say
h, 2h, 3h, . . . , Kh, and a number of “particles”, N . The continuous-time SMC algorithm will then proceed as follows:
(1)

(N )

(0) Initiate. Simulate x0 , . . . , x0

(i)

independently from p0 (x). Set w0 = 1 for i = 1, . . . , N . Set t = 0

(1) Propagate and Reweight. For i = 1, . . . , N run the CIS process for a time interval of length h, with initial
(i)

(i)

(i)

(i)

values xt and wt . Denote the output of the CIS process at time u = t + h by x̃u and w̃u .
(1)

(N )

(2) Resample. For i = 1, . . . , N sample ki from 1, . . . , N with probabilities proportional to |w̃u |, . . . , |w̃u |. Set
(i)

(k )

xu = x̃u i and
wu(i) =

N
w̃uki 1 X (j)
|w̃u |.
|w̃uki | N j=1

(3) Iterate. Set t = u. If t < Kτ go to step (1).
The resampling step is different from standard resampling used in SMC to allow for the possibility of negative weights.
The form of the weight after resampling is defined so that resampling is unbiased, and this requires the sign of the
weight assigned to any particle value to be unchanged. It is simple to extend the above use of resampling to allow
for lower-variance resampling schemes (Kitagawa, 1996; Carpenter et al., 1999; Doucet et al., 2000) in step (2), and
to allow resampling times to depend on the SMC output, for example to be times when the effective sample size of
the weights drops below some threshold (Liu and Chen, 1995).
The beneficial effect of resampling will be less when we have negative weights. For example, it is easy to show that
E(|Wt |) is unchanged by resampling. Thus resampling cannot counteract any increase in E(|Wt |) with t, and this
increase will necessarily imply deterioration of Monte Carlo performance with t. Thus the good long-term stability
properties that SMC with resampling often has (Del Moral and Guionnet, 2001; Douc et al., 2014) will not be possible
in the presence of negative weights. This issue with negative weights is well-known within related quantum Monte
Carlo methods (Foulkes et al., 2001), where it is termed the fermion-sign problem.
Fearnhead et al. (2016) suggest alternative resampling approaches for step (2) that can reduce E(|Wt |) at resampling
times, and thus may lead to long-term stability. Alternatively, we need to choose event rates in CIS to be sufficiently
large that negative weights are rare over the time-scales that we wish to run an SMC algorithm for.
11

3.3

CIS for Big Data: the SCALE Algorithm

Recently Pollock et al. (2016a) presented SCALE, an algorithm for sampling from a posterior distribution. The
original derivation of SCALE was based on constructing a killed Brownian-motion process whose quasi-stationary
distribution is the posterior distribution. The SCALE algorithm then samples from this quasi-stationary distribution.
A key property of SCALE is that it only needs to use a small sub-sample of the data at each iteration of the algorithm,
and thus it is suitable for large data applications.
Whilst the original derivation is very different, we show here that SCALE can be viewed as a CIS, or continuous-time
SMC, algorithm. Our setting is that we wish to sample from a posterior distribution which we will assume can be
written as
π(x) ∝

n
Y

πi (x),

i=1

where, to keep notation consistent with our presentation of CIS, x is the parameter vector. Here πi (x) is the likelihood
for the ith observation multiplied by the 1/nth power of the prior. As is common to Bayesian inference, the posterior
distribution is only known up to a constant of proportionality. We wish to develop a Monte Carlo algorithm for
sampling from this posterior that has good computational properties for large n.
The idea of SCALE and its link to CIS is as follows. We will use CIS to target a stochastic process that has π(x)
as its stationary distribution. To implement CIS we only need to know the Fokker-Planck equation for this process.
If we run CIS (or in practice a continuous-time SMC version) then after a suitable burn-in period this will gives us
weighted samples from π(x). A key property of the CIS algorithm is that the incremental weights depend on the
posterior only through log π(x). This is a sum, and it is easy to use sub-sampling to unbiasedly estimate this sum.
Thus to deal with potentially large data we will implement a random weight version of CIS (Fearnhead et al., 2008),
where we use sub-sampling to estimate the incremental weights. Using unbiased random weights leads to a valid
importance sampler, but one with an increased Monte Carlo error. The next ingredient to the SCALE algorithm
is to use control variates to reduce the variance of our sub-sampled estimates of the incremental weights, which in
turn helps to control the overall Monte Carlo error of the algorithm. Finally Pollock et al. (2016a) use ideas from
the Exact Algorithm to avoid the possibility of negative weights. We now detail each of these steps.
The first step to SCALE is the choice of a stochastic process that has π(x) as its asymptotical distribution. We need
specify this process through its Fokker-Planck equation. SCALE uses the stochastic process for which
!
d
d
X
1 X ∂ 2 pt (x)
1
∂ 2 π(x)
∂pt (x)
=
−
pt (x).
∂t
2 i=1 ∂x2i
2π(x) i=1 ∂x2i
It is trivial to see that π(x) is an invariant distribution for this stochastic process, as on substituting pt (x) = π(x) the
two terms on the right-hand side cancel. The actual underlying stochastic process can be interpreted as Brownian
motion with killing, conditioned on survival: see Pollock et al. (2016a) for more details.
If we implement CIS for this target process, and use Brownian motion as the proposal distribution, we have
( d
)

d
X (xi − yi )2
1
qs (x|y) = √
exp −
,
2s
2πs
i=1
and
d

∂qs (x|y)
1 X ∂ 2 qs (x|y)
=
.
∂s
2 i=1
∂x2i
Thus from (8), the function that determines the incremental weights is
"

2 #
d
d
1 X ∂ 2 π(x)
1 X ∂ 2 log π(x)
∂ log π(x)
ρ(x, y, s) = −
=−
+
.
2π(x) i=1 ∂x2i
2 i=1
∂x2i
∂xi
12

The right-hand side depends on π(x) only through derivatives of
log π(x) = constant +

n
X

log πi (x).

i=1

Importantly these derivatives do not depend on the unknown normalising constant of π(x). Furthermore as they are
sums, it is simple to unbiasedly estimate the derivatives. For example given j and k, two independent draws from
an uniform distribution on 1, . . . , n, we can estimate ρ(x, y, s) by



d 
∂ log πj (x)
1 X ∂ 2 log πj (x)
∂ log πk (x)
2
+n
−
n
.
2 i=1
∂x2i
∂xi
∂xi

(9)

We can reduce the variance of our estimate of ρ, and hence of the incremental weights, using control variates. Pollock
et al. (2016a) suggest using a pre-processing step that finds a values x̂ close to the posterior mode (for similar ideas,
see Bardenet et al., 2015). We precalculate the first and second derivates of log π(x) at x̂, and calculate the value,
ρ̂, of ρ at x̂. We then estimate of ρ at x as
d

nX
−
2 i=1







∂ log πj (x) ∂ log πj (x̂) ∂ log πj (x) ∂ log πj (x̂)
∂ 2 log πj (x) ∂ 2 log πj (x̂)
1 ∂ log π(x̂)
−
+n
−
−
+2
+ρ̂
∂x2i
∂x2i
∂xi
∂xi
∂xi
∂xi
n
∂xi
(10)

where, as before, j and k are independent draws from an uniform distribution on 1, . . . , n. The idea behind this
√
control-variate approach is that if x̂ is within 1/ n of the posterior mode, then with high-probability, at stationarity
√
the CIS process will be at an x value that is within 1/ n of x̂. If the first and second derivatives of the log πj s are
√
well-behaved this means that the terms in the square brackets will be Op (1/ n), and thus ρ will be O(n) with high
probability. Thus compares well with the naive subsampling estimator of ρ (9), which will be Op (n2 ).
Finally, to avoid negative weights, the SCALE algorithm then uses ideas from Beskos et al. (2008) and Burq and
Jones (2008), to simulate the proposal process in such a way that we know an upper and lower bound the process
takes within a given time-interval. With such bounds we can then choose our event-rate λ sufficiently high that
negative weights do not occur. This approach does come with a computational cost, as simulating Brownian motion
together with such a bound can be an order of magnitude, or more, slower than just simulating Brownian motion.
Below we investigate the feasibility of implementing a version of SCALE that allows for negative weights, but that
chooses event-rates to be sufficiently large that they are rare.

3.4

Extensions

Here we give two examples of how we can use the theory for PDPs to easily obtain generalisations of the basic CIS
and SCALE algorithms.

3.4.1

Alternative Proposal Distributions

With any importance sampling approach, the choice of proposal distribution can have a substantial impact on the
resulting Monte Carlo properties. The original derivation of CIS used the idea of a continuous-time stochastic process
that was being used as a proposal process. However our derivation of the CIS algorithm has not actually required
us to specify a proposal process, just a suitable family of transition densities, qs (x|y). This family needs to have
certain properties, such as being differentiable with respect to s, so that the incremental weight (8) can be calculated.
This appears to make it easy to consider alternative proposals, as we only need to specify an appropriate family of

13

densities. For example, in standard importance sampling applications it is often recommended to use heavy-tailed
proposals, so for CIS it is natural to consider transition densities that are t-distributed, such as
qs (x|y) ∝ s−d/2

1
(1 +

1
ν

Pd

i=1 (xi

− yi )2 /s)(ν+d)/2

,

(11)

for some appropriately chosen degrees of freedom ν > 0.

3.4.2

Alternative SCALE Processes

We can also develop alternatives to the SCALE algorithm that differ in terms of the underlying stochastic process
that they target. For example we can target the Langevin diffusion, for which the Fokker-Planck equation is


d
d 
∂pt (x)
1 X ∂ 2 pt (x) 1 X ∂
∂ log π(x)
−
=
p
(x)
.
t
∂t
2 i=1 ∂x2i
2 i=1 ∂xi
∂xi
Again, it is simple to see, by substituting pt (x) = π(x) on the right-hand side, that π(x) is the invariant distribution
for this process. If we implement CIS to sample from this Langevin diffusion, and use Brownian Motion as the
proposal distribution, we get that the incremental weights are 1 + ρ(x, y, s)/λ̃(s) with
ρ(x, y, s) = −


d 
1 X (yi − xi ) ∂ log π(x) ∂ 2 log π(x)
+
.
2 i=1
s
∂xi
∂x2i

We can unbiasedly estimate this using a sub-sample of size 1. For example, if j is drawn uniformly from {1, . . . , n},
one unbiased estimator is

d 
n X (yi − xi ) ∂ log πj (x) ∂ 2 log πj (x)
.
ρ̂(x, y, s) = −
+
2 i=1
s
∂xi
∂x2i
Alternatively we can develop lower-variance esimators using control variates. The incremental weight is easier to
estimate than the incremental weight in SCALE as it does not involve a square of the gradient of the log-posterior.
However, the (yi − xi )/s term has a variance, under the Brownian motion proposal, that is 1/s, and care is needed
to control this term for small inter-event times, s.

3.5

Example: Inference for Mixture Models

To give an indication of how the SCALE algorithm of Section 3.3 works, some of the issues with its implementation,
and the importance of using control-variates with the sub-sampling estimator or ρ, we will consider a simple example.
Assume we have IID data from a mixture distribution, where for j = 1 . . . , n

 N(0, 102 ) with probability p,
Yj ∼
 N(x, 12 ) otherwise.
Our interest is inference for x, and we assume a Gaussian prior with mean 0 and variance 4. In the following we
simulate data with x = 4.
Initially we implemented the SCALE algorithm without any subsampling on a data-set of size n = 150 and p = 0.95.
Figure 2 show the posterior distribution for x for this data set, the log-posterior and how ρ varies as a function of
x. Note that the average value of ρ at stationarity is 0, so values where ρ is greater than or less than 0 show regions
where weights of particles will increase or decrease relative to the average weight. The key point in these figures is
that while both the posterior and log-posterior are relatively well-behaved, ρ has multiple pronounced modes.
14

(b)

(c)

rho

−6

0.0

−10

−6

0.1

−8

−4

0.2

density

log−density

0.3

−2

−4

0.4

0

−2

0.5

(a)

−5

0

5

−5

0

x

5

x

−5

0

5

x

Figure 2: The posterior distribution (a); the log-posterior (b); and ρ (c) for a data set of size 150 from the mixture
model.

This multi-modality can cause issues with mixing. For example a particle that is currently at a value of −5 will have
to traverse a prolonged region where ρ is negative to reach the main region of posterior mass at values of x close to
4. As it traverses this region its weight will reduce substantially, and it is likely to be lost due to resampling. Thus
we could have relatively rare movement of particles from one mode of ρ to another.
In practice this means that initialisation of the SCALE algorithm can be crucial. We implemented SCALE with 200
particles over a time interval of length 100. We considered resampling at every integer time-point, and resampled
if the effective sample size of the weights was less than 100. We ran SCALE with a constant event rate of 12, and
observed no negative weights. The performance of SCALE appeared quite sensitive to the event rate, with rates of
10 and less giving noticeably worse performance. Figure 3 shows output from two runs of SCALE, one initialised
with particles drawn from the prior, and the other initialised with particles drawn uniformly on [−10, −5]. Figure
3 (a) and (c) show the evolution of the particles over the first 50 time-steps. They evolve according to Brownian
motion, but accrue or lose weight depending whether ρ is positive or negative at the particle values at event times.
Particles with low weights tend to be lost at resampling times. For the case where we initialise particles from the
prior, we see particles are quickly lost in regions aroung x = 0 where ρ is negative. More slowly they are also lost
from regions around x = −5 which is a local-maxima in ρ, as these particles, on average, attain a smaller weight
than those particles close to the posterior mode. By about time 15 the SCALE algorithm appears to have converged,
and Figure 3 (b) shows that it gives a good approximation to the true posterior distribution. By comparison, when
SCALE is initialised with no particles close to the posterior mode, the particles appear to get stuck close to local
mode of ρ at x = −5.
We now turn our attention to the possibly benefits of using subsampling to estimate ρ. Firstly to show the difference
between subsampling with and without control variates, Figure 4(a) shows the variance of our the sub-sampling

15

(b)

0.3
0.2
0.0

0.1

−5

x

0

density

0.4

5

0.5

(a)

10

20

30

40

50

−4

−2

0

2

Time

x

(c)

(d)

4

6

0.2
0.0

−5

0.1

x

0

density

0.3

5

0

0

10

20

30

40

50

−10

Time

−8

−6

−4

−2

0

2

x

Figure 3: Evolution of particles up to time 50 (left-hand column) and estimate of posterior (right-hand column)
from SCALE algorithm. Top row is for particles initiated from prior and bottom row for particles initiated uniformly
on [−10, −5]. Estimates of posterior are show as histograms based on the weighted particles from time 25 to 100,
and the true posterior is overlain.

16

Variance

0.04

0.06

10000
1000

Data Points

100

250
200
150

0.00

1

0

50

10

0.02

100

Standard Deviation

(c)

0.08

(b)

300

(a)

2

3

4

5

6

x

2e+03

1e+04

5e+04
n

5e+05

2e+03

1e+04

5e+04

5e+05

n

Figure 4: Variance of estimate of ρ using subsampling with (black full-line) and without (red dashed line) control
variates (plot a); the vertical line shows the value of x̂. Computational cost of estimating Wh , measured in terms of
number of data point accessed, as a function of n (plot b), and variance of Wh as a function of n (plot c). For these
latter two plots, results are shown for no subsampling (red dashed line) and subsampling with x̂ as the posterior
mode (black full-line) the mode plus posterior standard deviation (green dotted line) and mode plus three times the
standard deviation (blue dot-dash line). Results calculated from 2, 000 estimates of Wh for each method and each
value of n.

estimate of ρ for our data set of size 150. The control variate estimator (10) was implemented with x̂ set to be the
posterior mode. We see a substantially lower variance for x close to x̂ when we use control variates. Though as x
moves away from x̂ the variance increases, and eventually is worse than not using control-variates. The key to why
control variates works for large data is that if x̂ is close to the posterior mode then x will be very close to x̂ with
high-probabililty. In such cases the proportionate reduction in variance will be O(n) – and thus the gains of using
control variates will increase for larger data sets.
To gain insight into the benefit of using sub-sampling we looked at the variance of the weight at time h, Wh , of
a particle sampled from the true posterior, for different values of n both with and without subsampling. For a
fixed variance of Wh we wanted to see how many data-points need to be processed by each algorithm. So without
sub-sampling, each event requires access to all n data points, whereas with sub-sampling each event requires us to
access only 2 data points. We found that sub-sampling without control variates performed substantially worse than
the two alternatives, with increasingly poor performance as n increases. Thus we focus on comparing no subsampling
and subsampling with control variates.
√
As we increase n the posterior standard deviation will decrease at a rate 1/ n, and thus we choose h = 1/n so that
√
the distance moved by the particle will also be of order 1/ n. When using no-subsampling we chose the event rate
to be n/2. For sub-sampling with control variates we set x̂ to be the posterior mode and let the event rate depend
on the value of the particle at the most recent event, x0 say, and be of the form 2n + 4n2 (x0 − x̂)2 . These choices
gave variances of Wh that were similar for both cases and also for different n. Figure 4(b) shows how the number of
data-points accessed varied with n for the two methods and Figure 4(c) shows the estimates of the variance of Wh .
As n increases we have that the computational cost of using sub-sampling, as measured by the number of data-points
accessed, remains constant. By comparison, without sub-sampling, to maintain a fixed variance for Wh we need to
increase the computational cost linearly.
Finally we looked at how the choice of x̂ affected the performance of subsampling with control variates. We repeated

17

the above study but setting x̂ to be either one posterior standard-deviation or three posterior standard-deviations
from the mode. These correspond to values in the body of the posterior and in the tail of the posterior respectively.
From 4 (b) and (c) we see that as x̂ moves away from the posterior mode, the performance of the sub-sampler
decreases both in terms of computational cost and variance. However, for a fixed variance of Wh both methods have
a computational cost that is constant with n, as opposed to the linearly increasing cost we have when sub-sampling
is not used.

4

Continuous-Time MCMC

We now consider continuous-time versions of MCMC. These algorithms involve simulating a PDP process which
has a given target distribution, π(x), as its stationary distribution. Such algorithms were originally of interest as
they are non-reversible processes. As mentioned in the introduction, there is substantial evidence that non-reversible
MCMC algorithms are more efficient than standard, reversible MCMC. Intuitively this is because non-reversible
MCMC suppresses the random-walk behavour of reversible MCMC and thus can more rapidly explore the statespace. Furthermore it has been shown that continuous-time MCMC is suitable for using sub-sampling ideas, similar
to those in Section 3.3. Thus these methods are also promising for big-data applications of MCMC.

4.1

The Continuous-time limit of MCMC

To help build intuition for continuous-time MCMC, and to see how it links to discrete-time MCMC algorithms,
we will first derive a continuous-time algorithm as a limiting form for a simple non-reversible discrete-time MCMC
algorithm (Gustafson, 1998; Diaconis et al., 2000). This MCMC algorithm will target a joint distribution of (x, v),
where v can be viewed as a velocity. For our specific algorithm we will consider only velocities of a fixed, say unit,
speed, and hence v could equally be defined as a direction. Our MCMC will target a distribution π(x)pu (v) where
pu (v) will be the uniform distribution over all velocities with unit speed.
The MCMC algorithm will have two types of move. The first involves two deterministic proposals
(1a) Propose a move from (x, v) to (x + hv, −v). Accept this with the standard Metropolis-Hastings accept probability, which simplifies to


π(x + hv)
min 1,
π(x)
(1b) Move from (x0 , v0 ) to (x0 , −v0 ).
Both the moves in (1a) and (1b) are reversible, and can be shown to satisfy detailed balance. To make step (1a)
reversible we have to propose a move which flips the velocity, and hence in (1b) we flip the velocity back again. So
the net affect of applying both (1a) and (1b) is that the velocity is unchanged if we accept the proposed move in step
(1a) but flips if we reject the move. This is a standard approach in Hamiltonian Monte Carlo (Neal et al., 2011).
In fact this algorithm can be viewed as a type of Hamiltonian Monte Carlo move, but based on the dynamics of an
approximate potential for x which is uniform (and hence the velocity is not changed other than by the flip).
Whilst this moves keeps π(x)pu (v) invariant, it leads to a reducible Markov chain if the dimension of x is greater
than 1, as it only proposes moves along the direction given by v. Thus we need a second type of move to produce an
irreducible MCMC algorithm with the required asymptotic distribution. The second move we use is an update of v,
from some transition kernel that has pu (v) as its stationary distribution. We will imagine applying N transitions of
type 1 between each of these updates just of v.
18

Under this framework we can then consider letting h → 0 while keeping s = hN a constant. We will scale time so
that the ith MCMC transition will occur at time ih, and define (xt , vt ) to be the value of the state after the ith
MCMC transition for ih ≤ t < (i + 1)h.
Now for each move in step (1a) the rejection probability for small h is
max {0, 1 − exp[log π(x + hv) − log π(x)]} = max {0, 1 − exp[v · ∇ log π(x)h + o(h)]} = max {0, −v · ∇ log π(x)h}+o(h),
assuming that, for example, π(x) is twice differentiable.
Thus in our limit as h → 0, rejections in step (1a) will occur as events in a Poisson process of rate λ(xt , vt ) =
max{0, −vt · ∇π(xt )}. The dynamics between these events will be deterministic, with vt being constant and xt
changing as in a constant velocity model with velocity vt . At each event the velocity will just flip. Note that while
the process is moving to areas of higher probability density, as defined by π(x), the rate of the Poisson process will
be 0. Thus events will only occur if the process is moving to areas of lower probability mass.
This limiting process is just a PDP. The dynamics of the PDP are a constant velocity model between events, with
the velocity changing at event times. It is natural to consider a more general class of PDP processes, and see what
flexibility there is in choosing the distribution of the event times, and the distribution of the change of velocity at
events, so that we still have a process whose marginal stationary distribution for Xt is π(x).
Denote the state of our PDP by Zt = (Xt , Vt ). Our class of PDPs will have the following dynamics:
(i) The deterministic dynamics. For i = 1, . . . , d
(i)

(i)

dxt
dv
(i)
= vt , and t = 0.
dt
dt
The solution of these dynamics is given by (xt+s , vt+s ) = (xt + svt , vt ) for any s > 0.
(ii) The event rate. Events will occur at a rate, λ(zt ).
(iii) The transition distribution at events. At an event at time τ , xτ = xτ − and vτ is drawn from some
transition density q(·|xτ − , vτ − ).
The event rate and the transition density then need to be chosen so that π(x) is the marginal stationary distribution
of the PDP.

4.2

The Stationary Distribution of the PDP

A necessary condition for π(x) to be the marginal stationary distribution of our PDP is that it is the marginal of an
invariant distribution for the PDP. We will use the adjoint of the generator of our PDP to derive a condition on both
the event rate and the transition distribution at events for the PDP to have π(x) as the marginal of an invariant
distribution.
As above, let z = (x, v). Denote the invariant distribution of our PDP by p(z). We can factorise this as the product of
the marginal stationary distribution for x times the conditional for v given x, and we wish to have p(z) = π(x)p(v|x).
If A∗ is the adjoint of the generator of our PDP, as p(z) in an invariant distibution we have A∗ p(z) = 0. This gives
Z
− π(x)p(v|x)[v · ∇x log π(x) + v · ∇x log p(v|x) + λ(z)] + λ(x, v0 )q(v|x, v0 )π(x)p(v0 |x)dv0 = 0.
(12)
In the above ∇x denotes the vector of first partial derivative with respect to the components of x.
19

To date, all continuous-time MCMC algorithms have been designed so that under the invariant distribution v is
independent of x, and thus all components of ∇x log π(v|x) will be 0. If we wish to design such a process we need
to choose λ(x, v) and q(v0 |x, v) such that (by rearranging equation 12),
Z
pv (v)λ(x, v) − λ(x, v0 )q(v|x, v0 )pv (v0 )dv0 = −pv (v)v · ∇x log π(x).

(13)

for some distribution pv (v) for the velocity. The left-hand side is measuring the net probability flow out of states
with velocity v, this must offset the change in probability mass for V caused by the deterministic dynamics, which
is the term on the right-hand side.
Note that if we integrate (13) with respect to v, the left-hand side is 0. So we get E(V) · ∇ log π(x) = 0, where the
expectation is with respect to the invariant distribution for the velocity. As this will need to hold for all x, we can
see that the invariant distribution for all components of the velocity must have zero mean.
The actual processes we describe in the next section all allow velocities within some symmetrical set, and are
designed so that pv (v) is uniform on this set. They ensure (13) holds through deterministic dynamics at events.
They introduce a “flip” operator, Fx say, that can depend on x and which satistifies Fx (Fx (v)) = v. They then only
allow transitions between pairs of velocities, v and v0 that satisfy v0 = Fx (v) and, by definition of Fx , v = Fx (v0 ).
Under this constraint on the transitions at events we get a simple set of equations that we need the event rates to
satisfy. For any v, and with v0 = Fx (v), it is straightforward to show that (13) requires
λ(x, v) − λ(x, v0 ) = −v · ∇x log π(x).

(14)

for all x. Note that as this equation must also hold for v0 , we immediately see that only flip operators for which
Fx (v)∇x log π(x) = −v∇x log π(x) are allowable. The rates only depend on the target through the term ∇x log π(x),
which means that π(x) is only needed to be known up to proportionality. Also note that (14) does not uniquely
define the rates. If we have a set of rates, λ(x, v) that satisfy (14), then λ(x, v) + γ(x, v) will also satisfy (14) for
any positive function γ(x, v) for which γ(x, v) = γ(x, Fx (v)).
A natural choice of rates which satisfy (14) are those which are smallest. This will give λ(x, v) = max{0, −v ·
∇x log π(x)}. We will call these the canonical rates. Theoretical justification for the canonical rates when d = 1 is
given in Bierkens and Duncan (2016), who show that the asymptotic variance of Monte Carlo estimators is minimised
when using these rates.

4.2.1

Different continuous-time MCMC algorithms

We now describe a number of choices of flip operator, and the corresponding PDPs. We start with the limiting
process we derived in Section 4.1, then describe two continuous MCMC processes that have been recently proposed.
These three choices all lead to identical PDPs for a one-dimensional target, but differ in terms of how they extend
to higher dimensions. Each assume the target is defined on an unbounded domain; for extensions of these methods
to bounded domains see Bierkens et al. (2016). We will finish with some discussion of alternative schemes that are
possible.
Pure Reflection and Refresh
The continuous-time limit we derived in Section 4.1 corresponds to Fx (v) = −v, with the canonical rates. Such a
process we call a pure reflection process. For a multi-dimensional target distribution, this process would be reducible,
as it can only explore positions x that lie on a straight-line defined by the initial velocity. As such this is an example
where π(x) would be a marginal invariant distribution but not the marginal stationary distribution. To overcome
20

this we would need an additional move, which refreshes v. Such a refresh move would need to have pu (v) as its
stationary distribution. The times of refreshing could be either deterministic or random.
Bouncy Particle Sampler
The Bouncy Particle Sampler of Bouchard-Côté et al. (2015), based on an algorithm of Peters and de With (2012),
is an adaption of the pure reflection process, which minimises the change in velocity at each event. It does this by
defining Fx (v) to be
Fx (v) = v − 2

v · ∇ log π(x)
∇ log π(x).
∇ log π(x) · ∇ log π(x)

(15)

This flips the component of v that is in the direction of ∇ log π(x) but leaves the components of v that are orthogonal
to ∇ log π(x) unchanged. They again use the canonical rates. As with the pure reflection process this means that
events only occur if the PDP is moving to areas of lower probability mass according to π(x).
The original sampler of Peters and de With (2012) just implement that above sampler. However Bouchard-Côté
et al. (2015) shows that, for some targets, such a sampler can be reducible. The Bouncy Particle Sampler introduces
a refresh step, which occurs at events of an independent Poisson process of constant rate. Furthermore they show
that for any non-zero rate of this refresh process, the resulting sampling will have π(x) as its asymptotic distribution.
Zig-Zag Sampler
The Zig-Zag sampler (Bierkens and Roberts, 2015; Bierkens et al., 2016) considers a discrete set of velocities. If x
Pd
is d-dimensional, then v = i=1 θi ei , where each θi ∈ {−1, 1} and e1 , . . . , ed are a set of orthongal basis vectors for
Rd . The invariant distribution for v is defined as the uniform distribution over this set of 2d possible values.
The Zig-Zag sampler can be viewed as having d-distinct event types, each with its own rate, and each with its own
deterministic change to the velocity. The ith event will have a flip, F (i) , that switches θi to −θi , but keeps the velocity
in the other d − 1 directions unchanged. If we denote λi (x, v) to be the rate of events of type i, then this corresponds
Pd
to our general formulation of a PDP but with λ(x, v) = i=1 λi (x, v), and with the transition distribution at an
event being a discrete distribution over the d transitions that correspond to the d different flips. Flip i occurs with
probability λi (x, v)/λ(x, v). Subsituting this into (13) shows that we need to choose λi (x, v) so that
d n
d
o
X
X
∂ log π(x)
(i)
.
λi (x, v) − λi (x, F (v)) = −
θi
∂x(i)
i=1
i=1

Here we have assumed that x(i) is the component of x in direction ei . This can be achieved if we choose the rates
such that
∂ log π(x)
.
∂x(i)
As above, this does not uniquely define the rates, only the difference between rates for velocities that differ in terms
λi (x, v) − λi (x, F (i) (v)) = −θi

of their component in the ei direction.
It is a challenging goal to show that the Zig-Zag process is ergodic in full generality. So far it is established rigorously
in Bierkens et al. (2016) that the Zig-Zag process is ergodic in any of the following cases: (i) one-dimensional target
distributions, (ii) factorized target distributions, and (iii) switching rates that are positive everywhere (which can be
obtained by adding a constant ε > 0 to the canonical switching rates). Experiments suggest that ergodicity holds in
much more generality.
Note that the above argument easily generalises to allowing velocities of the form v =

Pm

i=1 θi ei ,

where the ei are

not constrained to be orthogonal, and we can even allow m > d directions. Whether there are advantages in using
such a set of possible velocities is not clear.
Alternatives
21

There is substantial extra flexibility in choosing the event rates and the type of transition at events beyond the
three examples we have detailed. For example we could consider transitions at an event that does not depend
on the current velocity. If we allow v to be any unit vector, then it is straightforward to show that choosing
λ(x, v) = max{0, −v · ∇x log π(x)}, and, at each event, sampling a new velocity from the distribution
q(v0 |x, v) ∝ max{0, v0 · ∇x log π(x)},
will lead to a PDP with invariant distribution that has π(x) as its marginal.
More substantial alternatives are also possible. For example, we could consider processes which allow the invariant
distribution of V to depend on x – something that Girolami and Calderhead (2011) has shown to be beneficial for
Hamiltonian Monte Carlo methods. For a proposed distribution π(v|x) we would then need to find a set of event
rates and transitions that satisfy (12).

4.3

Simulation and Use of Skeletons for Continuous MCMC

So far we have described a number of different PDPs that will have π(x) as their marginal invariant distribution.
For these to be useful in practice, we need to be able to simulate them efficiently. How to do this in practice will
depend on the form of π(x), but is likely to use the ideas briefly described at the end of Section 2.1. For further
detail see the discussion of this, and suggestions, in Bouchard-Côté et al. (2015) and Bierkens et al. (2016).
The output of simulating a PDP will be a set of event times and the values of the state at those event times. We
wish to use this output to obtain Monte Carlo estimates of expectations of functions of X, where X is distributed
according to π(x). Assume we have simulated the PDP for some time-interval T . We will discard the value of the
process in some burn-in period of length tb . Assume there were N events in the time-interval [tb , T ]. Denote these
as τi for i = 1, . . . , N , and let τ0 = tb and τN +1 = T .
There are two approaches to obtain a Monte Carlo estimate of

R

π(x)g(x)dx for some function g(x) of interest. The

first is to calculate the average of this function along the path of the PDP:
1
τN +1 − τ0

N Z
X
i=0

τi+1 −τi

g(xτi + svτi )ds.

0

Here each integral corresponds to the integrals of g(xt ) for t in [τi , τi+1 ], and uses the fact that for such a t,
xt = xτi + (t − τi )vτi .
The above approach is difficult if the integrals are not easy to evaluate. In this case we can resort to a standard
Monte Carlo approximation. Choose an interger M > 0, define h = (τN +1 − τ0 )/M and then use the Monte Carlo
estimator

M
1 X
g(xτ0 +jh ),
M j=1

where, as above, we can trivially calculate xτ0 +jh using the set of event times and the values of the PDP at those
event times.

4.4

Exact Approximation versions and Subsampling

Exact approximate algorithms (Andrieu and Roberts, 2009) are MCMC algorithms that use estimators of the target
distribution within the accept-reject step. If implemented correctly, and if these estimators are both positive and
unbiased, then it can be shown that the resulting MCMC algorithms are exact: in the sense they still have the target
22

distribution as their stationary distribution. It turns out that exact approximate versions of the continuous-time
MCMC algorithms detailed in the previous section are also possible.

4.4.1

Exact Approximation for Pure Reflection and Zig Zag

For concreteness and ease of presentation we will consider an exact approximate version of the Pure Reflection
process. Though the ideas we detail extend trivially to the Zig-Zag sampler (and see Bierkens et al., 2016, for more
details of an exact approximate version of Zig-Zag).
For the Pure Reflection process the requirement on the rates of events is that for any velocity v
λ(x, v) − λ(x, −v) = −v · ∇ log π(x).
For a given choice of rates, such as the canonical rates λ(x, v) = max{0, −v ·∇ log π(x)}, we would often use thinning
to simulate the event times (see Section 2.1). Thus if our current state is (xt , vt ) we would introduce a bound on
the event rate for s > 0
λ̃+ (s) ≥ λ(xt + svt , vt ),
simulate potential events at rate λ̃+ (s) and accept them with probability λ(x + sv, v)/λ̃+ (s). The time until the
next event is just the time until the first accepted event.
Now assume we have a estimator of −∇ log π(x), which we will denote U(x). This estimator is a random variable,
and examples of how it could be constructed are given below. We further introduce a random rate function
λ̂(x, v) = max{0, v · U(x)}.
This is just the canonical event rate, but replacing −∇ log π(x) with its unbiased estimator. The idea of an exactapproximate version is to simulate events using thinning, with a bound on the event rate that satisfies
λ̃+ (s) ≥ λ̂(xt + svt , vt )
almost surely, and where we accept points using the random acceptance probability λ̂(xt + svt , vt )/λ̃+ (s).
As the overall acceptance probability will be the expectation of the random acceptance probability, it is straightforward to show that simulating events in this way is equivalent to simulating events at a rate
λ(x, v) = E (max {0, v · U(x)}) ,

(16)

where expectation is with respect to the random variable U(x). Furthermore if we now calculate the difference in
rates, λ(x, v) − λ(x, −v), we have
λ(x, v) − λ(x, −v)

=

E (max {0, v · U(x)}) − E (max {0, −v · U(x)})

=

E (max {0, v · U(x)} − max {0, −v · U(x)})

=

E (max {0, v · U(x)} + min {0, v · U(x)}) = E (v · U(x)) .

Thus provided U(x) is an unbiased estimator of −∇ log π(x), the resulting process will have π(x) as its marginal
invariant distribution.
Note that using an unbiased estimator of −∇ log π(x) does not come without cost, as the resulting process will, in
general, be less efficient. This loss of efficiency comes first as the rates that are used for events, E(max{0, v · U(x)}),
will, in general, be larger than the canonical rates. The only exception being if U(x) has the same sign as −∇ log π(x)
23

with probability one. Using larger rates appears to reduce the rate of mixing of the process (see Bierkens and Duncan,
2016; Bierkens et al., 2016). A related issue is that the bound on the rates, used when implementing thinning, will
also tend to be larger. This will increase the cost of simulating the process. Intuitively we would expect both these
losses of efficiency to increase as the variability of our estimator increases.

4.4.2

Exact Approximation for the Bouncy Particle Sampler

The idea of an exact approximation version of the Bouncy Particle Sampler is slightly more complicated, due to the
fact that the flip operator used also depends on ∇ log π(x). To implement such an exact approximation version we
need to use the same estimate of ∇ log π(x) in the flip operator as was used in deciding to accept the event. So,
using the notation of the previous section, at a potential event at time s, simulated from a Poisson process with rate
λ̃+ (s) we will now:
(1) Simulate u, a realisation of U(xt + svt ).
(2) Accept the event with probability
λ̂(xt + svt , vt )
max{0, u · vt }
=
+
λ̃ (s)
λ̃+ (s)
(3) If we accept the potential event, this corresponds to an event at time t+s, with new state being xt+s = xt +svt
and
vt+s = vt − 2

vt · u
u.
u·u

Note that the same realisation, u, is used in both steps (2) and (3).
This leads to an algorithm where the transition at an event is random. We will assume that U is a discrete random
variable, which is consistent with the use of sub-sampling discussed in Section 4.4.3. It is straightforward to show
that the resulting PDP has events occurring at rate (16) as before, but with a transition probability mass function
at an event that is
q(v0 |x, v) = p(u|x)

max{0, u · v}
E (max {0, v · U(x)})

where p(u|x) is the probability mass of simulating u, and v0 = v − 2u(v · u)/(u · u). Now substituting these values
into the equation for the stationary distribution of the PDP (12), it is simple to show that the resulting process will
have an invariant distribution π(x)pv (v), where pv (v) is uniform on the set of velocities of fixed speed.

4.4.3

Use of Subsampling

An example of an exact-approximate version of these continuous-time MCMC algorithms arises if we use sub-sampling
of data points at each iteration when performing Bayesian inference in a big-data setting. As in Section 3.3, we will
consider a target density of the form
π(x) ∝

n
Y

πi (x).

i=1

In this case, the simplest unbiased estimator of −∇ log π(x) is just
− n∇ log πI (x),

(17)

where I is drawn uniformly from 1, . . . , n. However, to increase the efficiency of sub-sampling methods we would
want to try and minimise the variance of our estimator, for example by using control variates. As in Section 3.3,

24

one approach is to have a pre-processing step that finds x̂, a value of x that is near the posterior mode. We then
calculate and store ∇ log π(x̂), and use the following estimator
− ∇ log π(x̂) + n (∇ log πI (x̂) − ∇ log πI (x)) .

(18)

√
By a similar argument to that in Section 3.3, if x is within a distance of O( n) of x̂ and if π(x) is sufficiently smooth
we would expect the variance of this estimator to only increase linearly with n. By comparison the variance of the
simple estimator (17) will increase like O(n2 ).

4.5

Example: Mixture Model

To demonstrate some of the properties of these continuous-time MCMC algorithms we apply them to the simple
mixture model of Section 3.5. This involves inference for a univariate parameter, and for this case each of the three
versions of continuous-time MCMC introduced earlier are equivalent. Furthermore we do not need to introduce any
refreshing of the velocity to ensure ergodicity (Bierkens and Duncan, 2016). Our aim is to show how the continuoustime MCMC algorithms can be implemented, how the choice of the bounding process that simulates potential event
times can effect efficiency, and give some insight into how and when the subsampling ideas can lead to gains in
efficiency.
We will first look at implementing continuous-time MCMC using a global bound on the event rates. For this model,
Qn
if we write π(x) ∝ i=1 πi (x), where πi is the likelihood of the ith observation times the 1/nth root of the prior,
then

log πi (x) = log





p
1 2
1
1
exp −
yi + (1 − p) exp − (x − yi )2
− x2 .
10
200
2
8

As in Section 3.5 we will fix p = 0.95 in the simulations we present. We can bound |∇ log πi (x)| for each i, and this
bound increases with |yi |. If we let j be the observation with largest absoluate value, then the simplest global bound
on the event rates will be
λ+ = n max |∇ log πj (x)|.
x

(19)

We can then simulate the path of our continous-time MCMC algorithm by iterating the following steps. Assuming
we are currently at time t with state (xt , vt ):
(1) Simulate the time until the next putative event, s, a realisation of an exponential distribution with rate λ+ .
(2) Calculate xt+s = xt + svt , and the actual rate of an event at position xt+s :
λ(xt+s , vt ) = max{0, −vt ∇ log π(xt+s )}.
(3) With probability λ(xt+s , vt )/λ+ switch the sign of the velocity, vt+s = −vt and store the value (xt+s , vt+s ).
Otherise vt+s = vt .
This simulates using the canonical rates. To use the simplest version of sub-sampling we just replace the calculation
of the actual rate in step (2) by
λ(xs+t , vt ) = n max{0, −vt ∇ log πI (xs+t )},
for I sampled uniformly from {1, . . . , n}. Note that our choice of λ+ can still be used with sub-sampling, as it bounds
the above rate for all I and xs+t .
The above algorithm has similarities to one iteration of a standard MCMC algorithm. Steps (1) and (2) can
be viewed in terms of simulating a new state, and step (3) is a form of accept-reject step. However there are
25

fundamental differences. Firstly, the probability in step (3) depends on the target through ∇ log π(x), as compared
to the acceptance probability in MCMC which depend on π(x). Secondly, the algorithm moves from xs to xs+t
regardless of the outcome in step (3). Step (3) only affects the velocity component. Finally, as mentioned in Section
4.3 the use of the output is different. For continuous-time MCMC we have to take averages with respect to the
continuous-time path, or with respect the value of the process at equally-spaced time-points. By comparison MCMC
would average with respect to the value of the chain at the end of each iteration.
We now turn to how the use of sub-sampling impacts on the efficiency of continuous-time MCMC. For our above
implementation the average number of iterations needed to simulate a path over a time-interval of length T will be
T /λ+ for both the canonical and sub-sampling versions. Sub-sampling will involve a smaller cost in step (2) as the
rate depends on just a single data point rather than all n data point. However this computational saving comes
at the cost of an overall increased rate of switching velocity. This is shown in Figure 5, where we give examples of
∇ log π(x) for two simulated data sets, of size n = 150 and n = 1, 500 respectively, and each simulated with the true
value of x = 4. We also show the canonical rate of switching from a negative to a positive velocity, and the expected
rate of switching when we use sub-sampling. The canonical implementation has uniformly lower rates.
The impact of these different rates can be seen in Figure 6, where we show trace autocorrelation plots for analysing
the data set with n = 150. Using subsampling leads to paths of the sampler that switch velocity substantially
more frequently. As a result, the canonical implementation is more efficient in terms of suppressing random walk
behaviour, and this is seen in terms of better mixing and lower autocorrelation. The autocorrelation plots suggest we
need to run the sub-sampling version for roughly 5 times as long to obtain the same accuracy as using the canonical
implementation.
We can improve the computational efficiency of both these implementations of continuous-time MCMC through
using a lower bounding rate, λ+ . The possibility for lowering λ+ is greater, however, for the canonical version, and
this is a second advantage it has over using sub-sampling. For example, with an additional pre-processing cost, we
could choose
λ+ =

n
X
i=1

max |∇ log πi (x)|.
x

(20)

Such a choice can be used with sub-sampling if our estimate of the rate uses non-uniform sampling of data points:


∇ log πI (x)
λ(x, v) = λ+ max 0, −v
,
maxx |∇ log πI (x)|
where we sample I from 1, . . . , n, with value i having probability proportion to maxx |∇ log πi (x)|. For the canonical
implementation we can reduce the rate further if we are able to use the actual maximum of |∇ log π(x)|, but such
a choice is not valid with sub-sampling. For our example, using (20) rather than (19) will reduce the number of
iterations required by a factor of 5.3. If we used the actual maximum of |∇ log π(x)| for the canonical version, we
would reduce the number of iterations by a factor of nearly 30 when compared to using (19). Note that for all these
options for choosing λ+ , the underlying stochastic process we are simulating is unchanged – it is just the efficiency
of the simulation algorithm that is affected.
If we compare the best implementation of continous-time MCMC with subsampling, using (20), to the best version
of the canonical implementation we get that for the same accuracy the we would need just over 25 times as many
iterations using subsampling. Each iteration would be quicker, however, as it would need access just one, out of 150,
data-points.
To see any substantial gains from using subsampling, we need to have a lower variance estimator of ∇ log π(x), using,
for example, control variates (18). To implement this we need to upper bound our estimator of the rate. This is

26

8
6

1

4

Rates

0
−1

0

−3

2

−2

∇logπ(x)

−2

0

2

4

6

8

−2

0

2

6

8

4.0

4.5

40

10

20

Rates

30

5
0

10

−5

0

−10

∇logπ(x)

4
x

50

x

2.0

2.5

3.0

3.5

4.0

4.5

5.0

2.0

2.5

3.0

3.5

x

5.0

x

Figure 5: Plots of ∇ log π(x) for the mixture example (left-hand column), and rates at which the continuous-time
MCMC algorithm will switch from a negative to a positive velocity (right-hand column). For the latter plots we
show rates for the canonical process (full lines), simple sub-sampling (dashed lines) and sub-sampling with control
variates (dotted lines). The vertical dotted line shows the value of x̂. Top row if for 150 data points, and the bottom

0

200

400

600

800

1000

6
−2

0

x_t

2

4

6
4
2
x_t
0
−2

−2

0

x_t

2

4

6

row for 1,500 data points. Plots are restricted to areas of non-negligible posterior mass.

0

200

Series x.noss

400

600

800

1000

0

30

40

800

1000

1.0
0.8
0.0

0.2

0.4

ACF

0.6

0.8
0.6
ACF
0.4
0.2
0.0
20
Lag

600

Time

1.0

1.0
0.8
0.6
ACF
0.4
0.2

10

400

Series x.cv

Time

0.0
0

200

Series x.ss

Time

0

10

20
Lag

30

40

0

10

20

30

40

Lag

Figure 6: Trace plots (top row) and autocorrelation plots (bottom row) for three implementations of continuous-time
MCMC: canonical process (left-hand column); simple sub-sampling (middle column) and subsampling with control
variates (right-hand column). Auto correlation plots are values sampled every unit time-step from the continuous
sample paths.

27

possible for this application as the absoluate value of the second derivate of log πi (x) is bounded. Assume we can
find a bound, C say, then we use a bounding rate that depends of the form
λ+ (x) = |∇ log π(x̂)| + nC|x − x̂|.

(21)

To implement the resulting algorithm we can again iterate the three steps given above. The only changes are that
in step (1) we need to simulate the inter-event time from a point process with rate
λ̃+ (s) = |∇ log π(x̂)| + nC|xt + vt s − x̂|,
and in step (3) the probability of switching the velocity is
max{0, −v[∇ log π(x̂) + n(∇ log πI (xt+s ) − ∇ log πI (x̂))]}
.
|∇ log π(x̂)| + nC|xt+s − x̂|
We can get some insight into the advantage of using control variates by calculating the expected rate of switching the
velocity for the resulting algorithm and comparing with this rate for the other two implementation. This comparison
is shown in Figure 5. We see a much lower rate of switching when we use control variates if x is close to x̂ as
compared to the simple sub-sampling approach. However the rate is actually larger if x is far from x̂. Thus we see
the importance of x̂ being close to the mode of the posterior. This picture is the same for both small, n = 150, and
larger, n = 1, 500, data sets. However for larger data sets the posterior mass close to the posterior mode increases.
As such the amount of time that algorithm will be in regions where using control variates is better will increase as
we analyse larger data sets.
In Figure 6 we see output from the algorithm using control variates for n = 150. For such a small sample size, there
appears to be little advantage in using control variates. The mixing in the tails is poor, due to the large variability
of our esimators of the switching rate when x is not close to x̂. Note that we can avoid this issue by using a hybrid
scheme that estimates the rates using control variates when |x − x̂| is small, and uses simple sub-sampling otherwise.
We see advantages from using control variates as we analyse larger data sets. A comparison of our three implementations of continuous-time MCMC is shown in Table 1. Firstly note that for the canonical implementation, the
amount of time we need to run the continuous-time MCMC for decreases with sample size. This is as described in the
scaling limits discussed in Bierkens et al. (2016). The intuition is that for larger n the posterior is more concentrated,
and thus the underlying PDP process needs less time to explore the posterior. This property is also seen if we
use subsampling with control variates. Without control variates, the actual switching rates of the underlying PDP
increase quickly with n which slows down the mixing of the algorithm, and the amount of time we need to simulate
the underlying PDP for does not change much with n.
The actual computational cost also depends on the number of iterations of the algorithm per unit-time interval.
This increases with n for all implementations. However the overall computational cost, as measured by iterations
per ESS remains roughly constant when we use control variates. As the computation cost per iterations is O(1)
we see evidence that this algorithm has a computational cost that does not increase with n. By comparison the
number of iterations per ESS appears to increase roughly linearly with n if we use subsampling without control
variates. For the canonical implementation, even using the best possible global bound on the event rate, we have the
number of iterations per ESS remaining constant but the computational cost per iterations is O(n). Thus its overall
computational cost will increase linearly with n.

28

Canonical

Canonical

Bounding Rate

(20)

min |∇ log π(x)|

(20)

(21)

t per ESS

3.3

3.3

14

22

Iterations per unit time

57

11

57

210

Iteration per ESS

190

36

800

4,600

t per ESS

0.43

0.43

8.6

2.1

Iterations per unit time

570

15

570

1,000

Iteration per ESS

245

6.4

4,900

2,100

t per ESS

0.13

0.13

9.1

0.91

Iterations per unit time

5,600

100

5,600

3,800

730

13

51,000

3,500

n = 150

n = 1, 500

n = 15, 000

Iteration per ESS
Table 1:

Subsampling Control Variate

Comparison of different implementations of continuous-time MCMC: canonical, subsampling, and sub-

sampling with control variates; and how they vary as sample size, n, increases. Both canonical and subsampling
use a global bound on the event rate to simulate possible events, we give results for canonical using both (20) and
min |∇ log π(x)| as this bound. We give estimates of the time-length the MCMC algorithm needs to be run for each
effective sample size (ESS) and the average number of iterations per unit time-interval. The product of these is then
the number of iterations needed per ESS. Remember that the cost per iteration is O(n) for canonical, and O(1) for
subsampling and subsampling with control variates.

5

Discussion

We have shown how piecewise deterministic processes can be used to derive continuous-time versions of sequential
Monte Carlo and MCMC algorithms. These algorithms are fundamentally different from more standard discrete-time
versions. Currently only a few specific algorithms, from a much wider class of possibilities, have been suggested.
Whilst we have suggested a few extensions of existing methods, these just touch the surface of the range of developments that are possible. Whilst not discussed here, the continuous-time SMC methods seem particularly well-suited
for implementation on a distributed computing architecture, as evolution of particles can be carried out in parallel.
As well as such potential methodological developments, there are a wide-range of open theoretical questions. For
example, can we get results on how well continuous-time MCMC mixes? For such results to be practically meaningful they would need to account for the computational cost of simulating the underlying PDP, as opposed to just
measuring the mixing properties of the PDP itself. Can we characterise the situations where continuous-time MCMC
is more efficient than its discrete-time counterpart? Or understand which versions of continuous-time MCMC are
most efficient, and when?
We have also shown how subsampling ideas, where we approximate the gradient of the log-posterior using a small
sample of data points, can be used with these continuous-time methods. Unlike other subsampling approaches,
the methods still remain “exact”, in the sense that they still target the true posterior. Subsampling reduces the
computational cost per iteration but does lead to a increase in Monte Carlo error for a fixed number of iterations.
In the examples we have considered, it is only when using control variate ideas to reduce the variance of our subsampling estimator of the gradient of the log posterior, that we see any overall gain in efficiency of the algorithm.
Furthermore, when using suitable control variates, it appears possible to obtain algorithms whose computational
cost per effective sample size increases sub-linearly with the number of data points. This adds to existing evidence

29

of the importance of using control variates if we wish to have some form of super-efficiency for big data problems
(Bardenet et al., 2015).
Acknowledgements The authors think the Engineering and Physical Sciences Research Council for support through
grants EP/K014463/1 (i-Like) and EP/D002060/1 (CRiSM).

References
Andrieu, C. and Roberts, G. O. (2009). The pseudo-marginal approach for efficient Monte Carlo computations.
Annals of Statistics, 37:697–725.
Bardenet, R., Doucet, A., and Holmes, C. (2015).

On Markov chain Monte Carlo methods for tall data.

ArXiv:1505.02827.
Beskos, A., Papaspiliopoulos, O., and Roberts, G. O. (2006). Retrospective exact simulation of diffusion sample
paths with applications. Bernoulli, 12(6):1077–1098.
Beskos, A., Papaspiliopoulos, O., and Roberts, G. O. (2008). A factorisation of diffusion measure and finite sample
path constructions. Methodology and Computing in Applied Probability, 10:85–104.
Beskos, A. and Roberts, G. O. (2005). Exact simulation of diffusions. The Annals of Applied Probability, 15:2422–
2444.
Bierkens, J. (2015). Non-reversible Metropolis-Hastings. Statistics and Computing, 25:1–16.
Bierkens, J., Bouchard-Cote, A., Duncan, A., Doucet, A., Fearnhead, P., Roberts, G., and Vollmer, S. (2016).
Piecewise Deterministic Markov Processes for Scalable Monte Carlo on Restricted Domains. Submitted.
Bierkens, J. and Duncan, A. (2016). Limit theorems for the Zig-Zag process. ArXiv:1607.08845.
Bierkens, J., Fearnhead, P., and Roberts, G. (2016). The Zig-Zag Process and Super-Efficient Sampling for Bayesian
Analysis of Big Data. ArXiv:1607.03188.
Bierkens, J. and Roberts, G. (2015). A piecewise deterministic scaling limit of Lifted Metropolis-Hastings in the
Curie-Weiss model. ArXiv:1509.00302.
Bouchard-Côté, A., Vollmer, S. J., and Doucet, A. (2015). The Bouncy Particle Sampler: A Non-Reversible RejectionFree Markov Chain Monte Carlo Method. ArXiv:1510.02451.
Burq, Z. and Jones, O. (2008). Simulation of Brownian motion at first passage times. Mathematics and Computers
in Simulation, 77:64–71.
Carpenter, J., Clifford, P., and Fearnhead, P. (1999). An improved particle filter for non-linear problems. IEE
proceedings-Radar, Sonar and Navigation, 146:2–7.
Cinlar, E. (2013). Introduction to stochastic processes. Courier Corporation.
Davis, M. H. (1984). Piecewise-deterministic Markov processes: A general class of non-diffusion stochastic models.
Journal of the Royal Statistical Society. Series B, 46:353–388.
Davis, M. H. A. (1993). Markov models and optimization, volume 49 of Monographs on Statistics and Applied
Probability. Chapman & Hall, London.
30

Del Moral, P. and Guionnet, A. (2001). On the stability of interacting processes with applications to filtering and
genetic algorithms. Annales de l’IHP Probabilités et statistiques, 37(2):155–194.
Diaconis, P., Holmes, S., and Neal, R. (2000). Analysis of a nonreversible Markov chain sampler. The Annals of
Applied Probability, 10(3):726–752.
Douc, R., Moulines, E., Olsson, J., et al. (2014). Long-term stability of sequential Monte Carlo methods under
verifiable conditions. The Annals of Applied Probability, 24(5):1767–1802.
Doucet, A., Godsill, S. J., and Andrieu, C. (2000). On sequential Monte Carlo sampling methods for Bayesian
filtering. Statistics and Computing, 10:197–208.
Ethier, S. N. and Kurtz, T. G. (2005). Markov Processes: Characterization and Convergence (Wiley Series in
Probability and Statistics). Wiley-Interscience.
Fearnhead, P., Latuszynski, K., Roberts, G. O., and Sermaidis, G. (2016). Continuous-time Importance Sampling:
Monte Carlo Methods which Avoid Time-Discretisation Error.
Fearnhead, P., Papaspiliopoulos, O., and Roberts, G. O. (2008). Particle filters for partially-observed diffusions.
Journal of the Royal Statistical Society Series B, 70:755–777.
Foulkes, W., Mitas, L., Needs, R., and Rajagopal, G. (2001). Quantum Monte Carlo simulations of solids. Reviews
of Modern Physics, 73(1):33.
Girolami, M. and Calderhead, B. (2011). Riemann manifold Langevin and Hamiltonian Monte Carlo methods.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214.
Gustafson, P. (1998). A guided walk Metropolis algorithm. Statistics and Computing, 8(4):357–364.
Kitagawa, G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. Journal of
Computational and Graphical Statistics, 5:1–25.
Lewis, P. A. W. and Shedler, G. S. (1979). Simulation of nonhomogeneous Poisson processes by thinning. Naval Res.
Logist. Quart., 26(3):403–413.
Liu, J. S. and Chen, R. (1995). Blind deconvolution via sequential imputations. Journal of the American Statistical
Association, 90:567–576.
Liu, J. S. and Chen, R. (1998). Sequential Monte Carlo methods for dynamic systems. Journal of the American
Statistical Association., 93:1032–1044.
Lyne, A.-M., Girolami, M., Atchad, Y., Strathmann, H., and Simpson, D. (2015). On Russian Roulette Estimates
for Bayesian Inference with Doubly-Intractable Likelihoods. Statist. Sci., 30(4):443–467.
Ma, Y.-A., Chen, T., and Fox, E. (2015). A complete recipe for stochastic gradient MCMC. In Advances in Neural
Information Processing Systems, pages 2917–2925.
McGrayne, S. B. (2011). The theory that would not die: how Bayes’ rule cracked the enigma code, hunted down
Russian submarines, & emerged triumphant from two centuries of controversy. Yale University Press.
Neal, R. M. (1998). Suppressing random walks in Markov chain Monte Carlo using ordered overrelaxation. In
Learning in graphical models, pages 205–228. Springer.

31

Neal, R. M. (2004). Improving Asymptotic Variance of MCMC Estimators: Non-reversible Chains are Better.
Technical report, No. 0406, Department of Statistics, University of Toronto.
Neal, R. M. et al. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2:113–162.
Neiswanger, W., Wang, C., and Xing, E. (2013).

Asymptotically Exact, Embarrassingly Parallel MCMC.

arXiv:1311.4780.
Peters, E. A. J. F. and de With, G. (2012). Rejection-free Monte Carlo sampling for general potentials. Physical
Review E, 85(2):026703.
Pollock, M., Fearnhead, P., Johansen, A., and Roberts, G. O. (2016a). An Unbiased and Scalable Monte Carlo
Method for Bayesian Inference for Big Data. ArXiv:1609.03436.
Pollock, M., Johansen, A. M., and Roberts, G. O. (2016b). On the exact and ε-strong simulation of (jump) diffusions.
Bernoulli, 22(2):794–856.
Quiroz, M., Villani, M., and Kohn, R. (2015). Speeding up MCMC by efficient data subsampling. Riksbank Research
Paper Series, (121).
Robert, C. and Casella, G. (2011). A Short History of Markov Chain Monte Carlo: Subjective Recollections from
Incomplete Data. Statistical Science, 26(1):102–115.
Scott, S. L., Blocker, A. W., and Bonassi, F. V. (2016). Bayes and Big Data: The Consensus Monte Carlo Algorithm.
International Journal of Management Science and Engineering Management, 11(2):78–88.
Srivastava, S., Cevher, V., Tran-Dinh, Q., and Dunson, D. B. (2015). WASP: Scalable Bayes via barycenters of
subset posteriors. In AISTATS.
Welling, M. and Teh, Y. W. (2011). Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of
the 28th International Conference on Machine Learning (ICML-11), pages 681–688.

32

