Piecewise Deterministic Markov Processes for Continuous-Time Monte Carlo

1 Introduction

• Continuous-time MCMC and SMC are both based on simulating PDPs, which evolve deterministically between random event times and may also have random transitions at events.Only a finite number of events and transitions are simulated in any fixed time interval.Source: §1 · p. 2
• Existing continuous-time methods are specific examples of a potentially much wider class of PDP-based MCMC and SMC algorithms.The paper introduces relevant PDP theory informally to make developing such methods more accessible.Source: §1 · p. 2
• For product-form Bayesian posteriors, unbiased subsampling can estimate the gradient of the log posterior within continuous-time methods while preserving the true posterior target.The paper contrasts this with cited discrete-time subsampling methods, which target an approximation to the posterior.Source: §1 · p. 2
• The paper covers continuous-time SMC and MCMC, their PDP connections and target validity, and implementation using subsampling, including the role of low-variance estimators in efficiency.The two algorithm sections are written to be readable independently.Source: §1 · pp. 2–3

2 Piecewise Deterministic Markov Processes

• A PDP is a time-homogeneous continuous-time process whose state follows deterministic dynamics between random events, with transitions at events drawn from a specified kernel.Its definition also specifies an initial state distribution.Source: §2 · p. 3

2.1 Simulating a PDP

• To simulate a PDP, iterate by simulating the next event time, evolving the state to just before that event, and drawing the post-event state from the transition kernel.For a tractable rate, event times can be simulated directly; otherwise thinning uses an upper-bound rate and accepts candidate events probabilistically.Source: §2.1 · pp. 3–4

2.2 Analysing a PDP

• The generator and Fokker–Planck equation provide tools for checking weighted-sample validity in continuous-time importance sampling and invariant distributions in continuous-time MCMC.The paper introduces these tools as they are needed to analyze the algorithms reviewed later.Source: §2.2 · p. 4

2.2.1 The Generator

• The generator gives the time derivative of the expected value of a suitable function of the process, combining deterministic evolution with expected changes at events.Its form uniquely specifies the dynamics of the underlying continuous-time Markov process.Source: §2.2.1 · p. 5

2.2.2 The Forward Operator and Fokker-Planck Equation

• The forward operator is the adjoint of the generator, and applying it to the process density yields the Fokker–Planck equation.An invariant density satisfies A* p(z) = 0, where the probability flows associated with motion and events balance.Source: §2.2.2 · pp. 5–6

3 Continuous-Time Sequential Importance Sampling

• CIS targets distributions solving a known partial differential equation, including diffusion distributions whose density is unavailable, by using a tractable proposal process and importance weights.The proposal's transition density is assumed known and simple to simulate; weighted evaluations are intended to estimate target expectations unbiasedly.Source: §3 · pp. 6–7, §3.1 · p. 7

3.1 The CIS Algorithm

• CIS simulates proposal values at random event times and updates the importance weight by an incremental factor based on the proposal and target transition-density evolution.The process can be represented as a PDP with state (Y_t, W_t, S_t), where the components track the most recent proposal value, weight, and elapsed time since an event.Source: §3.1 · p. 8, §3.1.1 · p. 8

3.1.1 CIS as a PDP

• The CIS-as-PDP representation makes validity a requirement that the expected weighted proposal density equals the target density at time t.The target density is obtained as E(W_t q_{S_t}(x|Y_t)) = p_t(x).Source: §3.1.1 · p. 8

3.1.2 Obtaining the Incremental Weight

• The incremental weight is 1 + ρ(x, y, s)/λ̃(s), with ρ reflecting the difference between target and proposal transition-density evolution.A closer-matching proposal reduces weight variability; increasing the event rate lowers incremental-weight variance but increases event count and computational cost.Source: §3.1.2 · p. 10
• For some target–proposal combinations, ρ(x, y, s) can fall below −λ̃(s), producing negative weights and potentially undermining the finite absolute-weight moments needed for the informal derivation.Larger event rates reduce the probability of negative weights, but some naive CIS implementations may not satisfy E(|W_t|) < ∞.Source: §3.1.2 · pp. 10–11

3.2 Continuous-time Sequential Monte Carlo

• Continuous-time SMC jointly runs N CIS processes, propagates and reweights each for intervals of length h, and resamples at scheduled times.Resampling can also be triggered by the effective sample size; the resampling scheme preserves each particle’s weight sign.Source: §3.2 · p. 11
• With negative weights, resampling leaves E(|Wt|) unchanged, so growth in this quantity can still degrade Monte Carlo performance over time.Consequently, the long-term stability often associated with resampled SMC is unavailable in the presence of negative weights.Source: §3.2 · p. 11
• Alternative resampling may reduce E(|Wt|) at resampling times, while sufficiently large CIS event rates can make negative weights rare over the intended run.The resampling step must be unbiased and preserve particle-weight signs.Source: §3.2 · p. 11

3.3 CIS for Big Data: the SCALE Algorithm

• SCALE can be viewed as a CIS or continuous-time SMC algorithm for a posterior assumed to factor into observation-level terms.Each factor combines one observation’s likelihood with the 1/nth power of the prior.Source: §3.3 · p. 12
• CIS targets a stochastic process with π(x) as its stationary distribution; after burn-in, weighted samples from that process approximate π(x).The paper identifies a process with π(x) invariant and interprets it as Brownian motion with killing, conditioned on survival.Source: §3.3 · p. 12
• Subsampling can unbiasedly estimate the incremental weights, and control variates reduce their variance to help control Monte Carlo error.The incremental weights depend on the posterior through log π(x), a sum that can be estimated from a data subsample.Source: §3.3 · pp. 12–13
• The original SCALE construction uses bounds on the proposal path to choose an event rate that prevents negative weights, but this can make simulation at least an order of magnitude slower.A version allowing rare negative weights is also considered.Source: §3.3 · p. 13

3.4 Extensions

• The extensions illustrate a wider class of potential continuous-time Monte Carlo methods based on piecewise-deterministic processes.The paper presents two examples: alternative CIS proposals and alternative processes for SCALE.Source: §3.4 · p. 13

3.4.1 Alternative Proposal Distributions

• CIS can use alternative differentiable transition-density families, including t-distributed transitions with suitably chosen degrees of freedom.The family must permit calculation of the incremental weight.Source: §3.4.1 · pp. 13–14

3.4.2 Alternative SCALE Processes

• An alternative SCALE process can target the Langevin diffusion, whose posterior invariance is retained when used in CIS.Its incremental weight is easier to estimate than SCALE’s because it does not involve a squared log-posterior gradient, but a term with variance 1/s requires care for small inter-event times.Source: §3.4.2 · p. 14

3.5 Example: Inference for Mixture Models

• The example’s ρ has multiple pronounced modes despite a relatively well-behaved posterior and log-posterior, creating potential mixing problems.With particles initialized uniformly on [−10, −5], SCALE appears stuck near a local mode of ρ; prior initialization appears to converge by about time 15.Source: §3.5 · pp. 14–15
• Control variates substantially lower the variance of the ρ estimate near the posterior mode, but variance can become worse than without them farther away.When the reference point is close to the posterior mode, the proportionate variance reduction is O(n) with high probability.Source: §3.5 · p. 17
• Subsampling without control variates performs increasingly poorly as n grows, while moving the control-variate reference point away from the mode worsens cost and variance.Even with the reference point one or three posterior standard deviations from the mode, cost remains constant in n at fixed Var(Wh).Source: §3.5 · pp. 17–18
• At fixed Var(Wh), subsampling with control variates keeps data-access cost constant as n increases, while no subsampling requires linearly increasing cost.Figure 4 compares computational cost in plot (b) and estimated variance in plot (c).Source: §3.5 · p. 17, Figure 4 · p. 17

• Continuous-time MCMC simulates a piecewise-deterministic process designed to have the target distribution as its stationary distribution.These methods are non-reversible and can use subsampling ideas, making them promising for big-data applications.Source: §4 · p. 18
• Non-reversible MCMC is described as suppressing reversible MCMC’s random-walk behavior and potentially exploring the state space more rapidly.The passage notes substantial evidence that non-reversible MCMC can be more efficient than standard reversible MCMC.Source: §4 · p. 18

• A discrete-time proposal moves along a unit-speed velocity and flips it; separate velocity updates are needed for irreducibility when the position dimension exceeds one.The algorithm targets the joint distribution of position and velocity, with uniform velocity over directions of unit speed.Source: §4.1 · p. 18
• The paired discrete-time moves satisfy detailed balance, but the position proposal alone moves only along the current velocity and is reducible in dimensions above one.A velocity update with the uniform velocity distribution as its stationary distribution supplies the additional move.Source: §4.1 · p. 18
• As the step size tends to zero, proposal rejections become Poisson events with rate max{0, −v_t · ∇π(x_t)}, while position follows constant-velocity motion between events.Events occur only when the process moves toward lower probability mass; the derivation assumes, for example, that π(x) is twice differentiable.Source: §4.1 · p. 19
• The limiting process is a piecewise-deterministic process whose velocity changes at event times; more general event rates and velocity transitions can preserve π(x) as the marginal stationary distribution.The general formulation specifies deterministic dynamics, an event rate, and a transition distribution at events.Source: §4.1 · p. 19

4.2 The Stationary Distribution of the PDP

• The adjoint-generator condition links event rates and event transitions to an invariant distribution whose position marginal is π(x).For velocity independent of position, the rate and transition choices must balance probability flow against the deterministic dynamics; the velocity distribution must have zero mean.Source: §4.2 · pp. 19–20
• The canonical rates are the smallest rates satisfying the condition, λ(x, v) = max{0, −v · ∇x log π(x)}; in one dimension, their use minimizes estimator asymptotic variance.The rates depend on the target through its log-gradient, so the target need only be known up to proportionality.Source: §4.2 · p. 20

4.2.1 Different continuous-time MCMC algorithms

• The three described samplers coincide for one-dimensional targets but differ in higher dimensions, and the section assumes an unbounded target domain.Bounded-domain extensions are cited to Bierkens et al. (2016).Source: §4.2.1 · p. 20

Pure Reflection and Refresh

• Pure reflection uses canonical rates but is reducible for multidimensional targets because it stays on the line set by its initial velocity; velocity refreshment can address this.Refresh times may be deterministic or random, and the refresh distribution must have the desired velocity distribution as its stationary distribution.Source: Pure Reflection and Refresh · pp. 20–21

Bouncy Particle Sampler

• The Bouncy Particle Sampler reflects the velocity component along the target log-gradient, leaving orthogonal components unchanged, and adds independent Poisson refreshment.Bouchard-Cˆot´e et al. (2015) report that some targets make the sampler without refresh reducible, while any non-zero refresh rate gives the target as its asymptotic distribution.Source: Bouncy Particle Sampler · p. 21

Zig-Zag Sampler

• The Zig-Zag sampler uses d coordinate-wise velocity flips, each with its own event rate; the rates need not be uniquely determined.Rigorous ergodicity is established for one-dimensional targets, factorized targets, or switching rates positive everywhere; experiments suggest broader ergodicity.Source: Zig-Zag Sampler · p. 21
• The Zig-Zag framework can allow non-orthogonal velocity directions or more directions than dimensions, but the benefits are unclear.Source: Zig-Zag Sampler · p. 21

Alternatives

• An alternative draws a new unit velocity at each event using canonical rates and retains π(x) as the position marginal of an invariant distribution.The framework also allows velocity’s invariant distribution to depend on position, provided rates and transitions satisfy the corresponding condition.Source: Alternatives · p. 22

• After burn-in, the PDP skeleton consists of event times and states over the retained interval, which support Monte Carlo estimates of expectations under π(x).Source: §4.3 · p. 22
• A path-based estimate averages the function along the trajectory by integrating over each interval between consecutive events.Between events, x_t follows the linear path from the interval’s initial state and velocity.Source: §4.3 · p. 22
• When those integrals are difficult to evaluate, a Monte Carlo approximation samples equally spaced times across the retained interval and evaluates the process state there.The states at those times can be calculated from the event times and event states.Source: §4.3 · p. 22

4.4 Exact Approximation versions and Subsampling

• Exact-approximate versions of continuous-time MCMC are possible when target-distribution estimators are positive and unbiased.The section illustrates the approach with Pure Reflection and extends the ideas to Zig-Zag and the Bouncy Particle Sampler.Source: §4.4 · pp. 22–23, §4.4.1 · p. 23

4.4.1 Exact Approximation for Pure Reflection and Zig Zag

• An unbiased estimator of −∇log π(x) lets the Pure Reflection process retain π(x) as its marginal invariant distribution.In general, the estimator increases event rates and thinning bounds, which can reduce mixing and increase simulation cost; these losses are expected to grow with estimator variability.Source: §4.4.1 · pp. 23–24

4.4.2 Exact Approximation for the Bouncy Particle Sampler

• For the Bouncy Particle Sampler, using the same gradient-estimator realization for event acceptance and the velocity flip yields invariant distribution π(x)pv(v).Here pv(v) is uniform over velocities of fixed speed.Source: §4.4.2 · p. 24

4.4.3 Use of Subsampling

• Subsampling can estimate −∇log π(x) from data points, and reducing estimator variance is important for sampler efficiency.The simplest estimator samples an index uniformly from 1 to n; the section also discusses control variates.Source: §4.4.3 · p. 24
• A control-variate subsampling estimator near a posterior mode is expected to have variance growing linearly with n, versus O(n^2) for the simple estimator.The stated expectation applies when x is within O(√n) of the chosen near-mode value and π(x) is sufficiently smooth.Source: §4.4.3 · pp. 24–25

4.5 Example: Mixture Model

• The example compares canonical, simple-subsampling, and control-variate implementations; each proposal event advances x linearly and may switch velocity according to the actual-to-bound rate.The global event-rate bound generates putative event times; simple subsampling evaluates the rate using a sampled data point.Source: §4.5 · p. 25
• With n=150, subsampling switches velocity more frequently; canonical MCMC suppresses random-walk behavior better, mixes better, and has lower autocorrelation.The comparison is shown in Figure 6's trace and autocorrelation plots.Source: §4.5 · p. 26, Figure 6 · p. 27
• Using (20) reduces iterations 5.3-fold; using the canonical actual maximum reduces them nearly 30-fold versus (19), while λ+ choices leave the underlying process unchanged.The actual-maximum bound is not valid with subsampling.Source: §4.5 · p. 26
• Just over 25 times as many iterations are needed with the best subsampling implementation as with the best canonical implementation for equal accuracy, though each subsampling iteration accesses one of 150 data points.Source: §4.5 · p. 26
• Control variates lower switching rates near ˆx but raise them farther away; for n=150, tail mixing is poor, while a hybrid scheme can use simple subsampling in the tails.For larger data sets, posterior mass near the mode increases, expanding the regions where control variates have lower rates than simple subsampling.Source: §4.5 · p. 28
• With control variates, iterations per ESS remain roughly constant as n grows, while they rise roughly linearly without control variates and canonical computational cost rises linearly.The control-variate implementation has O(1) cost per iteration; canonical iterations cost O(n).Source: §4.5 · p. 28

5 Discussion

• Continuous-time SMC methods may suit distributed computing because particle evolution can be carried out in parallel.The discussion also identifies a much wider class of potential methods beyond the few algorithms suggested to date.Source: §5 · p. 29
• Open questions include how to assess MCMC mixing while accounting for simulation cost, when continuous-time MCMC is more efficient, and which versions work best.The authors note that mixing results would need to account for the computational cost of simulating the underlying process.Source: §5 · p. 29
• Subsampling preserves exact posterior targeting but increases Monte Carlo error at a fixed iteration count; observed efficiency gains required control variates.With suitable control variates, computational cost per effective sample size appears able to grow sub-linearly with data size.Source: §5 · p. 29
• The observed need for control variates for big-data efficiency aligns with prior evidence emphasizing their importance for super-efficiency (Bardenet et al., 2015).The authors report that their examples showed overall efficiency gains only when control variates reduced subsampling-estimator variance.Source: §5 · pp. 29–30