High-Dimensional Behaviour of Two Piecewise Deterministic Markov Processes

Efficiency Comparison & Asymptotic Variance Estimation

Hirofumi Shiba

Institute of Statistical Mathematics, Tokyo, Japan

10/13/2026

Scaling Analysis of Two PDMPs

• Section 1: \;What is PDMP?

  A new class of continuous-time Monte Carlo dynamics

• Section 2: \;What is scaling analysis?

  FECMC always performs better than BPS

• Section 3: \;Proof strategy

• Section 4: \;Asymptotic variance estimation (≒ ESS estimation)

  The PDMP framework allow a more efficient estimator

1 PDMP: A New Frontier of Monte Carlo

Markov Chain
(1953–)

Langevin Diffusion
(1978–)

PDMP
(2008–)

A Blog Entry on Bayesian Computation by an Applied Mathematician

$$

$$

1.1 Piecewise-Deterministic Markov Process

Langevin Diffusion

Randomized Hamiltonian Monte Carlo

1.2 Markov Chain Monte Carlo

computes space averages as time averages

1.3 Persistent Exploration with Hypocoercivity

Underdamped Langevin

Kinetic Langevin

\large \begin{cases} \dot{X}_t=V_t,\\ \dot{V}_t=-\big(\nabla U(X_t)+\gamma V_t\big)\,dt+\sqrt{2\gamma}\,dB_t. \end{cases}

Forward Event-Chain Monte Carlo

Velocity-Jump PDMP

\large \begin{cases} \dot{X}_t=V_{t-},\\ \dot{V}_t=\big(\operatorname{Jump}_{X_t}(V_{t-};\xi_t)-V_{t-}\big)\,dN_t. \end{cases}

2 A Scaling Analysis

We compare two famous PDMP methods under the following condition:

• ODE: Linear flow
• Jump: Reflection + Refreshment vs. Combination
• Target: High dimensional standard Gaussian U(x):=-\log\pi(x)=\|x\|^2/2.

2.1 BPS vs. FECMC: Which is Most Efficient?

Bouncy Particle Sampler
(Bouchard-Côté et al., 2018)

BPS with different hyperparameter \rho

Forward Event Chain Monte Carlo (Michel et al., 2020)

2.2 Jump Mechanisms in BPS vs. FECMC

Reflection

\Large +

Refreshment \rho

\Large\fallingdotseq

Stochastic Reflection

2.3 Empirical Comparison: BPS vs. FECMC

Effective Sample Size of U(x)=\|x\|^2/2, the negative log-density of 100 dimensional standard Gaussian, estimated over 1000 runs

2.4 Scaling Analysis of MCMCs

aims to identify the scaling factor \textcolor{#E95420}{d^{\textrm{\,some \,factor}}} that yields a diffusion limit as d\to\infty

\text{Plotting } \textcolor{#0096FF}{Y_t^{(d)}}=\frac{\|\textcolor{#0096FF}{X}_{\textcolor{#E95420}{d}\textcolor{#0096FF}{t}}^{\textcolor{#0096FF}{(d)}}\|^2-d}{\sqrt{d}} \text{ with } d=10^2,10^3,10^4:

2.5 Theorem 1: Diffusion Limits of the Rescaled Energy

dY_t^{\textcolor{#0096FF}{\text{B}}}=-\frac{\sigma^2_{\textcolor{#0096FF}{\text{B}}}(\rho)}{4}Y_t^{\textcolor{#0096FF}{\text{B}}}\,dt+\sigma_{\textcolor{#0096FF}{\text{B}}}(\rho)\,dB_t \sigma^2_{\textcolor{#0096FF}{\text{B}}}(\rho)=8\int^\infty_0e^{-\rho s}\operatorname{E}[R_0^{\textcolor{#0096FF}{\text{B}}}R_s^{\textcolor{#0096FF}{\text{B}}}]\,ds

dY_t^{\textcolor{#E95420}{\text{F}}}=-\frac{\sigma^2_{\textcolor{#E95420}{\text{F}}}(\rho)}{4}Y_t^{\textcolor{#E95420}{\text{F}}}\,dt+\sigma_{\textcolor{#E95420}{\text{F}}}(\rho)\,dB_t \sigma^2_{\textcolor{#E95420}{\text{F}}}(\rho)=8\int^\infty_0e^{-\rho s}\operatorname{E}[R_0^{\textcolor{#E95420}{\text{F}}}R_s^{\textcolor{#E95420}{\text{F}}}]\,ds

2.6 Theorem 2: Analytic Expressions of \sigma^2

\sigma^2_{\textcolor{#E95420}{\text{FECMC}}}(0)=\sqrt{\frac{32}{\pi}},\quad\sigma^2_{\textcolor{#0096FF}{\text{BPS}}}(0)=0, \begin{align*} \sigma^2_{\textcolor{#E95420}{\text{FECMC}}}(\rho)&=\sqrt{\frac{32}{\pi}}\biggr(1-\frac{\left(\rho^2-\rho\sqrt{\frac{\pi}{2}}+\Omega(\rho)\right)^2}{\rho^4\Omega(\rho)(2-\Omega(\rho))}\biggl),\quad\rho>0\\ \sigma^2_{\textcolor{#0096FF}{\text{BPS}}}(\rho)&=\frac{8}{\rho^4}\left(\rho^3-\rho^2\sqrt{\frac{8}{\pi}}+\rho-\sqrt{\frac{8}{\pi}}\frac{\left((1+\rho^2)\Omega(\rho)-\rho^2\right)^2}{\Omega(2\rho)}\right), \end{align*}

where \Omega(\rho)\coloneqq\sqrt{\frac{\pi}{2}}\rho\operatorname{erfcx}\left(\frac{\rho}{\sqrt{2}}\right)=\rho e^{\frac{\rho^2}{2}}\int^\infty_\rho e^{-\frac{t^2}{2}}\,dt satisfies \Omega(\rho)=1-\operatorname{E}[e^{-\rho\tau}] for \tau\sim\operatorname{Rayleigh}(1).

2.7 Theorem 2 (cont.): Plot of \sigma^2

While BPS atteins maximum at non-trivial value of \rho, FECMC achieves maximum at \rho=0

3 Proofs

• The hypocoercivity / slow-fast structure complicates the standard proof approach
  • The potential U(X_t)=\|X_t\|^2/2 evolves on a slow timescale.
  • Its time derivative \textcolor{#E95420}{R_t}:=(X_t|V_t)=dU(X_t)/dt fluctuates on a fast timescale.
• Spectral properties of \textcolor{#E95420}{R} are reflected in \sigma^2

3.1 The Non-Markovianity of the Rescaled Energy Y^d

Y^d_{t-} alone cannot determine the speed of Y^d_t at jumps.

3.2 Goal: Solve the Limiting Martingale Problem

For any \varphi\in C^\infty_c(\mathbb{R}), the following is a local martingale: M^\varphi_t:=\varphi(Y_t)-\varphi(Y_0)-\int^t_0\mathcal{L}\varphi(Y_s)\,ds.

Lifted Generator
(Ethier and Kurtz, 1986)

1. Find a Markovian tuple Z^d:=(Y^d,R^d,\cdots) with generator \mathcal{L}^d.
2. Guess the limiting generator \mathcal{L}\varphi=b\varphi'+\frac{c}{2}\varphi''.
3. Prove that the difference vanishes: \scriptsize \int^t_0|\mathcal{L}^d(\varphi\circ\mathrm{pr}_1)(Z^d_s)-\mathcal{L}\varphi(Y^d_s)|\,ds\xrightarrow{\operatorname{P}}0.

Semimartingale Characteristics (Jacod and Shiryaev, 2003)

1. Obtain a canonical semimartingale decomposition Y^d_t=Y^d_0+B^d_t+M^d_t.
2. Prove the convergence of semimartingale characteristics triplet: B^d,\quad C^d:=\langle M^d\rangle,\quad\nu^d.
3. Recover the generator via Ito’s formula \tiny \int^t_0(\mathcal{L}\varphi)(Y_s)\,ds =\varphi'(Y_t)\,dB_t +\frac12\varphi''(Y_t)\,dC_t

Same foundation, different language.

3.3 The Event-Time Skeleton Reveals Local Averaging

Original Potential Process \textcolor{#0096FF}{Y^d}

B^d_t=2\sqrt{d}\int^t_0R^d_{ds-}\,ds,%\quad R_t^d:=\brac{X_t^d,V_t^d}, C^d_t=0.

Difficult to see how C_t=\sigma^2t emerges.

approx.

Event-Time Skeleton \textcolor{#0096FF}{\overline{Y}^d}

\overline{B}^d_t=\sum_{n=1}^{N_t}\operatorname{E}[\Delta\overline{Y}^d_n|\overline{\mathcal{F}}^d_{n-1}] {\tiny \overline{C}^d_t=\sum_{n=1}^{N_t}\left(\operatorname{E}[(\Delta\overline{Y}^d_n)^2|\overline{\mathcal{F}}^d_{n-1}]-\operatorname{E}[\Delta\overline{Y}^d_n|\overline{\mathcal{F}}^d_{n-1}]^2\right).}

Leads directly to the limit C_t=\sigma^2t.

We then prove the asymptotic equivalence of these two.

3.4 \sigma^2 as the CLT Variance of Time-Derivative \textcolor{#E95420}{R}

From Theorem 1, we see \begin{align*} \frac{\sigma^2(\rho)}{4}=2\int^\infty_0e^{-\rho t}\operatorname{E}[\textcolor{#E95420}{R_0R_t}]\,dt=\lim_{T\to\infty}\frac{1}{T}\operatorname{Var}\left[\int^T_0\textcolor{#E95420}{R_t^{(\rho)}}\,dt\right]. \end{align*} Therefore, \sigma^2/4 is the asymptotic variance of the ergodic average of \textcolor{#E95420}{R_t^{(\rho)}}.

3.5 \sigma^2(\rho) and Resolvent Equations \textcolor{#E95420}{R}

From Theorem 1, we see \begin{align*} \frac{\sigma^2(\rho)}{4}=2\int^\infty_0e^{-\rho t}\operatorname{E}[\textcolor{#E95420}{R_0R_t}]\,dt=\lim_{T\to\infty}\frac{1}{T}\operatorname{Var}\left[\int^T_0\textcolor{#E95420}{R_t^{(\rho)}}\,dt\right]. \end{align*} Therefore, \sigma^2/4 is the asymptotic variance of the ergodic average of \textcolor{#E95420}{R_t^{(\rho)}}.

Pushing the calculation further using Fubini’s theorem, =2\operatorname{E}\biggl[\textcolor{#E95420}{R_0}\underbrace{\int^\infty_0e^{-\rho t}\operatorname{E}[\textcolor{#E95420}{R_t}|\textcolor{#E95420}{R_0}]\,dt}_{=:f_\rho(\textcolor{#E95420}{R_0})}\biggr]=2\operatorname{E}[\textcolor{#E95420}{R_0}f_\rho(\textcolor{#E95420}{R_0})], where f_\rho solves the resolvent equation of \textcolor{#E95420}{R} (\rho-\mathcal{L}_{R})f_\rho(x)=x where \mathcal{L}_{R} is the generator of \textcolor{#E95420}{R}. This can be explicitly solved. (cf. Bierkens and Lunel, 2022)

3.6 The Most Simple Case: \sigma^2(0) and Poisson Equations

Taking \rho\to0, the resolvent equation reduces to -\mathcal{L}_Rf_0(x)=x\quad\text{Poisson eq.} This can be solved in each sampler:

Sampler	Solution	\sigma^2(0)=8\mathbb{E}[R_0f_0(R_0)]
BPS	\displaystyle f_0^{\textcolor{#2780e3}{B}}(x)=\frac{1-x^2}{2}	0
FECMC	\displaystyle f_0^{\textcolor{#E95420}{F}}(x)=\frac14-\frac12\min(x,0)^2	\displaystyle\sqrt{\frac{32}{\pi}}

The limiting mixing speed of Y^d, governed by \sigma^2, is determined by the accumulated fluctuations of the fast process R^d.

Kipnis-Varadhan theory (Cattiaux et al., 2011) can systematically guarantee a CLT and that its variance is \sigma^2(0).

4 Asymptotic Variance Estimation for PDMP in High Dimensions

\frac{16}{\sigma^2}=\lim_{T\to\infty}T\operatorname{Var}\!\left[\frac1T\int_0^T \textcolor{#2780e3}{Y_t}\,dt\right], \sigma^2=4\lim_{T\to\infty}T\operatorname{Var}\!\left[\frac1T\int_0^T \textcolor{#E95420}{R_t}\,dt\right].

Estimators for quadratic variation \sigma^2 are equivalent to batch means estimators applied to time derivative R, providing a fast alternative to those for \textcolor{#0096FF}{Y_t^{(d)}}=\frac{\|\textcolor{#0096FF}{X}_{\textcolor{#E95420}{d}\textcolor{#0096FF}{t}}^{\textcolor{#0096FF}{(d)}}\|^2-d}{\sqrt{d}}.

4.1 ESS of \textcolor{#0096FF}{Y_t}

Informal Corollary (Relationship between \sigma^2 and ESS)

The effective sample size (ESS) for estimating \textcolor{#0096FF}{U} with a trajectory of length \textcolor{#E95420}{d}T is given by \operatorname{ESS}(\textcolor{#0096FF}{U})\fallingdotseq\frac{\sigma^2}{8}T\qquad(d\to\infty).

The Monte Carlo estimator (= ergodic average of \textcolor{#0096FF}{\{X_t\}}) \widehat{h}_T^d:=\frac{1}{T}\int^T_0\textcolor{#0096FF}{U(X_t^{(d)})}\,dt has the following variance ratio under the double limit: \lim_{T\to\infty}\lim_{d\to\infty}T\frac{\operatorname{Var}_\pi[\widehat{h}_{\textcolor{#E95420}{d}T}^d]}{\operatorname{Var}_\pi[\textcolor{#0096FF}{U}]}=\frac{8}{\sigma^2}\approx2.50\cdots. ESS equals the inverse of this variance ratio.

4.2 Theory Matches Practice: FECMC vs. BPS

ESS for U estimated over 1000 runs. Every bootstrap CI contains the theoretical limiting value {8}/{\sigma^2}.

4.3 \sigma^2 is unknown for General Targets

ESS for U estimated over 1000 runs.

4.4 Asymptotic Variance Estimation for Discrete-Time MCMC

A practical issue: Does this trajectory suffice for mean estimation? → A solution: asymptotic variance estimation

\text{Markov Chain CLT: }\quad\frac{1}{\sqrt{N}}\sum_{n=1}^NX_n\Rightarrow N(0,\sigma^2_{\text{asym}})\quad(N\to\infty). For mini-batches b=1,\cdots,B with length m:=N/B, \operatorname*{{Batch\; Means\; Estimator:}} \text{ }\quad\widehat{\sigma^2}_\text{asym}:=\frac{m}{B-1}\sum_{b=1}^B\left(\textcolor{#0096FF}{\overline{{X}}^{(b)}}-\textcolor{#0096FF}{\overline{{X}}}\right)^2, where \textcolor{#0096FF}{\overline{{X}}^{(b)}}=\frac{1}{m}\sum_{n=1}^m\textcolor{#0096FF}{X_{n+(b-1)m}} is the mean over the b-th batch.

4.5 Optimal Batch Size for Discrete-Time MCMC & PDMP

Due to the continuous-time nature, batch size selection becomes a nontrivial problem for PDMPs.

MSE Minimizing Batch Size for Classical MCMC (Liu et al., 2022)

m_{\text{opt}}=\left(\frac{\gamma_{\text{disc}}}{\sigma_{\text{asym}}^2}\right)^{\frac{2}{3}}N^{\frac{1}{3}} minimizes the MSE asymptotically under m\to\infty,\qquad N/m\to\infty, for a sufficiently mixing Markov chain \{\textcolor{#0096FF}{X_n}\}.

Bias Minimizing Batch Size for PDMP (S. & Kamatani 2026+)

Analogously, with \gamma_{\text{disc}}=2\sum_{k=0}^\infty k\operatorname{Cov}[\textcolor{#0096FF}{X_0},\textcolor{#0096FF}{X_k}] replaced by \gamma_{\text{cont}}=2\int^\infty_0t\operatorname{Cov}[\textcolor{#E95420}{X_0},\textcolor{#E95420}{X_t}]\,dt.

4.6 Diffusivity Estimation for High-Dimensional PDMPs

= Quadratic Variation (QV) Estimation under Misspecification

Plot of Y_t^d:=(|\textcolor{#E95420}{X^d}_{d\textcolor{#E95420}{t}}|^2-d)/\sqrt{d} when d=10^4. The trajectory converges to dY_t=-\frac{\sigma^2}{4} Y_t\,dt+\sigma\,dB_t as d\to\infty

The Realized Volatility estimator (Barndorff-Nielsen and Shephard, 2002) \widehat{Q}^d:=\frac{1}{T}\sum_{n=1}^N\biggr(Y_{\Delta n}^d-Y_{\Delta(n-1)}^d\biggl)^2,\quad T=\Delta N, converges to \sigma^2 under the conditions \Delta\to0,\quad d\to\infty,\quad\Delta d\to\infty.

4.7 Experiment: QV Estimator Improves with Dimension

Two estimators for the CLT variance 8/\sigma^2. slow corresponds to the Batch Means estimator, fast corresponds to the Realized Variation estimator

4.8 Experiment (cont.): QV Estimator Improvement

Two estimators for the CLT variance 8/\sigma^2. slow corresponds to the Batch Means estimator, fast corresponds to the Realized Variation estimator

4.9 Batch Size Selection for the QV Estimator

\fallingdotseq selecting the sampling step size \Delta under Misspecification (Aït-Sahalia et al., 2011)

Plot of Y_t^d:=(|\textcolor{#E95420}{X^d}_{d\textcolor{#E95420}{t}}|^2-d)/\sqrt{d} when d=10^4. Taking a closer look, the trajectory is an ODE, a `microstructure noise’

For too small / large \Delta, \widehat{Q}^{d}(\Delta)\approx0. The optimal choice seems to be \Delta=O(d^{1/2}).

Conclusion

1. Reducing redundant noise produces a more efficient algorithm: BPS → FECMC
   • FECMC has larger \sigma^2 because its noise in R^d is effective
   • Radial mixing in piecewise-linear PDMPs is essentially determined by R^d
2. Continuous-time MCMC allows more efficient estimators for asymptotic variance
   • potentially works whenever a diffusion limit exists
   • with a sensible choice of the batch size
3. Future direction: From diagnostics for U to global convergence guarantees.

References

Aït-Sahalia, Y., Mykland, P. A., and Zhang, L. (2011). Ultra High Frequency Volatility Estimation with Dependent Microstructure Noise. Journal of Econometrics, 160(1), 160–175.

Barndorff-Nielsen, O. E., and Shephard, N. (2002). Estimating Quadratic Variation using Realized Variance. Journal of Applied Econometrics, 17(5), 457–477.

Bierkens, J., Fearnhead, P., and Roberts, G. (2019). The Zig-Zag Process and Super-Efficient Sampling for Bayesian Analysis of Big Data. The Annals of Statistics, 47(3), 1288–1320.

Bierkens, J., and Lunel, S. M. V. (2022). Spectral analysis of the zigzag process. Annales de l’Institut Henri Poincaré, Probabilités Et Statistiques, 58(2), 827–860.

Bouchard-Côté, A., Vollmer, S. J., and Doucet, A. (2018). The Bouncy Particle Sampler: A Nonreversible Rejection-Free Markov Chain Monte Carlo Method. Journal of the American Statistical Association, 113(522), 855–867.

Cattiaux, P., Chafai, D., and Guillin, A. (2011). Central limit theorems for additive functionals of ergodic markov diffusions processes.

Deligiannidis, G., Paulin, D., Bouchard-Côté, A., and Doucet, A. (2021). Randomized Hamiltonian Monte Carlo as scaling limit of the bouncy particle sampler and dimension-free convergence rates. The Annals of Applied Probability, 31(6), 2612–2662.

Ethier, S. N., and Kurtz, T. G. (1986). Markov processes: Characterization and convergence. John Wiley & Sons, Inc.

Jacod, J., and Shiryaev, A. N. (2003). Limit Theorems for Stochastic Processes,Vol. 288. Springer Berlin, Heidelberg.

Liu, Y., Vats, D., and Flegal, J. M. (2022). Batch Size Selection for Variance Estimators in MCMC. Methodology and Computing in Applied Probability, 24(1), 65–93.

Michel, M., Durmus, A., and Sénécal, S. (2020). Forward Event-Chain Monte Carlo: Fast Sampling by Randomness Control in Irreversible Markov Chains. Journal of Computational and Graphical Statistics, 29(4), 689–702.

Appendix

PDMP as a New Discretisation Scheme

Many PDMP and MCMC behave similarly in high dimensions (Deligiannidis et al., 2021)

RHMC (d=2)

discretized by a symplectic integrator

e.g. HMC has O(d^{1.25}) complexity

BPS with Gaussian speed (d=10^3)

simulate piecewise linear trajectory

e.g. BPS with normal velocity has O(d^{1.5}) complexity

In Passing… Some Killer Applications of PDMP

x: CPU time, y: Estimate．(Bouchard-Côté et al., 2018)
Sparse Markov Random Field with d=10

O(d^1) Local Implementation

Exploiting sparsity, BPS atteins better scaling than HMC (MSE/sec.)

x: sample size n, y: log ESS/sec, d=16
(Bierkens et al., 2019)

Stochastic Gradient O(n^0)

Using an appropriate control variate, Zig-Zag attains O(1) complexity
in the limit of n\to\infty, it outperforms Langevin.

課題：非一様 Poisson 点過程のシミュレーション

PDMP シミュレーションの原則

方向転換地点 \textcolor{#78C2AD}{x_1,x_2,\cdots} を決めれば良い

＊これが難しい

\textcolor{#78C2AD}{x_1,x_2,\cdots} の決め方

確率密度 \pi に収束させたい場合，強度 m^{(i)}(t)=\biggr(-v\cdot\nabla\log\pi(\textcolor{#78C2AD}{x_{i-1}}+tv)\biggl)_+,\quad t\ge0, を持つ非一様 Poisson 点過程をシミュレートする．

最初の到着点を \textcolor{#78C2AD}{x_i}=\textcolor{#78C2AD}{x_{i-1}}+vT_1^{(i)} とする．

直観：１次元の場合

１次元の Gauss 分布を考える：

\pi(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}} \therefore\quad U(x):=-\log\pi(x)=\frac{x^2}{2\sigma^2}+\mathrm{const.}

負の対数尤度 U=-\log\pi は ポテンシャル ともいう．

\begin{align*} m^{(1)}(t)&=(v_0\cdot U'(\textcolor{#78C2AD}{x_0}+\textcolor{#E95420}{v_0}t))_+\\ &=(\textcolor{#78C2AD}{x_0}+\textcolor{#E95420}{v_0}t)\cdot1{\left\{\textcolor{#78C2AD}{x_0}+\textcolor{#E95420}{v_0}t>0\right\}},\quad t\ge0, \end{align*} を強度にもつ 非一様 Poisson 点過程 の最初の到着時刻 T_1^{(1)} について，\textcolor{#78C2AD}{x_1}=\textcolor{#78C2AD}{x_0}+\textcolor{#E95420}{v_0}T_1^{(1)} とすれば良い．