arXiv CSopen access

Windowed thinning and query complexity for the bouncy particle and Zigzag samplers

WINDOWED THINNING AND QUERY COMPLEXITY FOR THE BOUNCY PARTICLE AND ZIGZAG SAMPLERS JIANFENG LU AND YINCHEN LUO

A BSTRACT. Let µ( dx) ∝ e −U (x) dx on Rd , where U is m-strongly convex and L-smooth, and denote by κ = L/m the condition number. We consider windowed thinning, an exact simulation method for the bouncy particle sampler and the coordinate Zigzag process. The method divides a trajectory into deterministic windows and uses a gradient evaluation at the beginning of each window to

arXiv:2607.28413v1 [math.NA] 30 Jul 2026

construct a tractable local envelope for the event rate. Combining this construction with quantitative mixing estimates and finite-time bounds on the expected numbers of bounces and flips yields query complexity guarantees from a Gaussian cold start. For total-variation error ε, the expected query counts are O(κ1/2 d (d log κ + log 1ε )) gradient queries for the bouncy particle sampler and

O(κd 1/4 (d log κ + log 1ε )) full-gradient equivalents for Zigzag, where d coordinate-partial queries count as one equivalent.

1. I NTRODUCTION The bouncy particle sampler (BPS) and the Zigzag sampler offer an event-driven geometric approach to sampling. A particle follows a simple deterministic path and changes its velocity only at random event times. If those event times are generated exactly, the resulting continuoustime trajectory incurs no time-discretization error. This attractive feature leaves a less visible computational question: how many evaluations of the target gradient are needed to run the process until it is close to equilibrium? Answering this question requires controlling both the time to convergence and the cost of simulating the trajectory over that time. Our goal is to provide such a query complexity guarantee from an explicit Gaussian cold start. We consider the log-concave target measure µ( dx) = Z −1 e −U (x) dx,

x ∈ Rd ,

where Z is the normalizing constant and U ∈ C 2 (Rd ) is m-strongly convex and L-smooth: mI d ⪯ ∇2U (x) ⪯ LI d ,

(1)

x ∈ Rd .

We assume that m, L, and the minimizer x ⋆ of U are available, together with oracle access to either the full gradient ∇U (x) or an individual coordinate derivative ∂i U (x). Date: July 2026. 2020 Mathematics Subject Classification. 65C05, 60J25, 65C40, 65Y20. Key words and phrases. Piecewise deterministic Markov process, bouncy particle sampler, Zigzag process, Poisson thinning, query complexity, hypocoercivity. This work is supported in part by the National Science Foundation under grant DMS-2309378. We thank Sinho Chewi and Lihan Wang for helpful discussions. 1

2

JIANFENG LU AND YINCHEN LUO

Within this setting, we focus on two samplers based on piecewise deterministic Markov processes (PDMPs): BPS [BCVD18] and Zigzag [BFR19]. BPS reverses the component of its velocity in the local gradient direction, whereas Zigzag flips a single velocity coordinate at each event. Both processes evolve on Rd × Rd and preserve the product measure ρ ∞ ( dx dv) = µ( dx)ϕ( dv),

ϕ := N (0, I d ).

See Section 2 for a more detailed description. We initialize them from an explicit Gaussian distribution centered at x ⋆ , i.e., from a cold start. Turning either continuous-time process into an exact algorithm requires generating its bounce or flip times. Their event rates depend on ∇U along the evolving trajectory and cannot be evaluated continuously in the oracle model. Poisson thinning [LS79] replaces the true event rate by a tractable upper envelope: candidate events are proposed from the envelope and accepted with probability equal to the ratio of the true rate to the envelope. A valid envelope preserves exactness, but a loose one can generate many rejected proposals and therefore many unnecessary queries. Our main algorithmic contribution, windowed thinning, resolves this tradeoff locally. We divide the simulation horizon into deterministic windows and evaluate the gradient at the current position at the beginning of each window. The L-Lipschitz continuity of ∇U then supplies a local envelope for the event rate in terms of the cumulative distance travelled from this anchor. Short windows keep the envelope tight but require frequent anchor queries; long windows require fewer anchors but generate more rejected proposals. Balancing these two query counts leads to the window lengths used in our bounds. Combined with quantitative mixing estimates, windowed thinning gives the following expected query counts: ¡ ¡ ¢¢ BPS: O κ1/2 d d log κ + log 1ε , ¢¢ ¡ ¡ Zigzag: O κd 1/4 d log κ + log 1ε . Here κ := L/m is the condition number and ε denotes the desired total-variation accuracy. The BPS bound counts full-gradient queries, while the Zigzag bound is stated in full-gradient equivalents, with d coordinate-partial queries counted as one equivalent. Both algorithms simulate their underlying continuous-time processes exactly; the accuracy parameter only determines the simulation horizon. Table 1 compares our bounds with those of other high-accuracy samplers whose dependence on the target accuracy is polylogarithmic. All entries include the cost of starting from a feasible cold distribution; for the Metropolis-adjusted Langevin algorithm (MALA), this includes an algorithmic warm-start phase. The oracle and guarantee conventions are recorded in the table caption. Compared with MALA and first-order rejection sampling (FORS), while our bound for BPS has a smaller κ dependence, its dimension dependence is much worse. We also note that since Zigzag is a coordinate algorithm, it is possible to reduce the complexity under additional structural assumptions, such as assuming coordinate Lipschitz conditions with better constants; we do not pursue those in this work.

WINDOWED THINNING FOR BPS AND ZIGZAG

3

TABLE 1. Cold-start query complexities of high-accuracy samplers for mstrongly convex, L-smooth targets, stated for total-variation error at most ε. All entries include a feasible initialization. The MALA entry counts first-order queries returning both U and ∇U ; the FORS and PDMP entries count expected gradient queries, with Zigzag counts expressed in full-gradient equivalents. The e suppresses logarithmic factors in d , κ, and ε−1 . notation O Method

Oracle-query complexity

MALA with algorithmic warm start [AC24]

1/2 ) e O(κd 1/2 e O(κd )

Proximal sampler with FORS [CCDR26] Windowed-thinning BPS, this work Windowed-thinning Zigzag, this work

e 1/2 d 2 ) O(κ 5/4 ) e O(κd

To establish these query counts, we analyze convergence and simulation separately. Quantitative χ2 -contraction estimates, recalled in Section 2.3, first determine how long the process must run. Over this horizon, the thinning analysis counts the anchor evaluations and proposed events. The main difficulty is to control the event count directly from the nonstationary Gaussian cold start, without assuming stationarity or a warm start. Rather than tracking individual bounces and flips, we bound their expected total number through the time-integrated squared event rates. For each sampler, a Dynkin identity for a suitable observable bounds these integrals in terms of position- and velocity-dependent quantities, which are controlled using the convexity and smoothness of the potential together with the cold-start moment estimates. This turns a pathdependent counting problem into a moment estimate. The velocity-moment bounds also control the distance travelled within each window, and hence the expected number of rejected proposals. Balancing this rejection cost against the number of anchor evaluations determines the window length. Together, these estimates yield dimension-explicit query bounds over the full simulation horizon. The details are in Section 4. Related work. The event-driven bounce mechanism underlying BPS originated with Peters and de With [PdW12] and was developed into the BPS framework by Bouchard-Côté, Vollmer and Doucet [BCVD18]. Zigzag first appeared as a scaling limit of lifted Metropolis–Hastings in the Curie–Weiss model [BR17] and was later developed as a general-purpose sampler by Bierkens, Fearnhead and Roberts [BFR19]. Both belong to the broader PDMP Monte Carlo family surveyed in [FBPR18]. Later examples include the Coordinate Sampler, the Boomerang Sampler, and Forward Event-Chain Monte Carlo [WR20, BGKR20, MDS20]. Exact trajectory simulation requires sampling event times from inhomogeneous Poisson processes. Classical thinning [LS79] and adaptive thinning [Oga81] underlie the local construction used in the original BPS work [BCVD18]. Subsequent work develops analytic and numerical envelopes [SF23, CSR22, AK24]. Our work builds on this literature with a quantitative cold-start analysis.

4

JIANFENG LU AND YINCHEN LUO

Harris-type methods establish broad ergodicity results for these samplers. For BPS, Deligiannidis, Bouchard-Côté and Doucet proved geometric ergodicity and a central limit theorem under explicit tail and curvature conditions [DBCD19], while Durmus, Guillin and Monmarché treated broader classes of targets and velocity laws using coupling and quantitative minorization [DGM20]. For the Zigzag process, Bierkens, Roberts and Zitt proved total-variation convergence, and also established exponential convergence and a central limit theorem for empirical average under stronger tail assumptions [BRZ19]. From the scaling-limit perspective, BPS and Zigzag were analyzed in the limit of d → ∞ by Bierkens, Kamatani and Roberts [BKR22, BKR25]. Agrawal, Bierkens, Kamatani and Roberts study the transient motion from low-probability regions toward a typical set using fluid limits, with expected event counts as a proxy for computational cost [ABKR25]. Explicit nonasymptotic results quantify convergence, finite-time fluctuations, and query counts. Andrieu, Durmus, Nüsken and Roussel developed a general L 2 -hypocoercivity framework covering BPS and Zigzag [ADNR21]. Lu and Wang derived explicit dimension- and condition-numberdependent rates [LW22b], based on the space-time Poincare inequality approach [CLW23]; these are the continuous-time inputs used here. Their Zigzag complexity analysis combines such a rate with a global event envelope to obtain a high-probability query bound from a warm start [LW22a]. Beyond contraction rates, Birrell and Rey-Bellet derived Bernstein-type concentration inequalities for finite-time ergodic averages and uncertainty-quantification bounds [BRB25]. From the lifting viewpoint, Eberle and Lörler interpret BPS and Zigzag as nonreversible lifts and relate their relaxation times to those of the underlying reversible dynamics [EL26]. Recent works also consider convergence on the entropy level. Under convexity, a spatial log– Sobolev inequality and suitable growth conditions, Lu proves sharp exponential relative-entropy decay for underdamped Langevin dynamics from every initial law with finite relative entropy [Lu26]. Li and Lu establish a space–time log–Sobolev inequality and hypocoercive hypercontractivity, yielding sharp decay of Rényi divergences [LL26]. Unfortunately, for BPS and Zigzag with positive refreshment and compactly supported or standard Gaussian velocity laws, Monmarché and Wang show that, even for a standard Gaussian target, no nontrivial relative-entropy contraction holds [MW26], thus our analysis is based on L 2 -hypocoercivity instead. For comparison beyond PDMP samplers, Wu, Schmidler and Chen proved the optimal warmp 1/2 e e start complexity O(κd ) for MALA [WSC22] (the earlier bound was O(κd + κ3/2 d ) [DCWY19]). Altschuler and Chewi [AC24] obtain the same complexity from a feasible start by using discretized underdamped Langevin dynamics to construct a Rényi warm start before running MALA. Chen, Chewi, Daskalakis and Rakhlin [CCDR26] implement the proximal sampler using gradient-only FORS and match this expected complexity for the same target class. With polynomial dependence on the target accuracy, Altschuler, Chewi and Zhang [ACZ26] analyze a randomized-midpoint discretization of underdamped Langevin dynamics, obtaining e 5/6 d 5/3 ε−2/3 ) gradient evaluations from a Gaussian cold start. Combining their discretization O(κ e 5/6 d 1/3 ε−2/3 ) complexity analysis with the space–time log–Sobolev inequality [LL26] yields O(κ under additional regularity assumptions [Che26, Theorem 5.3.17].

WINDOWED THINNING FOR BPS AND ZIGZAG

5

Finally, we remark that using the hypocoercive hypercontractivity and the space–time log– e q (κ2/3 d 1/3 ) Sobolev inequality [Lu26, LL26], it is possible to design samplers with query complexity O to achieve ε2 error in Renyi divergence of order q [CCRZ26]. Notation. Throughout, |·| is the Euclidean norm on Rd , |·|1 is the ℓ1 norm, and (a)+ := max{a, 0}. For probability measures ν ≪ π, we write χ2 (ν∥π) :=

Z µ

¶2 dν − 1 dπ. dπ

e suppresses factors polylogarithmic in d , κ, and ε−1 . A gradient query evaluates The notation O ∇U (y), while a coordinate-partial query evaluates ∂i U (y); d coordinate-partial queries count as one full-gradient equivalent. The letter C denotes a universal constant whose value may change from line to line. Use of AI tools. The general algorithm design and analysis approach are created by the authors. LLM-based assistants were used in the preparation of this manuscript for drafting, language editing, and consistency checks. The authors are responsible for independently checking every statement and proof and for the correctness and integrity of the final manuscript. Upon finishing the draft, we have been informed by Fan Chen, Sinho Chewi and Yiping Lu that ChatGPT Pro 5.6 was able to produce the result for BPS (before our draft becomes public) and in fact establish a high probability guarantee. We decided not to include the stronger result in the manuscript. 2. P RELIMINARIES 2.1. Bouncy particle sampler. The bouncy particle sampler (BPS) is a piecewise deterministic Markov process on Rd × Rd with deterministic flow Ẋ t = Vt , V̇t = 0 between events. For ∇U (x) ̸= 0, define the reflection R x v := v − 2

(2)

〈v, ∇U (x)〉 ∇U (x), |∇U (x)|2

and set R x v = v if ∇U (x) = 0. At state (x, v), bounces occur with rate λBPS (x, v) := (〈v, ∇U (x)〉)+ , and each bounce replaces v by R x v. Independently, refreshment events occur at constant rate γBPS > 0, and each refresh replaces v by an independent draw from ϕ. Accordingly, the generator acts on sufficiently smooth f as (3) ¡ ¢ L BPS f (x, v) = v · ∇x f (x, v) + λBPS (x, v) f (x, R x v) − f (x, v) + γBPS

µZ Rd

¶ f (x, w)ϕ( dw) − f (x, v) .

The BPS process has only finitely many events on bounded time intervals and has invariant law ρ ∞ = µ ⊗ ϕ [BCVD18]. Note that R x is orthogonal, so it preserves |v| and the Gaussian law ϕ. It also reverses the sign of 〈v, ∇U (x)〉; consequently, the speed |Vt | changes only at refreshment times.

6

JIANFENG LU AND YINCHEN LUO

2.2. Zigzag sampler. The Zigzag sampler is a piecewise deterministic Markov process on Rd × Rd with the same deterministic flow Ẋ t = Vt , V̇t = 0 between events. Let F i v := v − 2v i e i denote the velocity obtained by flipping coordinate i . At state (x, v), the i th velocity coordinate flips with rate λZZ (x, v), and the total coordinate-flip rate is ΛZZ (x, v), where i λZZ i (x, v) := (v i ∂i U (x))+ ,

(4)

ΛZZ (x, v) :=

d X i =1

λZZ i (x, v).

At a coordinate-i event, v is replaced by F i v. Independently, refreshment events occur at constant rate γZZ > 0 and replace the entire velocity by an independent draw from ϕ. The corresponding generator is (5) L ZZ f (x, v) = v · ∇x f (x, v) +

d X i =1

¡ ¢ ZZ λZZ i (x, v) f (x, F i v) − f (x, v) + γ

µZ Rd

¶ f (x, w)ϕ( dw) − f (x, v) .

The Zigzag process likewise has only finitely many events on bounded time intervals and has invariant law ρ ∞ = µ ⊗ ϕ [BFR19, LW22a]. Each F i preserves |v j | for every j , and hence preserves both |v| and the Gaussian law ϕ; it also reverses the sign of v i ∂i U (x). Thus, as for BPS, the speed changes only at refreshment times. 2.3. Quantitative χ2 contraction estimates. Our analysis uses the following quantitative contraction estimates of Lu and Wang, stated with the refreshment rates explicit. Theorem 1 (Lu–Wang [LW22b], BPS). Assume (1) and set γBPS =

p d m. There exists a universal

constant K BPS ≥ 1 such that, for every initial law ρ 0 ≪ ρ ∞ with χ2 (ρ 0 ∥ρ ∞ ) < ∞, BPS initialized from ρ 0 satisfies ¶ µ ° ¡ ¢ 1 q t χ2 (ρ 0 ∥ρ ∞ ), t ≥ 0. χ2 Law(X t ,Vt ) ° ρ ∞ ≤ K BPS exp − BPS m d K p Theorem 2 (Lu–Wang [LW22a], Zigzag). Assume (1) and set γZZ = L. There exists a universal

(6)

constant K ZZ ≥ 1 such that, for every initial law ρ 0 ≪ ρ ∞ with χ2 (ρ 0 ∥ρ ∞ ) < ∞, Zigzag initialized from ρ 0 satisfies (7)

µ ° ¡ ¢ χ2 Law(X t ,Vt ) ° ρ ∞ ≤ K ZZ exp −

¶ m p t χ2 (ρ 0 ∥ρ ∞ ), K ZZ L

t ≥ 0.

Throughout, we use the explicit Gaussian cold start (8)

ρ 0 = N (x ⋆ , L −1 I d ) ⊗ N (0, I d ).

The Hessian bounds yield the density comparison (9)

dρ 0 ≤ κd /2 , dρ ∞

χ2 (ρ 0 ∥ρ ∞ ) ≤ κd /2 − 1.

Thus Theorems 1 and 2 apply. The exponentially large initial χ2 divergence enters the mixing horizons through its logarithm, and contributes a factor of d in the complexity bound.

WINDOWED THINNING FOR BPS AND ZIGZAG

7

3. W INDOWED THINNING AND COMPLEXITY THEOREMS The convergence estimates of Section 2 specify how long each process must run. This section combines those continuous-time guarantees with exact simulation algorithms to obtain end-toend complexity bounds. We first introduce windowed Poisson thinning and present its BPS and Zigzag implementations, then state the resulting end-to-end complexity theorems. The detailed query estimates are deferred to Section 4. Along each continuous-time process, events arrive at a time-varying Poisson rate. Thinning generates proposals from an upper rate, or envelope, determined by the trajectory observed so far. A proposal is accepted with probability equal to the ratio of the true rate to this envelope. If the envelope bounds the true rate at every time, the accepted events have exactly the desired law. Windowed thinning keeps this envelope local. We divide the simulation horizon into deterministic windows and evaluate the gradient at the particle’s current position at the beginning of each window. The L-Lipschitz continuity of ∇U then controls how far the true event rate can move from this anchored estimate as the particle travels. The window length sets a natural tradeoff: short windows use more anchor evaluations but keep the envelope tight, whereas long windows cause the envelope to loosen and generate more rejected proposals. We choose the window length by balancing these two effects. 3.1. Windowed thinning for BPS. Recall that for BPS the bounce rate is (〈Vs− , ∇U (X s− )〉)+ . Here and below, s− denotes the left limit, since an event may occur at time s. At the beginning of the window [t k , t k+1 ), we know the anchor gradient G k = ∇U (X tk ), but not its future values along the Rs path. Let D s := tk |Vr | dr denote the cumulative distance travelled from the anchor. Smoothness controls the deviation: |∇U (X s− ) − G k | ≤ L|X s− − X tk | ≤ LD s , which motivates the envelope (10)

λ̄s := (〈Vs− ,G k 〉)+ + L|Vs− |D s .

Indeed, it bounds the true rate at every time: λBPS (X s− ,Vs− ) = (〈Vs− , ∇U (X s− )〉)+ ≤ (〈Vs− ,G k 〉)+ + |〈Vs− , ∇U (X s− ) − G k 〉| ≤ (〈Vs− ,G k 〉)+ + L|Vs− |D s = λ̄s . Between state changes and window boundaries, the velocity is constant and D s is affine. Hence λ̄s is affine on each such interval, and its time integral is quadratic, so the next proposal time can be found by solving a quadratic equation. Algorithm 1 gives the resulting procedure, incorporating this envelope together with the refreshment clock and the next window boundary. In both algorithms, Exp(a) denotes the exponential distribution with rate a. The envelope bound above and standard Poisson thinning show that Algorithm 1 simulates BPS exactly. A proposal can occur only when λ̄t > 0, so its acceptance ratio is well defined. The complete proof that the construction is well defined and exact is given in Proposition 9.

8

JIANFENG LU AND YINCHEN LUO

Algorithm 1 Windowed thinning for BPS Require: horizon T > 0, window length τ > 0, refreshment rate γBPS , initial state (X 0 ,V0 ) 1: t ← 0; partition [0, T ] into windows [t k , t k+1 ), t k = kτ, the last window truncated at T 2: for each nonempty window [t k , t k+1 ) do 3:

query the anchor gradient G k ← ∇U (X tk ); set D tk ← 0

4:

while t < t k+1 do

5:

R tp draw E ∼ Exp(1); generate the next proposal time t p > t by solving t λ̄s ds = E , using (10); set t p = ∞ if there is no solution before t k+1

6:

draw t r ← t + Exp(γBPS ) independently

7:

t ← min{t p , t r , t k+1 }; advance X and D deterministically to time t

8:

if t = t p then

9:

query ∇U (X t − ); draw A ∼ Unif(0, 1)

10:

if A ≤ λBPS (X t − ,Vt − )/λ̄t then Vt ← R X t − Vt − , using the queried gradient

11: 12:

else Vt ← Vt −

13: 14: 15: 16:

end if else if t = t r then draw Vt ∼ ϕ, retaining G k and D

17:

end if

18:

end while

19: end for 20: return (X T ,VT )

3.2. Windowed thinning for Zigzag. Recall that for Zigzag, the flip rate of coordinate i is given by ¡ ¢ λZZ i (X s− ,V s− ) = Vi ,s− ∂i U (X s− ) + . At the beginning of the window [t k , t k+1 ), we know the anchor gradient G k = ∇U (X tk ), but not the Rs future coordinate derivatives along the path. As in the BPS case, let D s := tk |Vr | dr denote the cumulative distance travelled from the anchor. Then, for every coordinate i , |∂i U (X s− ) − G k,i | ≤ |∇U (X s− ) − G k | ≤ L|X s− − X tk | ≤ LD s . This motivates the coordinatewise envelopes and their sum (11)

¡ ¢ λ̄i ,s := Vi ,s−G k,i + + L|Vi ,s− |D s ,

Λ̄s :=

d X

λ̄i ,s .

i =1

Indeed, each envelope bounds its coordinate’s true rate: ¡ ¢ λZZ i (X s− ,V s− ) = Vi ,s− ∂i U (X s− ) + ¡ ¢ ≤ Vi ,s−G k,i + + |Vi ,s− | |∂i U (X s− ) − G k,i | ¡ ¢ ≤ Vi ,s−G k,i + + L|Vi ,s− |D s = λ̄i ,s .

WINDOWED THINNING FOR BPS AND ZIGZAG

9

Summing over the coordinates gives ΛZZ (X s− ,Vs− ) ≤ Λ̄s . The coordinate proposal clocks can therefore be combined into a single clock with rate Λ̄s . When it rings, coordinate i is selected with probability λ̄i ,s /Λ̄s , and only ∂i U (X s− ) needs to be queried to decide whether that flip is accepted. Between state changes and window boundaries, the velocity is constant and D s is affine. Hence the total envelope Λ̄s is affine on each such interval, and its time integral is quadratic, so the next proposal time can be found by solving a quadratic equation. Algorithm 2 gives the resulting procedure, incorporating this envelope together with the refreshment clock and the next window boundary. Algorithm 2 Windowed thinning for Zigzag Require: horizon T > 0, window length τ > 0, refreshment rate γZZ , initial state (X 0 ,V0 ) 1: t ← 0; partition [0, T ] into windows [t k , t k+1 ), t k = kτ, the last window truncated at T 2: for each nonempty window [t k , t k+1 ) do 3:

query the anchor gradient G k ← ∇U (X tk ); set D tk ← 0

4:

while t < t k+1 do

5:

draw E ∼ Exp(1); generate the next proposal time t p > t by solving

R tp t

Λ̄s ds = E , using

(11); set t p = ∞ if there is no solution before t k+1 6:

draw t r ← t + Exp(γZZ ) independently

7:

t ← min{t p , t r , t k+1 }; advance X and D deterministically to time t

8:

if t = t p then

9:

draw coordinate i with probability λ̄i ,t /Λ̄t

10:

query ∂i U (X t − ); draw A ∼ Unif(0, 1)

11:

if A ≤ λZZ (X t − ,Vt − )/λ̄i ,t then i Vt ← Fi Vt −

12: 13:

else Vt ← Vt −

14: 15: 16: 17:

end if else if t = t r then draw Vt ∼ ϕ, retaining G k and D

18:

end if

19:

end while

20: end for 21: return (X T ,VT )

The envelope bounds above and standard Poisson thinning show that Algorithm 2 simulates Zigzag exactly. A coordinate can be proposed only when λ̄i ,t > 0, so its acceptance ratio is well defined. The complete proof that the construction is well defined and exact is given in Proposition 11. 3.3. Main complexity theorems. For Algorithm 1, let Q ∇ (T ) denote the expected number of full-gradient queries up to time T . For Algorithm 2, let Q ∂ (T ) denote the expected number

10

JIANFENG LU AND YINCHEN LUO

of coordinate-partial queries and set Q eq (T ) := Q ∂ (T )/d for the corresponding number of fullgradient equivalents. These quantities measure oracle calls only; a direct implementation of Zigzag also uses O(d ) arithmetic operations per proposal to form and sample from the envelope weights. Combining the mixing estimates of Section 2 with the fixed-horizon query estimates proved in Section 4 gives end-to-end guarantees from the Gaussian cold start. We state the BPS and Zigzag results in parallel, giving the latter in both coordinate-partial and full-gradient-equivalent units. Theorem 3 (Cold-start BPS complexity). Suppose that (1) holds and initialize BPS according to (8). p Algorithm 1 is exact. Set γBPS = d m. For ε ∈ (0, 1/2), define     q s   BPS χ2 (ρ ∥ρ )  BPS  K 0 ∞ d γ  −1/2  , (12) TbεBPS := max 2K BPS log1 + . ,L  m 2ε 4L    At time TbεBPS , the joint χ2 error is at most 4ε2 . Consequently, the joint law and its position marginal have total-variation error at most ε. Choose τ = (Ld )−1/2 . The expected query count satisfies µ µ ¶¶ p 1 (13) Q ∇ (TbεBPS ) ≤ 5 Ld TbεBPS = O κ1/2 d d log κ + log . ε Theorem 4 (Cold-start Zigzag complexity). Suppose that (1) holds and initialize Zigzag according p to (8). Algorithm 2 is exact. Set γZZ = L. For ε ∈ (0, 1/2), define ( ) p µ ¶ K ZZ L K ZZ χ2 (ρ 0 ∥ρ ∞ ) ZZ −1/2 b (14) Tε := max log 1 + ,L . m 4ε2 At time TbεZZ , the joint χ2 error is at most 4ε2 . Consequently, the joint law and its position marginal have total-variation error at most ε. Choose τ = L −1/2 d −1/4 . The expected query counts satisfy µ µ ¶¶ p 1 Q ∂ (TbεZZ ) ≤ 5 L d 5/4 TbεZZ = O κd 5/4 d log κ + log , ε (15) µ µ ¶¶ p 1 Q eq (TbεZZ ) ≤ 5 L d 1/4 TbεZZ = O κd 1/4 d log κ + log . ε Remark 5 (Comparison with Lu–Wang). Lu and Wang [LW22a, Theorem 1] analyze Zigzag query complexity using the global envelope ¡ ¢ (v i ∂i U (x + t v))+ ≤ L|v i | |x| + t |v| , in the normalization x ⋆ = 0. Under their Assumption 2 and stated dimension-and-accuracy regime, an O(1)-χ2 warm start yields the high-probability bound O(κ2 d 3/2 (log 1ε )3/2 ) coordinatepartial queries, or O(κ2 d 1/2 (log(1/ε))3/2 ) full-gradient equivalents. Our locally anchored envelope, 5/4 e together with direct control of the expected event and proposal counts, gives O(κd ) full-gradient

equivalents from the explicit Gaussian cold start (8). We note that [LW22a] provides a highprobability warm-start guarantee under some additional assumptions, while our result only provides guarantee of expected complexity.

WINDOWED THINNING FOR BPS AND ZIGZAG

11

4. P ROOFS The convergence estimates above determine how long the process must be run, but not the cost of producing its trajectory. For the windowed thinning schemes, this cost comes from the gradient evaluations at the window anchors and at the proposed event times. The proposals themselves consist of accepted events, whose frequency is governed by the true PDMP intensities, and rejections, whose frequency is governed by the gap between those intensities and their windowed envelopes. The finite-horizon control needed for this analysis is common to both samplers. Bounces and coordinate flips preserve the speed, while refreshments redraw it, so the refreshment history provides a pathwise upper bound for the position and velocity. Besides controlling the envelope error, this bound supplies the integrability needed for the Dynkin identities used to estimate the accepted-event counts and rules out explosion of the proposal processes. We then treat BPS and Zigzag separately, according to the geometry of their jump mechanisms. 4.1. Common estimates. On [0, T ], for both processes, the speed can change only when the velocity is refreshed. Let NTref denote the number of refreshments and Z j the Gaussian velocity drawn at the j th refreshment. Define ( VTmax := max

)

|V0 |, max |Z j | , 1≤ j ≤NTref

with the inner maximum taken to be zero when NTref = 0. Lemma 6 (Trajectory and velocity-moment bounds). Let (X t ,Vt )t ≥0 be either BPS or Zigzag, and suppose that |X 0 | + |V0 | < ∞ almost surely. Then, for every T < ∞, (16)

sup |Vt | ≤ VTmax ,

0≤t ≤T

sup |X t − x ⋆ | ≤ |X 0 − x ⋆ | + T VTmax .

0≤t ≤T

If V0 ∼ N (0, I d ), then, for every T < ∞ and r > 0, E[(VTmax )r ] < ∞,

(17) and, for every t ≥ 0,

E|Vt |2 = d .

(18)

In the Zigzag case, one also has, for every t ≥ 0, 2 d (d − 1) ≤ d 2 . π If the process is initialized from the cold start (8), then, for every t ≥ 0, p ¡ ¢1/2 p (20) E|X t − x ⋆ |2 ≤ d /L + d t .

(19)

E|Vt |21 = d +

Proof. Bounces and coordinate flips preserve the velocity norm, so only refreshments can change it. This proves the first bound in (16); integrating Ẋ t = Vt proves the second. For (17), NTref is P Poisson and (VTmax )r ≤ |V0 |r + j ≤N ref |Z j |r . The refreshment draws are independent of the clock, T

so the expectation of the sum is ENTref E|Z1 |r < ∞, proving the claim.

12

JIANFENG LU AND YINCHEN LUO

For either process, every non-refreshment update preserves |V |, whereas each refreshment redraws V ∼ N (0, I d ). Since V0 ∼ N (0, I d ), |Vt | has the χd law at every deterministic time, which gives (18). For Zigzag, coordinate flips leave (|V1,t |, . . . , |Vd ,t |) unchanged. At each refreshment this vector is redrawn as the coordinatewise absolute value of an N (0, I d ) vector. Since V0 ∼ N (0, I d ), it has this distribution at every deterministic time t ≥ 0. Hence, with Z ∼ N (0, I d ), E|Vt |21 =

X i

EZi2 + 2

X

E|Zi |E|Z j | = d +

i<j

2 d (d − 1), π

which proves (19). Under the cold start,

¢1/2 p ¡ E|X 0 − x ⋆ |2 = d /L. Rt Moreover, since X t − x ⋆ = X 0 − x ⋆ + 0 Vs ds, we have à ¯Z ¯2 !1/2 ¯ t ¯ ¡ ¡ ¢ ¢ 2 1/2 2 1/2 ¯ E|X t − x ⋆ | ≤ E|X 0 − x ⋆ | + E ¯ Vs ds ¯¯ 0

Z t p ¡ ¢1/2 ds E|Vs |2 ≤ d /L + 0 p p = d /L + d t , where the last equality uses (18). This proves (20) and completes the proof.

□

Under the cold start, Lemma 6 bounds the position and velocity on every finite horizon by a random variable with moments of all orders. This domination also justifies Dynkin’s formula for the unbounded test functions used in the event-count arguments. We isolate the required localization in the following technical lemma, which will be invoked in both sampler-specific proofs below. Lemma 7 (Dynkin formula for observables of polynomial growth). Assume (1) and the cold start (8), and consider either BPS or Zigzag, with generator L . Let f ∈ C 1 (R2d ), and assume that both f and L f have polynomial growth; that is, for some C < ∞ and p ≥ 0, ¡ ¢p | f (x, v)| + |L f (x, v)| ≤ C 1 + |x − x ⋆ | + |v| . Then, for every T < ∞, E f (X T ,VT ) − E f (X 0 ,V0 ) =

Z T 0

E[L f (X t ,Vt )] dt .

Proof. Define stopping times σn := inf{t ≥ 0 : |X t − x ⋆ | + |Vt | ≥ n}. The jump rates are bounded on compact sets, and the Gaussian refreshment kernel preserves integrability of polynomial-growth observables, so the stopped process satisfies Dynkin’s formula: Z T ∧σn E f (X T ∧σn ,VT ∧σn ) − E f (X 0 ,V0 ) = E L f (X t ,Vt ) dt . 0

WINDOWED THINNING FOR BPS AND ZIGZAG

13

For a fixed horizon T , Lemma 6 gives ¡ ¢ sup |X t − x ⋆ | + |Vt | ≤ |X 0 − x ⋆ | + (T + 1)VTmax . 0≤t ≤T

The right-hand side has moments of every order, and σn > T eventually almost surely. The assumed polynomial growth bound therefore supplies a common integrable majorant for both sides of the stopped identity. Letting n → ∞ and applying dominated convergence proves the

□

claim.

4.2. BPS. By Cauchy–Schwarz, controlling the expected number of bounces reduces to controlling the time integral of the squared bounce rate. To estimate this integral using Dynkin’s formula, we seek an observable whose generator contains the squared bounce rate up to a spatial factor. For an observable ψ, the bounce contribution in Dynkin’s formula is ¡ ¢ λBPS (x, v) ψ(x, R x v) − ψ(x, v) , so the jump of ψ should produce the factor 〈v, ∇U (x)〉, since multiplication by λBPS (x, v) then gives λBPS (x, v)2 . For this purpose, we choose ψ(x, v) := 〈x − x ⋆ , v〉. When ∇U (x) ̸= 0, the reflection formula gives (21)

ψ(x, R x v) − ψ(x, v) = −2〈v, ∇U (x)〉

〈x − x ⋆ , ∇U (x)〉 . |∇U (x)|2

We center at x ⋆ so that the smoothness and convexity assumptions bound the spatial ratio on the right-hand side uniformly from below. This observation is the key input to the following estimate. Lemma 8 (Bound on the BPS bounce count from a cold start). Assume (1). For the cold start (8) and every T ≥ 0, Z T (22) 0

£ ¤ γBPS d E λBPS (X t ,Vt )2 dt ≤ Ld T + + 4

p Ld . 2

Hence, whenever T ≥ max{L −1/2 , γBPS /(4L)}, the expected number of bounces in [0, T ] satisfies Z T p £ ¤ (23) B BPS (T ) = E λBPS (X t ,Vt ) dt ≤ 2 Ld T. 0

Proof. Let us introduce the shorthand  〈x − x ⋆ , ∇U (x)〉   , ∇U (x) ̸= 0, |∇U (x)|2 r (x) :=  1/L, ∇U (x) = 0. With this notation, (21) can be written as ψ(x, R x v) − ψ(x, v) = −2〈v, ∇U (x)〉r (x). Strong convexity implies that ∇U (x) = 0 only at x ⋆ , where the bounce rate vanishes; the assigned value at x ⋆ is therefore not important.

14

JIANFENG LU AND YINCHEN LUO

We claim that (24)

L BPS ψ(x, v) = |v|2 − 2λBPS (x, v)2 r (x) − γBPS ψ(x, v).

Indeed, since ∇x ψ(x, v) = v, the transport contribution in (3) is v · ∇x ψ(x, v) = |v|2 . The bounce term follows from the calculation above. Finally, since ϕ is centered, Z Z ψ(x, w)ϕ( dw) = 〈x − x ⋆ , wϕ( dw)〉 = 0, Rd

Rd

so the refreshment contribution is −γBPS ψ(x, v). Combining the three contributions gives (24). The identity (24) involves the squared bounce rate weighted by r . To recover the unweighted rate, we need a uniform lower bound on this factor. We claim that 1 1 ≤ r (x) ≤ , L m

x ∈ Rd .

The claim holds at x ⋆ by the definition of r . Suppose now that x ̸= x ⋆ . Since ∇U (x ⋆ ) = 0, the fundamental theorem of calculus gives Z 1 ∇U (x) = A x (x − x ⋆ ),

A x :=

0

¡ ¢ ∇2U x ⋆ + s(x − x ⋆ ) ds.

A x is symmetric and satisfies mI d ⪯ A x ⪯ LI d by the bound on Hessian of U . Its spectrum is therefore contained in [m, L], and hence m A x ⪯ A 2x ⪯ L A x . Thus we arrived at the desired bound for r combined with the observation that 〈x − x ⋆ , ∇U (x)〉 〈x − x ⋆ , A x (x − x ⋆ )〉 r (x) = = . |∇U (x)|2 〈x − x ⋆ , A 2x (x − x ⋆ )〉 We may now apply Dynkin’s formula. The bounds λBPS (x, v) ≤ L|x − x ⋆ ||v| and r (x) ≤ 1/m show from (24) that |ψ| + |L BPS ψ| has polynomial growth. Therefore Lemma 7 applies and yields Z T £ ¤ Eψ(X T ,VT ) − Eψ(X 0 ,V0 ) = E L BPS ψ(X t ,Vt ) dt 0

Z T =

0

E|Vt |2 dt − 2

− γBPS

Z T 0

Z T 0

£ ¤ E λBPS (X t ,Vt )2 r (X t ) dt

Eψ(X t ,Vt ) dt ,

where the second equality uses (24). By (18), the first term on the right is d T . Rearranging the preceding identity therefore gives Z T Z T £ ¤ 2 E λBPS (X t ,Vt )2 r (X t ) dt = d T − γBPS Eψ(X t ,Vt ) dt 0 0 (25) − Eψ(X T ,VT ) + Eψ(X 0 ,V0 ). We next eliminate the time integral of ψ. The position is continuous, and Ẋ t = Vt between events. Therefore, for almost every t , d |X t − x ⋆ |2 = 2〈X t − x ⋆ ,Vt 〉 = 2ψ(X t ,Vt ). dt

WINDOWED THINNING FOR BPS AND ZIGZAG

15

Integration over [0, T ] and taking expectation therefore gives Z T ¢ 1¡ Eψ(X t ,Vt ) dt = E|X T − x ⋆ |2 − E|X 0 − x ⋆ |2 . 2 0 Consequently, the time-integral term in the rearranged Dynkin identity is Z T ¢ γBPS ¡ −γBPS Eψ(X t ,Vt ) dt = − E|X T − x ⋆ |2 − E|X 0 − x ⋆ |2 . 2 0 We now estimate this term together with the two endpoint terms on the right hand side of identity (25). Under the cold start, X 0 − x ⋆ and V0 are independent and centered, and hence E|X 0 − x ⋆ |2 =

Eψ(X 0 ,V0 ) = E〈X 0 − x ⋆ ,V0 〉 = 0,

d . L

Since E|X T − x ⋆ |2 ≥ 0, it follows that −γBPS

Z T 0

Eψ(X t ,Vt ) dt ≤

γBPS d . 2L

For the remaining endpoint term, Cauchy–Schwarz, (18), and (20) give −E〈X T − x ⋆ ,VT 〉 ≤ |E〈X T − x ⋆ ,VT 〉| ≤

q

s E|X T − x ⋆ |2 E|VT |2 ≤ 

 p d p d + d T  d = p + d T. L L

Inserting these three bounds into (25) yields Z T £ ¤ γBPS d d E λBPS (X t ,Vt )2 r (X t ) dt ≤ 2d T + +p , 2 2L 0 L and the lower bound r ≥ 1/L then gives p Z T £ BPS ¤ γBPS d Ld 2 E λ (X t ,Vt ) dt ≤ Ld T + + , 4 2 0 which proves (22). It remains to deduce the bounce-count estimate. Cauchy–Schwarz gives µ Z T ¶1/2 Z T BPS BPS BPS 2 B (T ) = Eλ (X t ,Vt ) dt ≤ T E[λ (X t ,Vt ) ] dt . 0

If T ≥ L

−1/2

BPS

and T ≥ γ

0

/(4L), then p Ld 1 ≤ Ld T, 2 2

γBPS d ≤ Ld T. 4

Thus (22) implies Z T

5 E[λBPS (X t ,Vt )2 ] dt ≤ Ld T. 2 0 Substitution into the Cauchy–Schwarz bound gives r p 5p Ld T ≤ 2 Ld T, B BPS (T ) ≤ 2 which proves (23).

□

16

JIANFENG LU AND YINCHEN LUO

The bounce count controls the accepted proposals. It remains to count the window anchors and the rejected thinning proposals. Proposition 9 (BPS windowed thinning: exactness and query count). Assume (1) and let T, τ > 0. Run Algorithm 1 with refreshment rate γBPS > 0 and a possibly random initial state (X 0 ,V0 ) satisfying |X 0 | + |V0 | < ∞ almost surely. Then Algorithm 1 is well defined and exactly simulates the BPS process initialized at (X 0 ,V0 ) with refreshment rate γBPS . Denote the process generated by the algorithm by (X t ,Vt ). It satisfies (26)

Q ∇ (T ) ≤

T + 1 + B BPS (T ) + LτT sup E|Vt |2 . τ 0≤t ≤T

For the cold start (8), choosing τ = (Ld )−1/2 gives p Q ∇ (T ) ≤ 5 Ld T

(27) whenever T ≥ max{L −1/2 , γBPS /(4L)}.

Proof. The calculation preceding (10) shows that the envelope dominates the BPS rate throughout each window. Let us show that the proposal process is nonexplosive. Since bounces preserve the speed and only refreshments can change it, for t ≤ T , |Vt − | ≤ VTmax ,

D t ≤ T VTmax ,

|X tk − x ⋆ | ≤ |X 0 − x ⋆ | + T VTmax .

Using |∇U (x)| ≤ L|x − x ⋆ |, we therefore obtain ¡ ¢ |G k | = |∇U (X tk )| ≤ L |X 0 − x ⋆ | + T VTmax . Substituting these bounds into (10) yields λ̄t ≤ |Vt − ||G k | + L|Vt − |D t ¡ ¢ ≤ LVTmax |X 0 − x ⋆ | + 2T VTmax < ∞. Thus the proposal intensity is bounded on [0, T ], which proves nonexplosion, and thus the simulation is well-defined and exact. We next count its gradient queries. There are at most T /τ + 1 nonempty windows and hence at most that many anchor queries. Since every proposal requires one further gradient query, we have Z T T Q ∇ (T ) ≤ + 1 + Eλ̄t dt τ 0 Z T £ ¤ T = + 1 + B BPS (T ) + E λ̄t − λBPS (X t − ,Vt − ) dt . τ 0 The last integral is the expected number of rejected proposals. To control this term, we compare the anchored part of the envelope with the true bounce rate. Since (·)+ is 1-Lipschitz, 0 ≤ λ̄t − λBPS (X t − ,Vt − ) = (〈Vt − ,G k 〉)+ − (〈Vt − , ∇U (X t − )〉)+ + L|Vt − |D t ≤ |〈Vt − ,G k − ∇U (X t − )〉| + L|Vt − |D t ≤ 2L|Vt − |D t .

WINDOWED THINNING FOR BPS AND ZIGZAG

17

The last inequality uses |G k − ∇U (X t − )| ≤ LD t following the smoothness of U . For t ∈ [t k , t k+1 ), Rt write s = t − t k . Since D t = tk |Vu | du, E[|Vt − |D t ] = ≤

Z t

E[|Vt − ||Vu |] du

tk Z tq tk

E|Vt − |2 E|Vu |2 du

≤ s sup E|Vu |2 . 0≤u≤T

If the kth nonempty window has length ℓk ≤ τ, its contribution to the rejected-proposal count is therefore at most Z tk +ℓk tk

£ ¤ E λ̄t − λBPS (X t − ,Vt − ) dt ≤ 2L

µ

sup E|Vu |2

¶Z ℓ

0≤u≤T

k

s ds 0

= Lℓ2k sup E|Vu |2 . 0≤u≤T

Summing over the windows and using Z T 0

P 2 k ℓk ≤ τ k ℓk = τT gives

P

£ ¤ E λ̄t − λBPS (X t − ,Vt − ) dt ≤ LτT sup E|Vu |2 . 0≤u≤T

Substitution into the query decomposition proves (26).

p For the cold start, (18) gives sup0≤t ≤T E|Vt |2 = d , while Lemma 8 gives B BPS (T ) ≤ 2 Ld T on

the stated horizon. With τ = (Ld )−1/2 , (26) becomes p p Q ∇ (T ) ≤ Ld T + 1 + 2 Ld T + L(Ld )−1/2 T d p p = 4 Ld T + 1 ≤ 5 Ld T.

□

This proves (27).

Proof of Theorem 3. The cold start satisfies the initial-state assumption in Proposition 9, and (12) includes the lower bounds required for (27). The proposition therefore shows that Algorithm 1 is well defined on [0, TbεBPS ], simulates BPS exactly, and satisfies (27). The choice of TbεBPS and Theorem 1 give ° ³ ´ ° χ2 Law(X TbBPS ,VTbBPS ) ° ρ ∞ ε ε r µ ¶ 1 m BPS 2 ≤ K BPS exp − BPS Tbε χ (ρ 0 ∥ρ ∞ ) ≤ 4ε2 , d K which gives the total-variation error ε by Cauchy–Schwarz inequality. By (9), 

 q ! Ã p BPS χ2 (ρ ∥ρ ) K 0 ∞  K BPS κd /4  log1 +  ≤ log 1 + 2ε 2ε µ ¶ 1 ≤ C d log κ + log . ε

18

JIANFENG LU AND YINCHEN LUO

p The other two entries in (12) are also bounded by d /m, since L −1/2 ≤ m −1/2 and γBPS /(4L) ≤ p κ−1 d /m. Therefore, absorbing some constants, s  µ ¶ 1 d TbεBPS = O  d log κ + log  . m ε Since

p p p Ld d /m = d κ, the query estimate becomes µ µ ¶¶ p 1 BPS BPS 1/2 b b Q ∇ (Tε ) ≤ 5 Ld Tε = O κ d d log κ + log , ε

□

which proves (13).

4.3. Zigzag. We follow a similar strategy as the analysis of BPS, while in the case of Zigzag, we would need to control the quadratic aggregated flip rate: S(x, v) :=

d X i =1

2 λZZ i (x, v) .

To estimate its time integral using Dynkin’s formula, we choose the observable q(x, v) := v · ∇U (x). A coordinate flip changes the corresponding summand of q so that the jump part of its generator produces S. The remaining terms are controlled by the Hessian bound and the moment estimates from Lemma 6. Lemma 10 (Bound on the Zigzag flip count from a cold start). Assume (1). Consider Zigzag initialized from the cold start (8), with arbitrary refreshment rate γZZ > 0. Then, for every T ≥ 0, Z T (28) 0

γZZ d ES(X t ,Vt ) dt ≤ Ld T + + 4

p Ld , 2

and hence, whenever T ≥ max{L −1/2 , γZZ /(4L)}, the expected total number of flips satisfies (29)

B ZZ (T ) =

Z T 0

p EΛZZ (X t ,Vt ) dt ≤ 2d L T.

Proof. We first compute the generator of q: (30)

L ZZ q(x, v) = v ⊤ ∇2U (x)v − 2S(x, v) − γZZ q(x, v).

Indeed, the transport contribution is v ⊤ ∇2U (x)v. For a flip of coordinate i , q(x, F i v) − q(x, v) = −2v i ∂i U (x),

ZZ 2 λZZ i (x, v)v i ∂i U (x) = λi (x, v) ,

so summing the flip terms gives −2S(x, v). The Gaussian velocity law is centered, and hence the refreshment term is −γZZ q(x, v). This proves (30).

WINDOWED THINNING FOR BPS AND ZIGZAG

19

Since ∇U (x ⋆ ) = 0, the L-smoothness and mean-value theorem give |∇U (x)| ≤ L|x − x ⋆ |. Consequently, |q(x, v)| ≤ |v| |∇U (x)| ≤ L|x − x ⋆ ||v|, S(x, v) ≤

d X

¡ ¢2 v i2 ∂i U (x) ≤ |v|2 |∇U (x)|2 ≤ L 2 |x − x ⋆ |2 |v|2 ,

i =1 2

|v ∇ U (x)v| ≤ ∥∇ U (x)∥op |v|2 ≤ L|v|2 . ⊤

2

These estimates show that q and L ZZ q have polynomial growth. We may thus apply Lemma 7. Using (30), Dynkin’s formula gives Eq(X T ,VT ) − Eq(X 0 ,V0 ) =

Z T 0

£ ¤ E Vt⊤ ∇2U (X t )Vt dt Z T

−2

0

ES(X t ,Vt ) dt − γZZ

Z T 0

Eq(X t ,Vt ) dt .

To express the last integral in terms of the potential, note that X t is continuous and Ẋ t = Vt for almost every t . The chain rule along each trajectory gives Z T Z T ∇U (X t ) · Vt dt = q(X t ,Vt ) dt . U (X T ) −U (X 0 ) = 0

0

Substituting this identity into Dynkin’s formula and rearranging yields Z T Z T £ ¤ 2 ES(X t ,Vt ) dt = E Vt⊤ ∇2U (X t )Vt dt 0

0

£ ¤ − γZZ E U (X T ) −U (X 0 ) − Eq(X T ,VT ) + Eq(X 0 ,V0 ). Since V0 is centered and independent of X 0 , Eq(X 0 ,V0 ) = 0. The Hessian bound and (18) imply Z T Z T £ ⊤ 2 ¤ E Vt ∇ U (X t )Vt dt ≤ L E|Vt |2 dt = Ld T. 0

0

Since U (X T ) ≥ U (x ⋆ ), the upper Hessian bound and the cold start give −E[U (X T ) −U (X 0 )] ≤ E[U (X 0 ) −U (x ⋆ )] ≤

L d E|X 0 − x ⋆ |2 = . 2 2

For the terminal term, Lemma 6 yields ¯ ¯ −Eq(X T ,VT ) ≤ ¯Eq(X T ,VT )¯ q ≤ L E|X T − x ⋆ |2 E|VT |2 s  p p p d ≤L + d T  d = L d + Ld T. L Inserting these bounds into the Dynkin identity gives Z T γZZ d p 2 ES(X t ,Vt ) dt ≤ 2Ld T + + L d, 2 0 which proves (28).

20

JIANFENG LU AND YINCHEN LUO

Finally, (ΛZZ )2 ≤ d S, so Cauchy–Schwarz on [0, T ] × Ω gives µ Z T ¶1/2 ZZ ZZ 2 B (T ) ≤ T E[Λ (X t ,Vt ) ] dt 0

µ ¶1/2 Z T ≤ dT ES(X t ,Vt ) dt . 0

Under the stated lower bound on T , the last two terms in (28) are at most Ld T and Ld T /2, respectively. Hence Z T 0

5 ES(X t ,Vt ) dt ≤ Ld T, 2

B

ZZ

r (T ) ≤

p 5 p d L T ≤ 2d L T, 2

□

which proves (29).

The flip count controls the accepted proposals. It remains to count the window anchors and the rejected coordinate proposals. Proposition 11 (Zigzag windowed thinning: exactness and query count). Assume (1) and let T, τ > 0. Run Algorithm 2 with refreshment rate γZZ > 0 and a possibly random initial state (X 0 ,V0 ) satisfying |X 0 | + |V0 | < ∞ almost surely. Then Algorithm 2 is well defined and exactly simulates the Zigzag process initialized at (X 0 ,V0 ) with refreshment rate γZZ . Denote the process generated by the algorithm by (X t ,Vt ). It satisfies ¶ µ ¶1/2 µ ¶1/2 µ T . + 1 + B ZZ (T ) + LτT sup E|Vt |21 sup E|Vt |2 (31) Q ∂ (T ) ≤ d τ 0≤t ≤T 0≤t ≤T If V0 ∼ N (0, I d ), then (32)

Q eq (T ) ≤

For the cold start (8), setting γZZ = (33)

p T B ZZ (T ) +1+ + L d τT. τ d

p L and choosing τ = L −1/2 d −1/4 gives p Q eq (T ) ≤ 5 L d 1/4 T

whenever T ≥ L −1/2 . Proof. The calculation preceding (11) shows that every coordinate envelope dominates its flip rate throughout each window. Let us show that the proposal process is nonexplosive. Since coordinate flips preserve the speed and only refreshments can change it, for t ≤ T , |Vt − | ≤ VTmax , D t ≤ T VTmax , |X tk − x ⋆ | ≤ |X 0 − x ⋆ | + T VTmax . p Moreover, |Vt − |1 ≤ d VTmax . Using |∇U (x)| ≤ L|x − x ⋆ |, we therefore obtain ¡ ¢ |G k | = |∇U (X tk )| ≤ L |X 0 − x ⋆ | + T VTmax . Substituting these bounds into (11) yields Λ̄t ≤ |Vt − ||G k | + L|Vt − |1 D t ¡ ¢ p ≤ LVTmax |X 0 − x ⋆ | + T VTmax + d LT (VTmax )2 < ∞.

WINDOWED THINNING FOR BPS AND ZIGZAG

21

Thus the total proposal intensity is bounded on [0, T ], which proves nonexplosion, and thus the simulation is well-defined and exact. We next count its coordinate-partial queries. There are at most T /τ + 1 nonempty windows and hence at most d (T /τ + 1) coordinate-partial anchor queries. Since every proposal requires one further coordinate-partial query, we have µ ¶ Z T T Q ∂ (T ) ≤ d +1 + EΛ̄t dt τ 0 ¶ µ Z T £ ¤ T + 1 + B ZZ (T ) + E Λ̄t − ΛZZ (X t − ,Vt − ) dt . =d τ 0 The last integral is the expected number of rejected coordinate proposals. To control this term, we compare the anchored part of each coordinate envelope with its true flip rate. Since (·)+ is 1-Lipschitz, 0 ≤ Λ̄t − ΛZZ (X t − ,Vt − ) =

d ¡¡ X ¢ ¡ ¢ ¢ Vi ,t −G k,i + − Vi ,t − ∂i U (X t − ) + + L|Vi ,t − |D t

i =1

≤

d X

|Vi ,t − | |G k,i − ∂i U (X t − )| + L|Vt − |1 D t

i =1

≤ 2L|Vt − |1 D t . The last inequality uses |G k,i − ∂i U (X t − )| ≤ |G k − ∇U (X t − )| ≤ LD t . For t ∈ [t k , t k+1 ), write s = t − t k . Rt Since D t = tk |Vu | du, Z t E[|Vt − |1 D t ] = E[|Vt − |1 |Vu |] du tk

≤

Z tq tk

µ ≤s

E|Vt − |21 E|Vu |2 du

sup E|Vu |21

¶1/2 µ

0≤u≤T

sup E|Vu |2

¶1/2 .

0≤u≤T

If the kth nonempty window has length ℓk ≤ τ, its contribution to the rejected-proposal count is therefore at most µ ¶1/2 µ ¶1/2 Z ℓ Z tk +ℓk k £ ¤ E Λ̄t − ΛZZ (X t − ,Vt − ) dt ≤ 2L sup E|Vu |21 sup E|Vu |2 s ds tk

0≤u≤T

= Lℓ2k

µ

sup E|Vu |21

0≤u≤T

0

0≤u≤T

¶1/2 µ

sup E|Vu |2

¶1/2 .

0≤u≤T

P P Summing over the windows and using k ℓ2k ≤ τ k ℓk = τT bounds the rejection term by µ ¶1/2 µ ¶1/2 2 2 LτT sup E|Vu |1 sup E|Vu | . 0≤u≤T

0≤u≤T

Substitution into the query decomposition proves (31). If V0 ∼ N (0, I d ), then (18) and (19) give µ ¶1/2 µ ¶1/2 p sup E|Vt |21 sup E|Vt |2 ≤ d d. 0≤t ≤T

0≤t ≤T

22

JIANFENG LU AND YINCHEN LUO

Combined with (31) now gives (32) as Q eq (T ) = d1 Q ∂ (T ). p For the cold start, γZZ = L and T ≥ L −1/2 imply T ≥ γZZ /(4L). Lemma 10 therefore gives p B ZZ (T )/d ≤ 2 L T . With τ = L −1/2 d −1/4 , (32) becomes p p p Q eq (T ) ≤ L d 1/4 T + 1 + 2 L T + L d 1/4 T p p ≤ 4 L d 1/4 T + 1 ≤ 5 L d 1/4 T. p The last inequality uses L d 1/4 T ≥ 1. This proves (33). □ Proof of Theorem 4. The cold start satisfies the initial-state assumption in Proposition 11, and (14) includes the lower bound required for (33). The proposition therefore shows that Algorithm 2 is well defined on [0, TbεZZ ], simulates Zigzag exactly, and satisfies (33). The choice of TbεZZ and Theorem 2 yield ° ³ ´ ° χ2 Law(X TbZZ ,VTbZZ ) ° ρ ∞ ε ε ¶ µ m ≤ K ZZ exp − p TbεZZ χ2 (ρ 0 ∥ρ ∞ ) ≤ 4ε2 . K ZZ L This gives the total-variation error ε by Cauchy–Schwarz inequality. By (9), Ã ! ¶ µ K ZZ κd /2 K ZZ χ2 (ρ 0 ∥ρ ∞ ) log 1 + ≤ log 1 + 4ε2 4ε2 µ ¶ 1 ≤ C d log κ + log . ε p p −1/2 The other entry in (14) is also bounded by L/m, since L ≤ L/m. Therefore, absorbing some constants, TbεZZ = O Since

Ãp µ ¶! L 1 d log κ + log . m ε

p p L( L/m) = κ, the query estimates become µ µ ¶¶ p 1 Q ∂ (TbεZZ ) ≤ 5 L d 5/4 TbεZZ = O κd 5/4 d log κ + log , ε ¶¶ µ µ p 1 Q eq (TbεZZ ) ≤ 5 L d 1/4 TbεZZ = O κd 1/4 d log κ + log . ε

□

This proves (15). R EFERENCES

[ABKR25] Sanket Agrawal, Joris Bierkens, Kengo Kamatani, and Gareth O. Roberts. Transient regime of piecewise deterministic Monte Carlo algorithms, 2025. arXiv:2509.16062. [AC24]

Jason M. Altschuler and Sinho Chewi. Faster high-accuracy log-concave sampling via algorithmic warm starts. Journal of the ACM, 71(3):24:1–24:55, 2024.

[ACZ26]

Jason M. Altschuler, Sinho Chewi, and Matthew S. Zhang. Shifted composition IV: Toward ballistic acceleration for log-concave sampling. In Proceedings of the 58th Annual ACM Symposium on Theory of Computing, pages 1739–1750. ACM, 2026.

[ADNR21] Christophe Andrieu, Alain Durmus, Nikolas Nüsken, and Julien Roussel. Hypocoercivity of piecewise deterministic Markov process–Monte Carlo. The Annals of Applied Probability, 31(5):2478–2517, 2021.

WINDOWED THINNING FOR BPS AND ZIGZAG

[AK24]

23

Charly Andral and Kengo Kamatani. Automated techniques for efficient sampling of piecewise-deterministic Markov processes, 2024. arXiv:2408.03682.

[BCVD18] Alexandre Bouchard-Côté, Sebastian J. Vollmer, and Arnaud Doucet. The bouncy particle sampler: A nonreversible rejection-free Markov chain Monte Carlo method. Journal of the American Statistical Association, 113(522):855–867, 2018. [BFR19]

Joris Bierkens, Paul Fearnhead, and Gareth Roberts. The Zig-Zag process and super-efficient sampling for Bayesian analysis of big data. The Annals of Statistics, 47(3):1288–1320, 2019.

[BGKR20] Joris Bierkens, Sebastiano Grazzi, Kengo Kamatani, and Gareth Roberts. The boomerang sampler. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 908–918, 2020. [BKR22]

Joris Bierkens, Kengo Kamatani, and Gareth O. Roberts. High-dimensional scaling limits of piecewise deterministic sampling algorithms. The Annals of Applied Probability, 32(5):3361–3407, 2022.

[BKR25]

Joris Bierkens, Kengo Kamatani, and Gareth O. Roberts. Scaling of piecewise deterministic Monte Carlo for anisotropic targets. Bernoulli, 31(3):2323–2350, 2025.

[BR17]

Joris Bierkens and Gareth O. Roberts. A piecewise deterministic scaling limit of lifted Metropolis–Hastings in the Curie–Weiss model. The Annals of Applied Probability, 27(2):846–882, 2017.

[BRB25]

Jeremiah Birrell and Luc Rey-Bellet. Concentration inequalities and UQ bounds for hypocoercive MCMC samplers. Quarterly of Applied Mathematics, 2025. Electronically published November 14, 2025.

[BRZ19]

Joris Bierkens, Gareth O. Roberts, and Pierre-André Zitt. Ergodicity of the zigzag process. The Annals of Applied Probability, 29(4):2266–2301, 2019.

[CCDR26] Fan Chen, Sinho Chewi, Constantinos Daskalakis, and Alexander Rakhlin. High-accuracy sampling for diffusion models and log-concave distributions, 2026. arXiv:2602.01338. [CCRZ26] Fan Chen, Sinho Chewi, Alexander Rakhlin, and Matthew S. Zhang. Private communication. 2026. [Che26]

Sinho Chewi. Log-concave sampling. https://chewisinho.github.io/main.pdf, 2026. Book draft, accessed on Jul 30, 2026.

[CLW23]

Yu Cao, Jianfeng Lu, and Lihan Wang. On explicit L 2 -convergence rate estimate for underdamped Langevin dynamics. Archive for Rational Mechanics and Analysis, 247(5):90, 2023.

[CSR22]

Alice Corbella, Simon E. F. Spencer, and Gareth O. Roberts. Automatic Zig-Zag sampling in practice. Statistics and Computing, 32:107, 2022.

[DBCD19] George Deligiannidis, Alexandre Bouchard-Côté, and Arnaud Doucet. Exponential ergodicity of the bouncy particle sampler. The Annals of Statistics, 47(3):1268–1287, 2019. [DCWY19] Raaz Dwivedi, Yuansi Chen, Martin J. Wainwright, and Bin Yu. Log-concave sampling: Metropolis–Hastings algorithms are fast. Journal of Machine Learning Research, 20(183):1–42, 2019. [DGM20]

Alain Durmus, Arnaud Guillin, and Pierre Monmarché. Geometric ergodicity of the bouncy particle sampler. The Annals of Applied Probability, 30(5):2069–2098, 2020.

[EL26]

Andreas Eberle and Francis Lörler. Non-reversible lifts of reversible diffusion processes and relaxation times. Probability Theory and Related Fields, 194(1–2):173–203, 2026.

[FBPR18]

Paul Fearnhead, Joris Bierkens, Murray Pollock, and Gareth O. Roberts. Piecewise deterministic Markov processes for continuous-time Monte Carlo. Statistical Science, 33(3):386–412, 2018.

[LL26]

Bowen Li and Jianfeng Lu. Space-time log-Sobolev inequality and hypocoercive hypercontractivity for underdamped Langevin dynamics, 2026. arXiv:2605.25083.

[LS79]

P. A. W. Lewis and G. S. Shedler. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly, 26(3):403–413, 1979.

[Lu26]

Jianfeng Lu. A sharp hypocoercive entropy decay estimate for underdamped Langevin dynamics, 2026. arXiv:2605.01933.

[LW22a]

Jianfeng Lu and Lihan Wang. Complexity of Zigzag sampling algorithm for strongly log-concave distributions. Statistics and Computing, 32:48, 2022.

24

JIANFENG LU AND YINCHEN LUO

[LW22b]

Jianfeng Lu and Lihan Wang. On explicit L 2 -convergence rate estimate for piecewise deterministic Markov processes in MCMC algorithms. The Annals of Applied Probability, 32(2):1333–1361, 2022.

[MDS20]

Manon Michel, Alain Durmus, and Stéphane Sénécal. Forward event-chain Monte Carlo: Fast sampling by randomness control in irreversible Markov chains. Journal of Computational and Graphical Statistics, 29(4):689–702, 2020.

[MW26]

Pierre Monmarché and Lihan Wang. On the entropic convergence for piecewise deterministic samplers: Speedup and obstruction, 2026. arXiv:2606.26086.

[Oga81]

Yosihiko Ogata. On Lewis’ simulation method for point processes. IEEE Transactions on Information Theory, 27(1):23–31, 1981.

[PdW12]

E. A. J. F. Peters and G. de With. Rejection-free Monte Carlo sampling for general potentials. Physical Review E, 85(2):026703, 2012.

[SF23]

Matthew Sutton and Paul Fearnhead. Concave-convex PDMP-based sampling. Journal of Computational and Graphical Statistics, 32(4):1425–1435, 2023.

[WR20]

Changye Wu and Christian P. Robert. Coordinate sampler: A non-reversible Gibbs-like MCMC sampler. Statistics and Computing, 30(3):721–730, 2020.

[WSC22]

Keru Wu, Scott Schmidler, and Yuansi Chen. Minimax mixing time of the Metropolis-Adjusted Langevin Algorithm for log-concave sampling. Journal of Machine Learning Research, 23(270):1–63, 2022.

D EPARTMENT OF M ATHEMATICS , D EPARTMENT OF P HYSICS , AND D EPARTMENT OF C HEMISTRY, D UKE U NIVERSITY Email address: [email protected] D EPARTMENT OF M ATHEMATICS , D UKE U NIVERSITY Email address: [email protected]