Zig-Zag Sampler: A PDMP MCMC Technique

• Zig-Zag sampler is a PDMP method that uses an augmented state space combining position and discrete velocity to achieve exact continuous-time sampling.
• It simulates deterministic linear motion between random event times where coordinate-specific velocity flips are triggered by gradient information, eliminating discretization error.
• Enhancements like sub-sampling, control variates, and extensions for complex target geometries improve scalability and efficiency in large-scale Bayesian computation.

Searching arXiv for core and recent Zig-Zag sampler papers to ground the article. Searching for “Zig-Zag sampler” on arXiv. The Zig-Zag sampler is a non-reversible, continuous-time Markov chain Monte Carlo method based on a piecewise deterministic Markov process (PDMP) on an augmented state space that combines a position variable with a discrete velocity vector. For targets on $\mathbb{R}^d$ of the form $\pi(x)\propto \exp\{-U(x)\}$, the process evolves by linear motion between random event times and flips one velocity coordinate at a time with rates determined by the gradient of the potential. This construction yields exact continuous-time dynamics, no time-discretisation error, and an invariant marginal equal to the desired target distribution. In the PDMP literature, Zig-Zag and the Bouncy Particle Sampler are the two original versions on which many later samplers are based (Bierkens et al., 2016, Wu et al., 2018).

1. Core PDMP formulation

For the standard Zig-Zag sampler, the target density is written as

$\pi(x)=\frac{1}{Z}\exp(-U(x)),$

with $U:\mathbb{R}^d\to\mathbb{R}$ differentiable, and the augmented state space is

$E=\mathbb{R}^d\times\{-1,1\}^d.$

A state is $z_t=(x_t,v_t)$, where $x_t\in\mathbb{R}^d$ is position and $v_t\in\{-1,1\}^d$ is velocity. Between event times, the dynamics are deterministic: $\frac{dx_t}{dt}=v_t,\qquad \frac{dv_t}{dt}=0.$ Hence the trajectory is piecewise linear in position and piecewise constant in velocity (Sutton et al., 2022).

The event structure is coordinate-wise. For $i=1,\dots,d$, the canonical switching rate is

$\pi(x)\propto \exp\{-U(x)\}$0

equivalently $\pi(x)\propto \exp\{-U(x)\}$1. At an event of type $\pi(x)\propto \exp\{-U(x)\}$2, only the $\pi(x)\propto \exp\{-U(x)\}$3-th velocity component flips sign. In generator form,

$\pi(x)\propto \exp\{-U(x)\}$4

where $\pi(x)\propto \exp\{-U(x)\}$5 flips the $\pi(x)\propto \exp\{-U(x)\}$6-th velocity coordinate (Bierkens et al., 2021, Chevallier et al., 2021).

A more general admissible family of rates is

$\pi(x)\propto \exp\{-U(x)\}$7

with $\pi(x)\propto \exp\{-U(x)\}$8. The key balance relation is

$\pi(x)\propto \exp\{-U(x)\}$9

which ensures that the augmented invariant law has density proportional to $\pi(x)=\frac{1}{Z}\exp(-U(x)),$0 in position and uniform law on $\pi(x)=\frac{1}{Z}\exp(-U(x)),$1 in velocity (Bierkens et al., 2016). The process is therefore exact in the sense that the $\pi(x)=\frac{1}{Z}\exp(-U(x)),$2-marginal of stationarity is the desired target distribution.

2. Event-driven simulation and asymptotic theory

Simulation proceeds by following the deterministic flow until the next event. In the standard skeleton construction, for each coordinate $\pi(x)=\frac{1}{Z}\exp(-U(x)),$3, one simulates the next event time from the inhomogeneous Poisson rate $\pi(x)=\frac{1}{Z}\exp(-U(x)),$4, takes the minimum over coordinates, updates the position by linear motion, and flips the corresponding velocity component. The continuous trajectory is then reconstructed between successive events (Sutton et al., 2022).

A central practical consequence of the PDMP formulation is that Zig-Zag can be simulated without time-discretisation error. Event times may be sampled directly when the integrated hazard is tractable, or by Poisson thinning using explicit upper bounds on the rates. This exact event-driven character distinguishes Zig-Zag from Hamiltonian Monte Carlo, which uses a discrete approximation of Hamiltonian dynamics and a Metropolis correction, and from standard Metropolis–Hastings, which is discrete-time and reversible (Bierkens et al., 2016).

The one-dimensional process admits explicit limit theorems. With state $\pi(x)=\frac{1}{Z}\exp(-U(x)),$5 and generator

$\pi(x)=\frac{1}{Z}\exp(-U(x)),$6

where

$\pi(x)=\frac{1}{Z}\exp(-U(x)),$7

central limit theorems hold for ergodic averages under Lyapunov-type conditions, and the asymptotic variance can be characterized explicitly. The same analysis also yields a diffusion limit when the switching rate tends to infinity: after rescaling time, the position process converges to an Itô diffusion whose drift and diffusion coefficients depend on $\pi(x)=\frac{1}{Z}\exp(-U(x)),$8 and $\pi(x)=\frac{1}{Z}\exp(-U(x)),$9 (Bierkens et al., 2016).

This asymptotic theory makes clear that switching intensity is not merely a technical device. Excess switching can preserve invariance, but it also changes diffusivity and asymptotic variance. A plausible implication is that practical performance depends on balancing direction persistence against the cost and variance inflation induced by frequent velocity flips (Bierkens et al., 2016).

3. Bayesian computation, sub-sampling, and large-sample scaling

In Bayesian applications with i.i.d. data, the posterior potential typically decomposes as

$U:\mathbb{R}^d\to\mathbb{R}$0

so $U:\mathbb{R}^d\to\mathbb{R}$1 is a sum of data contributions. Canonical Zig-Zag evaluates this full gradient and therefore has $U:\mathbb{R}^d\to\mathbb{R}$2 cost per switching-rate evaluation. The multi-dimensional Zig-Zag paper introduced a sub-sampling construction in which one replaces the full gradient by an unbiased estimator inside the switching rule, while preserving the exact posterior as invariant distribution; this was termed an exact approximate scheme (Bierkens et al., 2016).

The same work introduced control variates centered at a reference point $U:\mathbb{R}^d\to\mathbb{R}$3, typically near a posterior mode. After an initial pre-processing step, the resulting Zig-Zag process can be super-efficient: essentially independent samples from the posterior are obtained at a computational cost which does not depend on the size of the data (Bierkens et al., 2016). The large-sample scaling analysis later made this statement precise. In the transient phase, Zig-Zag trajectories are approximated by an ODE whose drift points in the direction of decreasing KL-divergence between the assumed model and the true distribution. In the stationary phase, different Zig-Zag variants admit weak limits after appropriate rescaling, and the analysis estimates that for data size $U:\mathbb{R}^d\to\mathbb{R}$4, suitably controlled Zig-Zag with sub-sampling has $U:\mathbb{R}^d\to\mathbb{R}$5 cost for an essentially independent sample, representing an $U:\mathbb{R}^d\to\mathbb{R}$6 speed-up over the canonical version (Agrawal et al., 2024).

The same scaling analysis also clarifies a limitation. Sub-sampling without effective control variates increases switching rates through Jensen’s inequality and can slow convergence and poor mixing of the process, even though each event is cheap (Agrawal et al., 2024). This sharpens an important methodological point: the scalability of Zig-Zag in big-data regimes is not a generic consequence of sub-sampling alone, but depends on variance reduction and on how the transient and stationary regimes interact.

4. Variants for complex target geometries and state spaces

A major line of development extends Zig-Zag beyond smooth densities on Euclidean space. For piecewise-smooth densities and bounded domains, the main difficulty is defining the process when the deterministic trajectory hits a discontinuity. A general boundary condition can be stated in terms of a boundary kernel that preserves the correct flux across discontinuity surfaces; with this condition, one can construct Zig-Zag samplers for densities with discontinuities or compact support, including truncated targets (Chevallier et al., 2021).

For sparse Bayesian models with atomic mass at zero, Sticky Zig-Zag augments the state space so that each coordinate has two copies of zero, $U:\mathbb{R}^d\to\mathbb{R}$7 and $U:\mathbb{R}^d\to\mathbb{R}$8, and introduces sticky boundaries on coordinate hyperplanes. When a coordinate hits zero with velocity pointing inward, it freezes there and later thaws at rate $U:\mathbb{R}^d\to\mathbb{R}$9. The resulting invariant measure has the mixed continuous–atomic form

$E=\mathbb{R}^d\times\{-1,1\}^d.$0

which is appropriate for spike-and-slab posteriors and other sparse targets (Bierkens et al., 2021).

The Zig-Zag construction has also been extended to hybrid discrete/continuous state spaces. In phylogenetic applications based on the Kingman coalescent, the discrete component is a ranked tree topology and the continuous component contains branch lengths and mutation parameters. Between boundary hits, the process follows ordinary Zig-Zag dynamics in the continuous coordinates; at a boundary, a kernel $E=\mathbb{R}^d\times\{-1,1\}^d.$1 performs a topology change together with the corresponding velocity update. This produces a non-reversible sampler on tree space that avoids the analytically and implementationally complex boundary crossings encountered by hybrid-space HMC (Koskela, 2020).

In function-space Bayesian inference, an infinite-dimensional Zig-Zag sampler has been constructed on $E=\mathbb{R}^d\times\{-1,1\}^d.$2, where $E=\mathbb{R}^d\times\{-1,1\}^d.$3 is a separable Hilbert space and $E=\mathbb{R}^d\times\{-1,1\}^d.$4 is a smoother velocity space. Its generator is

$E=\mathbb{R}^d\times\{-1,1\}^d.$5

with

$E=\mathbb{R}^d\times\{-1,1\}^d.$6

Under suitable assumptions on $E=\mathbb{R}^d\times\{-1,1\}^d.$7, the process exists, is non-explosive, and preserves the target measure, although the paper does not prove an explicit infinite-dimensional ergodicity theorem for Zig-Zag (Dobson et al., 2022).

5. Tempering, Gibbs blocking, and related PDMP constructions

Continuous tempering provides one route to overcoming the tendency of gradient-based samplers to become trapped in local modes or to explore heavy tails slowly. In the continuously tempered Zig-Zag sampler, the state is augmented with an inverse temperature $E=\mathbb{R}^d\times\{-1,1\}^d.$8 and an additional velocity component. The extended target interpolates between an easy base density $E=\mathbb{R}^d\times\{-1,1\}^d.$9 at $z_t=(x_t,v_t)$0 and the posterior $z_t=(x_t,v_t)$1 at $z_t=(x_t,v_t)$2, with a marginal over $z_t=(x_t,v_t)$3 that consists of a continuous component on $z_t=(x_t,v_t)$4 and a point mass at $z_t=(x_t,v_t)$5. When the trajectory first hits $z_t=(x_t,v_t)$6, the $z_t=(x_t,v_t)$7-velocity is set to zero, the process remains on the posterior for an exponential holding time with rate

$z_t=(x_t,v_t)$8

and then re-enters the tempered region with negative $z_t=(x_t,v_t)$9-velocity. The resulting process is $x_t\in\mathbb{R}^d$0-ergodic and yields exact posterior samples whenever $x_t\in\mathbb{R}^d$1 (Sutton et al., 2022).

Another generalization is the Gibbs zig-zag sampler. Here the parameter is partitioned as $x_t\in\mathbb{R}^d$2, with a Zig-Zag PDMP applied to $x_t\in\mathbb{R}^d$3 conditional on $x_t\in\mathbb{R}^d$4, and Gibbs-style or Metropolis–Hastings updates applied to $x_t\in\mathbb{R}^d$5 at Poisson event times of rate $x_t\in\mathbb{R}^d$6. If the Zig-Zag part preserves $x_t\in\mathbb{R}^d$7 and the kernel $x_t\in\mathbb{R}^d$8 preserves $x_t\in\mathbb{R}^d$9, then the combined generator $v_t\in\{-1,1\}^d$0 preserves $v_t\in\{-1,1\}^d$1. Under positivity and smoothness assumptions, the process is ergodic and satisfies a central limit theorem for time averages (Sachs et al., 2020).

The Coordinate Sampler is a closely related PDMP explicitly presented as a variant of the Zig-Zag sampler. It replaces the velocity space $v_t\in\{-1,1\}^d$2 by $v_t\in\{-1,1\}^d$3, so only one coordinate moves between events. The paper proves theoretical validation and geometric ergodicity for a class of targets, and reports that the induced chain is more efficient than Zig-Zag in terms of effective sample size in several numerical examples (Wu et al., 2018). This comparison is important because it shows that Zig-Zag is not the unique ballistic PDMP construction; coordinate-wise alternatives may dominate in some regimes.

6. Empirical behavior, transient efficiency, and limitations

Empirically, Zig-Zag and its extensions have been effective in settings where reversible chains mix slowly. Continuously tempered Zig-Zag reduces root-mean-square errors and improves thinning efficiency on multimodal Gaussian mixtures, spike-and-slab transdimensional examples, and a continuous relaxation of a Boltzmann machine, while preserving direct posterior sampling at $v_t\in\{-1,1\}^d$4 (Sutton et al., 2022). In sequential Monte Carlo, a discretized Zig-Zag refinement step embedded in a composite MH kernel improves mean squared error and acceptance ratios in high-dimensional filtering examples, although that construction is explicitly a discretization rather than the exact continuous-time sampler (Han et al., 2021).

Recent transient-regime analysis further refines the picture. For convex targets under steep scaling, early Zig-Zag dynamics admit a fluid-limit description. For Gaussian targets, the transient cost measured by expected jump count is comparable to random-walk Metropolis; for convex heavy-tailed families with subquadratic growth, PDMP methods can be more efficient. For Zig-Zag specifically, under a diagonal-dominance condition, the transient choice of direction on tangency coordinates coincides with the solution of a box-constrained quadratic program involving the Hessian, while outside that regime the process can exhibit snapping, that is, many rapid flips near tangency surfaces (Agrawal et al., 19 Sep 2025). This directly qualifies a common simplification: Zig-Zag is not uniformly superior to reversible MCMC, and its transient efficiency depends strongly on target geometry.

Several limitations recur across the literature. Choice of control variates or tempering weights is problem-specific; poor calibration can produce too much time at low temperature or too much excess switching (Sutton et al., 2022, Agrawal et al., 2024). Infinite-dimensional Zig-Zag requires a smoother velocity space $v_t\in\{-1,1\}^d$5, and the paper notes that different choices of velocity magnitudes can lead to different invariant measures on $v_t\in\{-1,1\}^d$6 unless the support of the velocity law is restricted (Dobson et al., 2022). High-dimensional scaling theory comparable to that available for HMC remains incomplete in a number of settings (Sutton et al., 2022).

Taken together, these developments place the Zig-Zag sampler within a broad research program on non-reversible PDMP Monte Carlo. Its defining mechanism remains simple—piecewise linear motion with gradient-based coordinate flips—but that mechanism has supported exact sub-sampling for big data, tempered exploration of multimodal posteriors, sparse and transdimensional extensions, hybrid discrete/continuous inference, and function-space formulations. The resulting literature presents Zig-Zag less as a single algorithm than as a family of velocity-jump constructions for target distributions whose geometry can be exploited directly in continuous time (Bierkens et al., 2016, Dobson et al., 2022).