Generalizing Parallel Replica Dynamics: Trajectory Fragments, Asynchronous Computing, and PDMPs
David Aristoff


Connectivity of hyperplane sections of domains
    
Matteo Varbaro Supported by the program for abroad research activity of the University of Genoa, Pr. N, 100021-2016-MV-ALTROUNIGE_001, and by the GNSAGA-INdAM 

Dip. di Matematica
 Università  di Genova


<varbaro@dima.unige.it> 
    
==========================================================================================================================================================================================================================================




 
We study the Parallel Replica 
Dynamics in a general 
setting. We introduce 
a trajectory fragment framework that can 
be used to design and prove 
consistency of Parallel 
Replica algorithms 
for generic Markov 
processes. We use our 
framework to formulate 
a novel condition that 
guarantees an 
asynchronous algorithm 
is consistent. Exploiting 
this condition and 
our trajectory fragment 
framework, we
present new 
synchronous and 
asynchronous Parallel Replica 
algorithms for piecewise 
deterministic Markov processes.




Parallel Replica Dynamics, long time dynamics, stationary distributions, asynchronous computing, piecewise deterministic Markov processes




  65C05, 65C20, 65C40, 
  65Y05, 82C80




§ INTRODUCTION


Many problems in applied sciences 
require the sampling of 
complex probability distributions. 
In computational chemistry – 
which is the main setting 
of this article – such 
distributions can arise from 
stochastic models of molecular 
dynamics <cit.> or chemical reaction networks <cit.>, 
while obstacles to efficient 
sampling include high dimensionality 
and metastability, the latter being 
the tendency to become stuck in certain 
subsets of state space <cit.>. 
Some attempts to surmount 
these difficulties have been based on importance sampling and 
stratification <cit.>, interacting particles <cit.>, coarse graining and preconditioning <cit.>, accelerated dynamics <cit.> 
and nonreversibility <cit.>.

This article concerns the 
Parallel Replica Dynamics
(ParRep) <cit.>, 
an accelerated 
dynamics method designed to 
overcome metastability. ParRep has 
two distinct advantages over many other 
enhanced sampling methods. First, 
it computes correct dynamical <cit.> as well as 
stationary or equilibrium <cit.> quantities associated with a stochastic process. 
And second, ParRep is very general: it only 
requires mild assumptions on the underlying process. 
Indeed, though originally 
intended for Langevin dynamics <cit.>, straightforward 
extensions of ParRep 
to discrete and 
continuous time 
Markov chains have appeared in <cit.>.


The goal of this article 
is as follows. First, we introduce  a new mathematical framework that may be used to design, and 
prove consistency of, ParRep algorithms for Markov processes satisfying a few mild assumptions. Second, we use our framework 
to obtain valuable insights 
about 
asynchronous computing. 
In particular, we present specific, novel 
conditions that ensure an asynchronous ParRep algorithm 
is consistent.
Lastly, we construct 
ParRep algorithms 
for piecewise deterministic 
Markov processes (PDMPs) in both the 
synchronous and asynchronous setting, leaning on our new 
framework to demonstrate their consistency. 
Asynchronous ParRep 
algorithms must be carefully designed, since as we show below, 
inaccuracies can arise
when the speed of computing paths of the underlying process is 
coupled to the process itself.

PDMPs are emerging as an useful 
tool in fields as 
diverse as applied probability <cit.>, computational chemistry <cit.>, machine 
learning <cit.>, 
and big data <cit.>. 
As indicated by the name, PDMPs
move along deterministic paths in between 
random jump times. In the context of chemical reaction networks,
PDMPs called hybrid models can be obtained 
by approximating
fast reactions by a deterministic 
flow, and representing slow reactions with an 
appropriate Poisson process <cit.>. The 
resulting PDMPs can be metastable <cit.>, 
making direct simulation
unattractive. Several 
PDMP-based algorithms have 
also been proposed for sampling from distributions known up to normalization – like the Boltzmann distribution or the posterior distribution in Bayesian 
analysis – including
Event chain Monte Carlo <cit.>, the Zig-Zag process <cit.>, 
and the Bouncy Particle Sampler <cit.>. 
Below, we give a general argument 
suggesting these PDMPs also 
become metastable under certain 
conditions. 


This article is organized as follows.  Section <ref> 
defines notation that 
we use throughout. In Section <ref>, 
we formally define metastability 
in terms of quasistationary distributions.  
We describe ParRep in more detail, and explain what we mean 
by a consistent ParRep 
algorithm, in Section <ref>. 
In Section <ref> 
we outline a general mathematical 
framework for ParRep, and in 
Section <ref> we study 
synchronous and 
asynchronous computing. 
Section <ref> 
serves as a brief introduction 
to PDMPs, while 
Section <ref> 
outlines several ParRep 
algorithms for PDMPs that 
are based on our framework 
from Section <ref>. 
A numerical example is in 
Section <ref>.
All proofs are in Section <ref>.



§ NOTATION


Throughout,
X(t)_t ≥ 0 is 
a time 
homogeneous Markov process, either discrete or continuous 
in time, with values in a standard Borel state space; 
U is a subset 
of state space; and g 
is a real-valued 
function defined on 
state space. 
Without 
explicit mention we assume
all sets are measurable and all functions are 
bounded and measurable.
We write 
X(t) to refer 
to the process X(t)_t ≥ 0 at 
time t. We 
denote various expectations and probabilities by
𝔼 and ℙ, 
with the precise meaning being 
clear from context. 
We write ℒ  
for the probability law of a 
random object, with ℒ 
above an equals sign indicating
equality in law. We say a random 
object is a copy of another random 
object if it has the same law as that object. 
When we say a collection of random 
objects is independent we mean 
these objects are mutually independent unless otherwise specified.
We define a∧ b = min{a,b} and 
a ∨ b =max{a,b}, 
and write ⌊ s⌋ 
for the greatest integer less than or equal to s.




§ METASTABILITY

Informally, U is a metastable 
set for X(t)_t ≥ 0 if 
X(t)_t ≥ 0 tends to reach 
a local equilibrium in U much faster than it escapes from U. Local equilibrium 
can be understood in terms of quasistationary distributions (QSDs):


Fix a subset
U of state space, and consider

    T = inf{t≥ 0:X(t) ∉ U},

the 
first time X(t)_t ≥ 0 escapes U.
A QSD ρ of X(t)_t ≥ 0 in 
U satisfies ρ(U) = 1 and

    ρ(A) = ℙ(X(t) ∈ A|ℒ(X(0)) = ρ, T>t)

for every t ≥ 0 and A ⊆ U. 


Note that ρ is supported in U. Equation (<ref>) states 
that if X(0) is distributed as ρ and X(t)_t ≥ 0
does not escape from U by time t, then X(t) is distributed 
as ρ.
Throughout, we will assume the QSD of X(t)_t ≥ 0 in 
U exists, is unique, and is the long time distribution of X(t) conditioned to never escape U. 
That is, we assume that for any initial distribution of X(0) supported in U,

    ρ(A) = lim_t →∞ℙ(X(t) ∈ A|X(s) ∈ U  for s ∈ [0,t])   ∀ A ⊆ U.

The QSD ρ can then be sampled as follows: choose a time  T_corr^ρ(U) 
for relaxation to ρ. 
Start X(t)_t ≥ 0 in U, 
and if it escapes from
U before time t = T_corr^ρ(U), restart it in 
U. Repeat this until 
a trajectory of X(t)_t ≥ 0 
remains in U for a consecutive time 
interval of length T_corr^ρ(U). 
This trajectory's terminal position is 
then a sample of ρ. For more 
details on the QSD see for 
instance <cit.>. 
For conditions ensuring existence of and 
convergence to the QSD for general 
Markov processes, see <cit.>.

A more formal definition of metastability is: a set U is metastable for X(t)_t ≥ 0
if the time scale to reach ρ is small 
compared to the mean time to escape from U starting at ρ. In 
some cases these times can be written
in terms of the eigenvalues 
of the adjoint, L^*, of the generator L of X(t)_t ≥ 0,
with absorbing boundary conditions on the complement of U.
See <cit.> 
and <cit.> 
for the corresponding 
spectral analysis for
overdamped Langevin dynamics and finite state space discrete and continuous time Markov chains, and <cit.> 
for an application of these 
ideas to choosing T_corr^ρ(U).




§ PARALLEL REPLICA DYNAMICS


ParRep can boost the efficiency of 
simulating metastable processes <cit.>.
Currently, implementations have 
been proposed only for Langevin 
or overdamped Langevin dynamics <cit.> 
and discrete or continuous time Markov chains <cit.>. However, the 
generality of ParRep allows 
for extensions to any metastable
time homogeneous strong Markov 
process with càdlàg paths, in cases where the QSD exists and
metastable sets can be identified. 
We make this precise in 
the next section.

ParRep algorithms are based 
on two basic steps: 
5pt


  * A step in which X(t)_t ≥ 0 is allowed to reach the QSD in some metastable set U, 
using direct or serial simulation – called the decorrelation step;

  * A step generating an escape event from U, starting from the QSD, using parallel simulation – called the parallel step.

5pt

By escape event 
we mean the random pair (T,X(T)), where 
T is the time for 
X(t)_t ≥ 0 
to escape from U when 
X(0) is distributed as the QSD 
in U, and X(T) is the 
corresponding
escape point. The parallel 
step efficiently computes an escape 
event starting from the QSD 
via a sort of time parallelization.


The decorrelation step, as 
it uses only serial 
simulation, is exact. 
By exact we mean there 
is zero error – except for 
the inevitable error 
in simulating X(t)_t≥ 0 
arising from numerical discretizations, 
which we will ignore.
Our analysis 
will therefore focus on the 
parallel step. 

The parallel step is sometimes divided 
into two subrouties: first, 
a routine that generates 
independent samples 
of the QSD in U – called dephasing – and 
second, a routine that 
uses copies of X(t)_t ≥ 0 
starting from these QSD samples to generate 
an escape event of X(t)_t ≥ 0 
from U. Below, we will mostly omit discussion of the dephasing 
routine, and we will not discuss the error 
associated with imperfect convergence 
to the QSD in the dephasing and decorrelation 
steps, as these points have been previously studied in <cit.>.


We say the parallel step of a ParRep algorithm is 
consistent when:
5pt


  * The parallel step generates  
escape events from each metastable set U with the correct probability law – see Theorem <ref> below; 

  * The parallel step produces 
correct mean contributions to 
time averages in each metastable set U – see Theorem <ref> below.

5pt

By correct we mean exact 
provided the QSD sampling has zero error. 
A consistent ParRep algorithm 
defines a 
coarse dynamics, that is,
a dynamics that is correct 
on the quotient space obtained 
by considering each metastable 
set as a single point <cit.>. A 
consistent ParRep algorithm also 
defines stationary averages that are 
correct for functions 
defined on the original 
uncoarsened state 
space <cit.>. ParRep produces only a coarse dynamics because the parallel step does not resolve the
exact behavior of X(t)_t ≥ 0. The parallel step is faithful enough to X(t)_t≥ 0, however,
to produce correct stationary averages on the original uncoarsened state space <cit.>. 











Previous analyses of ParRep 
have relied on the structure 
of X(t)_t ≥ 0 and the particular algorithms 
studied <cit.>. We introduce a new
framework below that allows us 
to study the consistency of 
any ParRep algorithm.  Our 
analysis is inspired by 
ParSplice, a recent implementation 
of ParRep employing asynchronous 
computing <cit.>. 
Our framework provides 
explicit conditions that 
ensure an asynchronous ParRep algorithm 
is consistent. In particular, it shows a certain class 
of asynchronous ParRep 
algorithms is consistent 
provided the wall-clock time 
to simulate a step of 
X(t)_t ≥ 0
is not coupled to its 
position in state space; see 
Section <ref> 
below for precise statements.




§ TRAJECTORY FRAGMENTS


We now formalize 
conditions which lead 
to consistency of ParRep. 
Our arguments are based on what 
we call trajectory 
fragments. The 
fragments are copies of the underlying process 
satisfying the dependency conditions of 
Assumption <ref> below. 
In practice, the trajectory fragments 
may be computed asynchronously 
in parallel. We discuss this in the 
next section. 







Let X(t)_t ≥ 0 have 
càdlàg paths and 
the strong Markov property. 
Assume X(t)_t ≥ 0 has a QSD ρ
in an open set U and that T = inf{t ≥ 0:X(t) ∉ U} is finite almost surely. 
Let
(X_m(t)_t≥ 0,T_m)_m ≥ 1 be copies of (X(t)_t ≥ 0,T) such that:

    conditional 
    on X_m(0), X_m(t)_t ≥ 0 is 
    independent of (X_k(t)_t ≥ 0)_1 ≤ k < m; 
       ℒ(X_m(0)|T_m-1>t_m-1,…,T_1> t_1) = ρ for  m ≥ 2,   ℒ(X_1(0)) = ρ,

where t_m>0 are deterministic times satisfying ∑_m=1^∞ t_m = ∞.






The X_m(t)_0 ≤ t ≤ t_m 
from Assumption <ref> 
are the trajectory 
fragments. We will 
refer to X_m(t)_t > t_m 
as a fragment's irrelevant future. The 
reason for this choice 
of words is that the 
output of a general 
parallel step, described 
in Algorithm <ref>, is the same no matter how the X_m(t)_t ≥ 0 
are defined for times t>t_m.













Algorithm <ref>  
outlines a general parallel step. As discussed above, this parallel step 
can be combined with a decorrelation step 
to compute a coarse dynamics or 
a time average of a function g.
The idea behind 
Algorithm <ref> 
is simple – we  
imagine concatenating 
fragments 
whose starting points 
and terminal points are distributed 
as the QSD ρ, thus  
obtaining an artificial
long trajectory. See Figure <ref> below.
One must be careful, however, 
in treating
dependencies of the fragments. 
The dependencies described 
in Assumption <ref> lead 
to a consistent Algorithm <ref> 
in the sense of Theorems <ref>-<ref>. 
More general dependencies 
can violate consistency, 
as we will discuss in 
the next section.
 













Our next two results demonstrate 
consistency of Algorithm <ref> under 
the conditions in 
Assumption <ref>. 
Theorem <ref> 
states that Algorithm <ref> produces
the correct escape law 
from U starting at the 
QSD in U, 
while Theorem <ref> says that Algorithm <ref> produces 
the correct mean contribution to time averages.


Suppose that Assumption <ref> holds. Let X(t)_t ≥ 0 be such that ℒ(X(0)) = ρ, and set T = inf{t>0:X(t) ∉ U}. Then in Algorithm <ref>,  
    (T_par,X_par) ℒ= (T,X(T)).










Suppose that Assumption <ref> holds. Let X(t)_t ≥ 0 be such that ℒ(X(0)) = ρ, and set T = inf{t>0:X(t) ∉ U}. Then in Algorithm <ref>, in the discrete time case, 

    g_par := 𝔼(∑_m=1^L∑_t=0^T_m ∧ t_m-1 g(X_m(t))) = 𝔼(∑_t=0^T-1g(X(t))),

while in the continuous time case, 

    g_par := 𝔼(∑_m=1^L∫_0^T_m ∧ t_m g(X_m(t)) dt) = 𝔼(∫_0^Tg(X(t)) dt).



Recall that the gain in ParRep is 
from parallel computations in the parallel 
step. In our trajectory 
fragment framework,
the basic idea is that the work 
to compute the fragments X_m(t)_0 ≤ t ≤ t_m 
can be spread over multiple 
processors.  See <cit.> for related results 
in special cases. 
The parallel step is more 
efficient than direct, or serial, simulation 
provided the computational 
effort to sample the QSD is 
small relative to the effort 
to simulate an escape from U 
via serial simulation.

We actually do not 
need to assume that U is open and
that X(t)_t ≥ 0 has the strong 
Markov property and 
càdlàg paths to prove 
consistency of the parallel 
step in Algorithm <ref>. 
Indeed, Algorithm <ref>
is consistent in the
sense of Theorems <ref>-<ref>  
whenever X(t)_t ≥ 0 is a time 
homogeneous Markov process and (<ref>)-(<ref>) 
hold. 
However, to combine the parallel 
step with a decorrelation step 
to obtain a coarse dynamics or 
stationary average, 
we want càdlàg paths 
to ensure the escape time from an 
open set U is a 
stopping time, and we need 
the strong Markov property 
so that we can start afresh 
at these stopping times. 



 





§ SYNCHRONOUS 
AND ASYNCHRONOUS COMPUTING
 
Recall that the speedup in ParRep 
comes from computing the trajectory fragments 
X_m(t)_0 ≤ t ≤ t_m partly 
or fully in parallel. These 
fragments must be ordered, 
via the index m ≥ 1, to 
obtain the long trajectory 
pictured in Figure <ref>. 
Below, we explore two possible
ways to 
order the fragments, 
depending on whether 
we want to employ
synchronous 
or asynchronous 
computing. In the 
former case, we have 
in mind
a computing environment 
consisting of R 
processors that are 
nearly synchronous. 
In the latter case 
we consider an 
arbitrary number of 
processors 
that potentially 
have widely different 
performance. 

Below, we will consider only 
fragments of constant time 
length, t_m ≡Δ t. 
For synchronous 
computing, following 
ideas
from <cit.>, we 
consider the 
ordering of trajectory 
fragments in Proposition <ref> below. 





Suppose Y^r(t)_t ≥ 0, 
r=1,…,R, are independent 
copies of X(t)_t ≥ 0 
with ℒ(X(0)) = ρ. 
Let m_k = ⌊ (k-1)/R⌋ and r_k = k-R⌊ (k-1)/R⌋, and for k ≥ 1 define 
trajectory fragments X_k(t)_0 ≤ t ≤Δ t by

    X_k(t) = Y^r_k(m_kΔ t + t),      0 ≤ t ≤Δ t.

Then Assumption <ref> holds with an appropriate 
definition of the trajectory
fragments' irrelevant 
futures. Thus the conclusions of Theorem <ref> and <ref> hold.

Figure <ref> shows 
the trajectory fragments defined in (<ref>).



In asynchronous 
computing, perhaps the most natural ordering 
is the wall-clock ordering: 
the fragments are ordered 
according to the wall-clock 
time that their starting points 
are computed. See Figure <ref>.
When does the wall-clock ordering satisfy Assumption <ref>? 
Note that (<ref>) simply says that each 
fragment evolves forward in time independently 
of the preceding fragments and 
their irrelevant futures. 
This condition is easy to establish 
with
an appropriate choice of the 
irrelevant futures.
Ensuring (<ref>) holds is more subtle. We will 
show, however, that 
if the wall-clock time it takes 
to compute each fragment depends on processor 
variables, but not on the fragments themselves, 
then the wall-clock time ordering is independent of 
the fragments and (<ref>) holds.



We will distinguish 
between a wall-clock time and 
a physical time, where the former 
is self-explanatory and the
latter refers to the time index 
t of a copy of X(t)_t ≥ 0. 
Let Y^r(t)_t ≥ 0 be independent 
copies of 
X(t)_t ≥ 0 starting at 
independent samples of the 
QSD ρ. The wall-clock time 
ordering of fragments satisfies 
Assumption <ref> above if (i) the wall-clock times are independent 
of the physical times, 
(ii) the wall-clock time 
to compute copy Y^r(t)_t ≥ 0 
is an increasing function of 
the physical time, and 
(iii) two processors 
never finish at exactly the same 
wall-clock time, so that the wall-clock 
times can be given a 
unique ordering. Write t_wall^r(m) 
for the wall-clock time it takes to 
compute Y^r(t)_t ≥ 0 up to 
physical time t = mΔ t. Proposition <ref> below makes 
the claims above precise:







 
Suppose Y^r(t)_t ≥ 0, 
r=1,…,R, are independent 
copies of X(t)_t ≥ 0 
with ℒ(X(0)) = ρ. 
Assume (t_wall^r(m)_m ≥ 0)^1 ≤ r ≤ R are nonnegative random numbers such that: 
5pt


  (i)
(t_wall^r(m)_m ≥ 0)^1 ≤ r ≤ R is independent of (Y^r(t)_t ≥ 0)^1≤ r ≤ R; 
2pt

  (ii) Almost surely, t_wall^r(m) ≤ t_wall^r(n) when 
m ≤ n and 1 ≤ r ≤ R;
 2pt

  (iii) Almost surely, 
there is a unique sequence
 (r_k,m_k)_k ≥ 1 such that:
 
    (r_k,m_k)_k ≥ 1 has range {1,…,R}×{0,1,2,…}    (surjectivity),
       t_wall^r_1(m_1)<t_wall^r_2(m_2)<t_wall^r_3(m_3)<…        (monotonicity).


5pt
For k ≥ 1 define trajectory fragments X_k(t)_0 ≤ t ≤Δ t by

    X_k(t) = Y^r_k(m_kΔ t + t),      0 ≤ t ≤Δ t.

Then Assumption <ref> holds 
with an appropriate definition 
of the trajectory fragments' irrelevant futures.
Thus the conclusions of 
Theorems <ref> and <ref> hold.


Figure <ref> 
shows the trajectory fragments X_m(t)_0 ≤ t ≤ t_m defined in (<ref>). 

Assumptions (ii) and (iii) are 
quite natural, but assumption
(i) can 
fail in many very ordinary settings. We 
sketch an example explaining 
how this could happen. 
Suppose
U = (0,1) ⊆ℝ 
and say X(t)_t ≥ 0 obeys some one dimensional stochastic differential equation. 
Suppose we use an integrator 
for X(t)_t ≥ 0 that is 
slow near 0 but fast near 
1. Then the wall-clock ordering 
will likely put trajectory fragments 
that are near 1 
ahead of those near 0. This 
bias in the ordering 
would in turn create a bias 
toward escaping through 1: 
in Algorithm <ref>, we 
would expect that ℙ(X_par= 1) > ℙ(X_T = 1), where X_T is the correct escape point. We 
construct a specific 
example demonstrating 
this bias in Remark <ref> 
in Section <ref> below.

The speed of integrators does 
commonly depend on position in 
space, particularly when 
the time step varies to account 
for numerical stiffness <cit.>. This is an important caveat 
to keep in mind for asynchronous 
algorithms. This  
issue has not been explored 
much in the literature; see 
however brief discussions in <cit.> and <cit.>.

The setting of Parsplice <cit.> is slightly 
different from the above. In ParSplice, 
a splicer tells a producer 
to generate fragments among several 
metastable sets. The splicer 
distributes the fragments 
according to where it speculates 
that they will be needed.
These 
fragments are given a 
label
as soon as they are assigned, 
and this label never changes. The 
labels are assigned 
in wall-clock time 
order. Thus label 
i is less than label j if and only if 
the splicer tells 
the producer to generate
fragment i before it 
tells the producer
to generate fragment j.
When the splicer tells the producer 
to generate a fragment 
in a particular metastable 
set U, it takes as 
its starting point the 
terminal point of the fragment 
in U with the smallest label. 
Crucially, this label 
is smallest among all fragments 
in U and not just among 
fragments in U which have been 
fully computed at the current 
wall-clock time. 
Thus, the ordering of fragments 
in U, fixed by the splicer, can 
be seen as independent of 
the fragments themselves, 
and the arguments above demonstrate 
consistency of ParSplice.





§ PDMPS



The remainder of this article will focus on 
applying our ideas above to PDMPs. 
We begin with a brief informal
description of PDMPs. 
A PDMP is a 
càdlàg process consisting 
of a deterministic dynamics 
interrupted by jumps at 
random times; formally, a PDMP in ℝ^d has a generator L of the form

    Lf(z) = ∂_Γ(z) f(z) + λ(z) ∫ (f(z')-f(z))Q(z,dz')

acting on suitable f:ℝ^d →ℝ. Here λ(z) is
the jump rate at z ∈ℝ^d, a Markov kernel
Q(z,dz') describes the jump distribution, and Γ defines 
the deterministic flow 
    ∂_t ψ(t,z) = Γ(ψ(t,z)),    ψ(0,z) = z.



Write θ_0 + … + θ_n-1 for the nth jump time, so that
θ_n-1 is the 
holding time before the nth  
jump, and 
write ξ_n for the position immediately after the 
nth jump, with ξ_0 the initial position. Then 
the PDMP generated by (<ref>) is described by ψ together 
with 
(ξ_n,θ_n)_n ≥ 0; 
we call the latter the skeleton chain of 
Z(t)_t ≥ 0. Note that the 
skeleton chain is a
time homogeneous Markov chain. 


For convenience we describe a 
way to simulate a PDMP described 
by (<ref>) in
Algorithm <ref> below. In 
the algorithm, we 
abuse notation by writing 
θ_n, ξ_n, and 
Z(t) for particular 
realizations of these 
random objects.





Sampling the times θ_n is a nontrivial task, but there 
are efficient methods 
based on Poisson 
thinning <cit.> and 
identifying certain critical points along the flow direction <cit.>. See also <cit.> for other 
methods to simulate a PDMP, including some 
based on time discretization. 
We will 
always assume our initial 
points (ξ_0,θ_0) 
are chosen so that θ_0 
satisfies (<ref>) 
for n = 0, so that we can skip the first step in 
Algorithm <ref>.














 §.§ Example: linear flow


Consider a PDMP with deterministic paths that are lines in ℝ^d-1 corresponding 
to a finite collection of velocity vectors d_i ∈ℝ^d-1, i ∈ℐ⊆ℕ. Its generator L is defined on suitable functions f:ℝ^d-1×ℐ→ℝ by

    Lf(x,i) = d_i ·∇ f(x,i) + ∑_j iλ_j(x,i)(f(x,j)-f(x,i)),

where λ_j(x,i)≥ 0 for j  i. Suppose we want to sample the probability density

    Z^-1e^-V(x),     Z = ∫ e^- V(x) dx,

where V:ℝ^d-1→ℝ is 
smooth and grows 
sufficiently fast at ∞ 
so that Z<∞. For the PDMP generated by (<ref>) to have an invariant probability density independent of i and proportional to (<ref>), 
the jump rates must satisfy

    ∑_j i(λ_i(x,j) - λ_j(x,i)) = -  d_i ·∇  V(x).

See Remark <ref> in Section <ref> below, 
and <cit.> for a 
similar calculation.

Event 
Chain Monte Carlo and the Zig-Zag process fit into this 
framework. And while these methods 
were designed for efficient sampling, 
we argue that they may be 
limited in certain situations. 
To see why, note that (<ref>) says that at x, the rate into d_i minus the rate out of d_i equals minus 
the gradient of V in direction d_i. Thus the PDMP is likely to change 
directions when it moves up a steep slope of V. This suggests the PDMP 
can struggle to escape from a 
basin of attraction of V, 
defined as the set 
of initial conditions x(0) for which dx(t)/dt = -∇ V(x(t)) has a unique long-time limit.




§ PARREP FOR PDMPS


In this section Z(t)_t ≥ 0 is 
a PDMP with stationary 
distribution π, and 
f is a real-valued function 
defined on the state space 
of Z(t)_t ≥ 0.
Below we outline some 
ParRep algorithms 
for estimating coarse dynamics 
as well as stationary averages, 
with a focus on the latter. The 
stationary average of f is 

    ⟨ f ⟩ = ∫ f(x)π(dx).


Algorithms <ref>  and <ref> 
below are ParRep algorithms based on the skeleton chain 
and the continuous time PDMP, respectively. 
Algorithms <ref> and <ref> 
are parallel steps for synchronous 
computing,
while 
Algorithms <ref> and <ref> 
are for asynchronous computing. 
Algorithms <ref> and <ref>, 
which are essentially 
extensions to PDMPs of 
algorithms recently proposed 
for continuous time Markov chains <cit.>, use the ordering 
of trajectory fragments defined in 
Proposition <ref>. 
Algorithms <ref> and <ref> employ the 
wall-clock time ordering of 
fragments from Proposition <ref>. 
We prove consistency 
of all of our parallel steps 
via our trajectory 
fragment framework.



We do not 
attempt to prove existence, uniqueness 
or convergence to the QSD for general PDMPs. Instead we 
refer the reader 
to recent articles <cit.> 
for conditions which 
ensure convergence 
to a unique QSD. 
From those works,
under appropriate assumptions,
one can establish  convergence 
to a unique 
QSD in D ×ℐ⊆ℝ^d
for a PDMP generated by (<ref>).
For 
instance, exponential 
convergence 
is guaranteed if D⊆ℝ^d-1 is an open
connected bounded domain 
and there exist 
m,M so that 0<m≤λ_j(x,i)≤ M for 
all x ∈ D and i,j ∈ℐ: see <cit.>, pg. 261.
Similar arguments can be made for the 
QSD of the skeleton chain.
Even without theoretical 
guarantees, in practice, 
one can empirically 
validate convergence 
to the QSD using 
certain diagnostics; 
see for instance <cit.>. 




 §.§ Skeleton chain-based parrep algorithm


Let 𝒲 be the 
collection of metastable sets for 
the skeleton chain 
(ξ_n,θ_n)_n ≥ 0. For 
instance, if
Z(t)_t ≥ 0 has generator 
similar to the form (<ref>) 
and we want to sample from 
the distribution (<ref>),
it is natural to define 𝒲 
in terms of basins of attraction of V, 
in which case elements of 𝒲 may be identified on the 
fly by gradient descent <cit.>. See Section <ref> for an example of metastable sets defined this way.

(ξ_n,θ_n)_n ≥ 0
has a QSD ν = ν_W in each W∈𝒲 satisfying 

    ν(A) = lim_n →∞ℙ((ξ_n,θ_n) ∈ A| (ξ_m,θ_m) ∈ W, 0≤ m ≤ n)   ∀ A ⊆ W.

 


For simpler notation, we
do not explicitly indicate the dependence of ν on W. 











Algorithms <ref> and <ref> are parallel steps 
designed for synchronous 
and asynchronous computing, respectively. See Figure <ref> 
for a diagram of both parallel steps. 
The first step in Algorithm <ref> and 
Algorithm <ref> – 
called dephasing in the literature <cit.> – involves generating 
R independent samples from the QSD ν in W. 
These QSD samples may be 
obtained in a variety of ways. One option is to do rejection sampling 
using independent copies of the skeleton chain: whenever a copy 
escapes from W, start it afresh in W until each
copy has remained in W for a long enough consecutive time. Another possibility 
is based on the Fleming-Viot branching process <cit.>: when a copy escapes from W, 
restart it at the current position of a 
copy still in W chosen at random. For more discussion see <cit.>. 







Recall the speedup from ParRep comes from 
the parallel step. The speedup – 
the factor by which ParRep reduces the wall-clock computation time, compared to 
serial simulation of a trajectory 
of the same physical time – can be a factor of up to R, the number of copies or replicas <cit.>, 
when Algorithm <ref> is used to simulate the skeleton chain. See Figure <ref>.
The parallel step is consistent no matter the 
choice of W, but if W is not 
metastable there 
may be no gain in efficiency, as too much computation time 
will be spent sampling the QSD.






Theorem <ref> 
gives conditions that establish 
consistency of  Algorithm <ref> and <ref>. 
For Algorithm <ref>,
the crucial condition is that the wall-clock times it takes 
for the processors to compute steps of the 
skeleton chains are independent of those 
chains. Whether this holds true will depend on the algorithm used to simulate the PDMP. 
If it is a time discretization-based algorithm, or an implementation 
of Algorithm <ref> based 
on Poisson thinning, then the computational effort to 
obtain one step of the skeleton 
chain can be larger in regions 
in state space with lower jump rates.
For CTMCs simulated via the SSA/Gillespie  
algorithm <cit.>, the effort 
to simulate one step of the 
skeleton chain may be
essentially independent 
of the position of the chain. 


150pt
2pt
(i) Let ν be the QSD of (ξ_n,θ_n)_n ≥ 0 in some W∈𝒲, suppose that ℒ(ξ_0,θ_0) = ν, and define M = inf{n≥ 0:(ξ_n,θ_n) ∉ W}. Then in Algorithm <ref>,

    𝔼(f_par) = 𝔼(∑_n=0^M-1∫_0^θ_nf(ψ(t,ξ_n)) dt )
 
and

    (ξ_par,θ_par) ℒ= (ξ_M,θ_M).

5pt
(ii) 
Suppose  
(t_wall(m)_m ≥ 0)^1 ≤ r ≤ R, the wall-clock times from Algorithm <ref>, satisfy
the assumptions 
of Proposition <ref> when (Y^r(t)_t ≥ 0)^r=1,…,R 
equals ((ξ_n^r,θ_n^r)_n ≥ 0)^r=1,…,R. Adopt the 
assumptions in (i) above. Then (<ref>)-(<ref>) hold.


The times T_corr^ν(W) in Algorithm <ref> may be 
 chosen on the fly, or they may be 
 set at the beginning of simulations. Choosing an appropriate value may be done using various
 convergence diagnostics or 
 a priori information; see <cit.> for details.


Consistency of 
the parallel steps, 
together with exactness 
of the decorrelation step,
show that Algorithm <ref> 
produces correct stationary 
averages, provided 
some mild recurrence 
assumptions hold <cit.>. The 
reason is 
essentially the law of 
large numbers: for 
computations of stationary 
averages, due 
to repeated visits 
to each metastable set, 
in the parallel steps
it is enough to get contributions
f_par 
to f_sim with
the correct average value 
along with
escape events with 
the correct law.

We do not attempt here to prove 
ergodicity using this argument, 
but mention it
has been studied previously 
in <cit.>. 
Our numerical simulations in 
Section <ref> below 
also support its validity. One interesting aspect of the 
parallel step is that  
the averaging over independent copies or replicas  
can be considered a bonus, as it likely lowers 
the variance of the estimate 
f_sim/T_sim≈⟨ f ⟩ of the stationary average, compared to an estimate
from a serial trajectory of physical 
time length T_sim.







 §.§ Continuous time PDMP-based algorithm


Let 𝒱 be the 
collection of metastable sets for 
Z(t)_t ≥ 0. As above, if
Z(t)_t ≥ 0 has a generator 
similar to (<ref>) 
and we want to sample 
from the distribution (<ref>), 
the elements of 𝒱 
can be defined in terms of the basins of attraction of V. 
We will require a 
time interval 
Δ t>0, which 
is not necessarily 
a time step for discretizing the PDMP. 
For instance, Δ t could 
be a polling time for resynchronizing parallel processors. 





 






We will adopt the following assumption.

Z(t)_t ≥ 0
has a QSD μ = μ_W in each W∈𝒱 
satisfying 

    μ(A) = lim_t →∞ℙ(Z(t) ∈ A|Z(s) ∈ W, 0 ≤ s ≤ t)  ∀ A ⊆ W.

 

We do not explicitly indicate the dependence of μ on W. 
Notice the QSD of the PDMP is different from that of its skeleton chain in general.



Algorithms <ref> and <ref> are parallel steps 
designed for synchronous 
and asynchronous computing, respectively; see Figure <ref>.
The first step in Algorithm <ref> and 
Algorithm <ref> – the dephasing step – involves generating 
R independent samples from the QSD μ in W. 
Note that this is the QSD of the PDMP in W, 
not the QSD of its skeleton chain. The 
QSD samples can be obtained exactly as 
described in the previous section, but with the PDMP taking 
the place of the skeleton chain.







The speedup from the parallel 
step can be up to a factor of R, the number of copies or replicas <cit.>, 
provided the underlying PDMP 
simulation algorithm is based on time discretization. If the PDMP simulation algorithm is based on computing the skeleton chain, then 
the speedup in Algorithm <ref> 
may be reduced. 
This can be mitigated, however,
by using Algorithm <ref> instead.
The parallel step is consistent for any W, with a speedup if W is metastable for the PDMP.

Theorem <ref> 
gives conditions that establish 
consistency of the parallel steps,
Algorithm <ref> and <ref>. 
For the asynchronous parallel step, 
the crucial condition essentially 
says that the wall-clock time it takes to compute 
a Δ t time interval of Z(t)_t≥ 0 is independent of its position. This is reasonable if Z(t)_t≥ 0 is 
simulated via a time discretization  
technique with a fixed time step. It may not be reasonable 
if a skeleton chain-based technique, 
like Algorithm <ref>, is used instead.






150pt
2pt
(i) Let μ be the QSD of Z(t)_t ≥ 0 in some W ∈𝒱, 
suppose that ℒ(Z(0)) = μ, 
and define τ = inf{t≥ 0:Z(t) ∉ W}. Then in Algorithm <ref>,

    𝔼(f_par) = 𝔼(∫_0^τ f(Z(t)) dt)
 
and

    (T_par,Z_par) ℒ= (τ,Z(τ)).

5pt
(ii) 
Suppose 
(t_wall(m)_m ≥ 0)^1 ≤ r ≤ R, the wall-clock times from Algorithm <ref>, satisfy
the assumptions 
of Proposition <ref> when (Y^r(t)_t ≥ 0)^r=1,…,R 
equals (Z^r(t)_t ≥ 0)^r=1,…,R. Adopt the 
assumptions in (i) above. Then (<ref>)-(<ref>) hold.



Algorithm <ref> generates correct stationary 
averages by the same 
argument as in the previous 
section. 
It is worth mentioning that
Algorithms <ref> 
and <ref>
have a property not shared by 
Algorithms <ref> and <ref>: the 
escape events (T_par,Z_par) in these 
parallel steps
have the correct law for the PDMP. 
This allows us to use Algorithm <ref> to 
compute the dynamics of Z(t)_t≥ 0. 
More precisely, Algorithm <ref> 
leads to a PDMP dynamics 
that is correct on the 
quotient space obtained 
by considering each W ∈𝒱 
as a single point. Note 
that Algorithm <ref> 
cannot be used in this 
way, as it generates 
dynamics of the skeleton 
chain and not the PDMP.

 




§ NUMERICS





Here we test our algorithms 
above on a toy 
PDMP model, our aim 
being to illustrate  
Algorithms <ref> 
and <ref>. 
We will use these algorithms 
to sample the stationary 
average of a function f with 
respect to the Boltzmann 
density
π = Z^-1e^-β V, 
where f and V are defined 
below and β > 0 
is inverse temperature. The 
toy model is a two-dimensional 
version of a PDMP that may be 
defined in 
an arbitrary dimension d-1, as 
follows. Let d_0,…,d_N-1∈ℝ^d-1 be direction vectors such that 
d_0 + … + d_N-1 = 0. Let 
ℤ_N denote the integers 
modulo N, consider the indices of the d_k's as elements of ℤ_N, and for k,ℓ∈ℤ_N 
define

    F_k,ℓ(x)    = β (d_k + … d_k+ℓ)·∇ V(x).



Consider the PDMP with generator defined by

    Lg(x,k) = d_k ·∇ g(x,k) 
    + [g(x,k-1)-g(x,k)]max_0 ≤ℓ≤ N-1F_k,ℓ(x)

for suitable 
g:Ω×ℤ_N →ℝ 
where either Ω = ℝ^d-1 
or Ω is a cube in ℝ^d-1
with periodic boundaries. This is the generator 
for a PDMP that, when moving in direction d_k at point x, switches to direction 
d_k-1 with rate max_0 ≤ℓ≤ N-1F_k,ℓ(x). The 
resulting process can be seen 
as a rejection-free or “lifted” 
version of the sequential Metropolis 
algorithm <cit.>, historically the first 
nonreversible sampling algorithm <cit.> for sampling the Boltzmann distribution. Straightforward calculations show this PDMP has invariant density proportional 
to π; 
see Remark <ref> 
in Section <ref>.

We consider the case where 
state space is Ω = [0,1]^2 with 
periodic boundaries, N = 4, 



d_0 = (1,0), d_1 = (-1,0),  
d_2 = (0,1), d_3 = (0,-1), 
and
the potential energy
V is
pictured in Figure <ref>.
Specifically 

    V(x,y) = cos(4π x)+cos(4π y)+1/5sin(2π x)+1/5sin(2π y).

We define
𝒲 and 𝒱 
using the basins of attraction 
W_i, i=1,…,4 defined 
as
the four squares of 
equal side length 1/2 inside [0,1]^2. See Figure <ref>. 
Thus with x and k the position and direction variables, respectively, of the skeleton chain and PDMP, and θ the jump time variable of the skeleton chain,

    𝒲   ={{(x,k,θ):x ∈ W_i, k ∈{0,1,2,3}, θ>0}: i=1,2,3,4},
    𝒱   = {{(x,k):x ∈ W_i, k ∈{0,1,2,3}}:i=1,2,3,4}.

That is, the skeleton chain or PDMP is in a given set in 𝒲 or 
𝒱 at a 
particular time if and only if its position 
variable belongs to a given 
W_i at that time.



We tested Algorithm <ref> 
and <ref> with 
the synchronous parallel steps 
Algorithm <ref> 
and Algorithm <ref>, 
respectively. We used 
both algorithms to estimate the stationary 
average
⟨ f ⟩ where f(x,k) = 1_x ∈ W_1, 
the characteristic function of the deepest basin of V.
We used up to R = 100 replicas and decorrelation 
times that 
were the same in each basin, T_corr^ν≡ T_corr^ν(W_i) and T_corr^μ≡ T_corr^μ(W_i), i=1,…,4. We 
used Algorithm <ref>
 with time step δ t=10^-2
to simulate the PDMP. 
In Algorithm <ref> we took 
Δ t = δ t =10^-2.
The results 
are in Figures <ref>, <ref>
and <ref>.





To analyze our results, we 
defined an idealized speedup factor as follows. 
Let T^R be an idealized 
wall-clock time for a 
simulation of Algorithm <ref> or <ref> using the parallel steps 
Algorithm <ref> 
and <ref>, respectively, up to a 
fixed time T_stop. The 
idealized wall-clock time T^R is
obtained by  
assuming that we use R parallel processors with zero 
communication cost, such that
on each processor, one 
step of the skeleton chain 
is computed in wall clock 
time 1. Writing T^1 for the wall-clock time corresponding to  
1 processor or direct serial simulation, we define

    time speedup = T^R/T^1.

We include computation time from the 
dephasing step – i.e. the QSD sampling 
step in Algorithms <ref> 
and <ref> – as 
part of T^R. We assume 
this dephasing is done using 
a Fleming-Viot-based technique 
as described above, with R copies 
of the underlying skeleton chain or PDMP. 




Processor communication, 
which we do not account 
for, 
of course takes a toll on the time 
speedup. However, 
the processor communication cost 
is small compared 
to the rest of the computational effort
if the sets in 𝒲 
and 𝒱 are 
significantly metastable. Thus, our time speedup 
gives a reasonable picture of 
the gain that can be expected.


The time speedup depends on the parameters in Algorithms <ref> and <ref>. 
If all other parameters are held constant,
the time speedup 
increases with R or β, due 
to increasing parallelization or metastability, 
respectively (Figure <ref>), while the 
time speedup decreases 
with T_corr^ν and 
T_corr^μ, due to increased 
effort to sample the QSD (Figure <ref>). 
With increasing metastability, 
the relative computational effort to sample 
the QSD decreases in comparison with 
the effort to simulate an 
escape from a metastable set. 
Since the latter is done in 
parallel, increasing the metastability 
leads to a larger time speedup. 








Figure <ref> shows 
that the approximation 
f_sim/T_sim approaches 
the stationary average 
⟨ f⟩ as the 
decorrelation times
T_corr^ν and T_corr^μ 
increase, as expected. As 
discussed above, appropriate 
values of these QSD sampling times 
depend on the degree 
of metastability.
Note that the approximations
f_sim/T_sim are quite 
good even for small QSD sampling times. 
This was true not just for f(x,k) = 1_x ∈ W_1 
but for a
variety of other functions. This 
feature, of reasonable accuracy in ParRep
even for relatively small decorrelation 
times, was observed before 
in <cit.>. 
Here this may be a result of the momentum-like direction variables, which 
can make the PDMP unlikely to immediately 
escape from a metastable set just 
after entering.











§ PROOFS



Our first two results 
below, Proposition <ref> and <ref>, 
establish the memoryless 
distribution of the escape 
time starting from the QSD, and the independence of 
the escape time and escape 
point, for general 
Markov processes. See for instance <cit.> for details. 
We include proofs 
for completeness.


Let X(t)_t ≥ 0 have a
QSD ρ in U, and suppose 
ℒ(X(0)) = ρ. Suppose
T= inf{t≥ 0: X(t) ∉ U} 
is finite almost surely. 
Then T
has a memoryless distribution; 
that is, T is either  
exponentially or geometrically distributed.



The definition (<ref>) of the QSD 
together with the Markov property show that
ℒ(X(t+s)_s ≥ 0 | T > t) = ℒ(X(s)_s ≥ 0).
Thus, ℙ(T > t+s|T> t) = ℙ(T > s) for any s,t≥ 0.
When T is finite-valued, the only distribution satisfying this 
is

    ℙ(T>t) = e^-λ t.

We use (<ref>) to indicate either the (continuous) exponential distribution 
with parameter λ > 0, 
or the (discrete) geometric distribution 
with parameter p = 1-e^-λ.

 
 
 Let X(t)_t ≥ 0 have a 
 QSD ρ in U, and suppose 
 ℒ(X(0)) = ρ. 
Suppose T = inf{t ≥ 0:X(t) ∉ U} 
is finite almost surely. Then 
 T and X(T) are 
 independent.
 
 
 
 By Proposition <ref>, ℙ(T>t) = e^-λ t. For A in the complement of U,
 
    ℙ(T ∈ ((n-1)t,nt],X(T) ∈ A) 
       = ℙ(T > (n-1)t, 
    X(nt ∧ T) ∈ A) 
       = ℙ(X(nt∧ T) ∈ A|T>(n-1)t)ℙ(T>(n-1)t)
       = ℙ(X(T) ∈ A,T ≤ nt|T>(n-1)t)ℙ(T>(n-1)t)
       = ℙ(X(T)∈ A,T ≤ t)e^-λ (n-1)t,

where the last step uses (<ref>). 
Summing over n ≥ 1 establishes the result:

    ℙ(X(T) ∈ A) 
       = ℙ(X(T) ∈ A, T≤ t)1/1-e^-λ t
       = ℙ(X(T) ∈ A, T≤ t)/ℙ(T ≤ t).


 
 
Below, we write MGF 
for the moment generating 
function of a random variable.
The results in 
Proposition <ref> below 
hold in both continuous and 
discrete time. To connect
the discrete and continuous time cases, 
we write 1-p = e^-λ 
where p ∈ (0,1) is the geometric parameter 
and λ>0 is the exponential rate. 


Let 
t_1,t_2,… be nonnegative deterministic times such that ∑_m=1^∞ t_m = ∞. Let τ_1,τ_2,… be 
random variables such that ℙ(τ_1>t) =  e^-λ t and 

    ℙ(τ_m>t|τ_m-1>t_m-1,…,τ_1 >t_1) = e^-λ t,    m≥ 2.
 
Let 
L= inf{m ≥ 1:τ_m ≤ t_m}.
Then 
    ℙ(t_1+… + t_L-1+τ_L>t)= e^-λ t.




Let s_m = t_1+… +t_m and s_0 = 0. 
By (<ref>) and induction,

    ℙ(τ_m>t_m,…,τ_1>t_1)= e^-λ (t_1+… + t_m) = e^-λ s_m.

Using (<ref>) again, in the continuous case,

    𝔼(e^uτ_m1_τ_m ≤ t_m|τ_m-1>t_m-1,…,τ_1>t_1) 
       = ∫_0^t_m e^usλ e^-λ s ds =  λ/u-λ(e^(u-λ)t_m-1),

while in the discrete case, where e^-λ = 1-p, 

    𝔼(e^uτ_m1_τ_m ≤ t_m|τ_m-1>t_m-1,…,τ_1>t_1) 
       = ∑_s=1^t_m e^use^-λ (s-1)(1-e^-λ) = e^u-e^u-λ/1-e^u-λ(1-e^(u-λ)t_m).

Note also that 

    {L = m} = {τ_m ≤ t_m,τ_m-1 > t_m-1,…,τ_1 > t_1}.

Consider the continuous case. Combining (<ref>), (<ref>) and (<ref>) gives

    𝔼(e^uτ_m1_L=m)
       = 𝔼(e^uτ_m1_τ_m≤ t_m1_τ_m-1>t_m-1,…,τ_1>t_1)
       = 𝔼(e^uτ_m1_τ_m≤ t_m|τ_m-1>t_m-1,…,τ_1>t_1)ℙ(τ_m-1>t_m-1,…,τ_1>t_1)
       = λ/u-λ(e^(u-λ)t_m-1) e^-λ s_m-1.

We now see that s_L-1+τ_L has the MGF of an exponential(λ) random variable:

    𝔼(e^u(s_L-1+τ_L))    = ∑_m=1^∞𝔼(e^u(s_L-1+τ_L)1_L=m)  
       = ∑_m=1^∞ e^us_m-1𝔼(e^uτ_m1_L=m)
       = λ/u-λ∑_m=1^∞(e^(u-λ)s_m-e^(u-λ)s_m-1)
       = λ/λ - u    if u<λ.

Similarly, in the discrete case, combining (<ref>), (<ref>) and (<ref>) gives

    𝔼(e^uτ_m1_L=m)
       = 𝔼(e^uτ_m1_τ_m≤ t_m1_τ_m-1>t_m-1,…,τ_1>t_1)
       = 𝔼(e^uτ_m1_τ_m≤ t_m|τ_m-1>t_m-1,…,τ_1>t_1)ℙ(τ_m-1>t_m-1,…,τ_1>t_1)
       = e^u-e^u-λ/1-e^u-λ(1-e^(u-λ)t_m) e^-λ s_m-1.

This shows again that s_L-1+τ_L has the MGF of a geometric(p) random variable, via

    𝔼(e^u(s_L-1+τ_L))    = ∑_m=1^∞𝔼(e^u(s_L-1+τ_L)1_L=m)  
       = ∑_m=1^∞ e^us_m-1𝔼(e^uτ_m1_L=m)
       = e^u-e^u-λ/1-e^u-λ∑_m=1^∞ (e^(u-λ)s_m-1-e^(u-λ)s_m) 
       =e^u-e^u-λ/1-e^u-λ    if u<λ
       = pe^u/1-(1-p)e^u     if u< -log(1-p).





Let A be a subset of the complement of U. Note that, due to (<ref>), 
the events {X_m(T_m) ∈ A,T_m≤ t} and
{T_m-1>t_m-1,…,T_1>t_1} are independent conditional on 
{X_m(0)= x}. Using this, (<ref>) and 
Proposition <ref>,

    ℙ(X_m(T_m) ∈ A,T_m ≤ t|T_m-1 > t_m-1,…,T_1 > t_1) 
       = ∫ℙ(X_m(T_m) ∈ A,T_m ≤  t|X_m(0) = x,T_m-1 > t_m-1,…,T_1 > t_1)
                   ×ℙ(X_m(0) ∈ dx|T_m-1 > t_m-1,…,T_1 > t_1) 
       = ∫ℙ(X_m(T_m) ∈ A,T_m ≤  t|X_m(0) = x)ρ(dx) 
       = ℙ(X(T) ∈ A,T ≤  t) 
       = ℙ(X(T) ∈ A)ℙ(T ≤ t).

Taking A as the complement of U in (<ref>), and using Proposition <ref>, 

    ℙ(T_m > t|T_m-1 > t_m-1,…,T_1 > t_1) =  ℙ(T>t) = e^-λ t,   some λ>0.

By (<ref>), ℙ(T_1>t) = e^-λ t. 
Thus by Proposition <ref>, ℒ(T_par) = ℒ(T). Notice 

    {L = m} = {T_m ≤ t_m,T_m-1>t_m-1,…,T_1>t_1}.

From (<ref>), (<ref>) and (<ref>) we have

    ℙ(X_par∈ A, T_par > t|L=m) 
       =ℙ(X_m(T_m) ∈ A, s_m-1+T_m > t|L=m)
       =ℙ(.X_m(T_m) ∈ A, T_m > t-s_m-1|T_m ≤ t_m,T_m-1 > t_m-1,…,T_1 > t_1)
       =  ℙ(.X_m(T_m) ∈ A,T_m∈ (t-s_m-1,t_m] |T_m-1 > t_m-1,…,T_1 > t_1)
                   ×ℙ(T_m ≤ t_m|T_m-1 > t_m-1,…,T_1 > t_1)^-1
       =ℙ(X(T) ∈ A)ℙ(T ∈(t-s_m-1,t_m])ℙ(T ≤ t_m)^-1
       = ℙ(X(T) ∈ A)ℙ(T>t-s_m-1|T≤ t_m).
 
From (<ref>) we conclude

    ℙ(X_par∈ A, T_par > t) 
       = ∑_m=1^∞ℙ(X_par∈ A,T_par>t|L=m)ℙ(L=m)
       = ℙ(X(T) ∈ A)∑_m=1^∞ℙ(T > t-s_m-1|T≤ t_m)ℙ(L=m).

In the last display, taking t =0 shows ℙ(X_par∈ A)= ℙ(X(T) ∈ A), while taking A as the complement of U shows ℙ(T_par>t) = ∑_m=1^∞ℙ(T > t-s_m-1|T≤ t_m)ℙ(L=m). Thus, (<ref>) 
shows that ℒ(X_par) = ℒ(X(T)) 
and X_par, T_par are independent. 




We consider only the discrete time case, since 
the arguments in the continuous time case are 
analogous. 
Let r ∈ℕ be fixed. Observe that

    ∑_t=0^r-1ℙ(T>t) = 
    ∑_t=0^r-1𝔼(1_T>t) 
    = ∑_t=0^∞𝔼(1_T∧ r>t) 
    = 𝔼(T ∧ r).

By the preceding display and the 
the definition (<ref>) of 
the QSD ρ,

    𝔼(∑_t=0^T∧ r-1g(X(t)))    = 
    ∑_s=1^∞∑_t=0^s-1𝔼(g(X(t))1_T ∧ r = s) 
       = ∑_t=0^∞𝔼(g(X(t))1_T ∧ r > t) 
       = ∑_t=0^r-1𝔼(g(X(t))1_T>t) 
       = ∑_t=0^r-1𝔼(g(X(t))|T>t)ℙ(T>t)
       = (∫ g dρ)𝔼(T ∧ r).

Since 
{L ≥ m} = {T_m-1>t_m-1,…,T_1>t_1},
using (<ref>), (<ref>) and (<ref>) we get 

    𝔼(.∑_t=0^T_m ∧ t_m-1 g(X_m(t))|L ≥ m) 
       = ∫𝔼(.∑_t=0^T_m ∧ t_m-1 g(X_m(t))|X_m(0) = x,L≥ m) ℙ(X_m(0) ∈ dx|L ≥ m)
       = ∫𝔼(.∑_t=0^T_m ∧ t_m-1 g(X_m(t))|X_m(0) = x)ρ(dx)
       = 𝔼(∑_t=0^T ∧ t_m-1 g(X(t))) = (∫ g dρ)𝔼(T∧ t_m)
       = (∫ g dρ)𝔼(T_m ∧ t_m|L ≥ m)

where the last line of (<ref>) follows from taking g ≡1 in the first four lines of (<ref>). By Theorem <ref>, ℒ(T) = ℒ(T_par)=  ℒ(∑_m=1^L T_m ∧ t_m) and thus

    𝔼(T)    = 𝔼(∑_m=1^L T_m ∧ t_m)
       = ∑_n=1^∞∑_m=1^n𝔼(T_m ∧ t_m1_L=n) 
       = ∑_m=1^∞𝔼(T_m ∧ t_m 1_L ≥ m) 
       = ∑_m=1^∞𝔼(T_m ∧ t_m |L ≥ m)ℙ(L ≥ m).

Now by (<ref>) and (<ref>),

    𝔼(∑_m=1^L∑_t=0^T_m ∧ t_m-1 g(X_m(t)))    = 
    ∑_n=1^∞∑_m=1^n𝔼(1_L=n∑_t=0^T_m ∧ t_m-1 g(X_m(t))) 
       =  ∑_m=1^∞𝔼(1_L≥ m∑_t=0^T_m ∧ t_m-1 g(X_m(t))) 
       = ∑_m=1^∞𝔼(.∑_t=0^T_m ∧ t_m-1 g(X_m(t))|L≥ m)ℙ(L≥ m) 
       = (∫ g dρ)∑_m=1^∞𝔼(T_m ∧ t_m|L ≥ m)ℙ(L≥ m) 
       = (∫ g dρ)𝔼(T).

Letting r →∞ in (<ref>), using dominated convergence, and 
comparing with (<ref>), 
    𝔼(∑_m=1^L∑_t=0^T_m ∧ t_m-1 g(X_m(t))) = 𝔼(∑_t=0^T-1g(X(t)))

as desired.


Below we will need 
the following basic facts. 

Let ℱ, 
𝒢,
ℋ, and 𝒦 be σ-algebras.

    
   (i) Let σ(𝒢,ℋ) be the σ-algebra 
    generated by 𝒢 and ℋ. Suppose that 
    ℱ,
    σ(𝒢,ℋ) 
    are independent conditional on 𝒦, and that 𝒢, ℋ are independent conditional on 𝒦. 
    Then ℱ, 𝒢, ℋ are mutually 
    independent conditional on 𝒦.
    
   (ii) Suppose ℋ⊆𝒢. If ℱ 
    and 𝒢 are independent, 
    then ℱ 
    and 𝒢 are independent 
    conditional on ℋ.



Throughout let A ∈ℱ, 
B ∈𝒢, and 
C ∈ℋ, and write 
    1_S for the 
    characteristic or indicator function of a set S. Consider (i). Since A ∈ℱ, B ∩ C ∈σ(𝒢,ℋ), and ℱ,  σ(𝒢,ℋ)
 are independent conditional on 𝒦, 𝔼(1_A 1_B ∩ C|𝒦) 
    = 
    𝔼(1_A|𝒦) 𝔼(1_B ∩ C|𝒦) almost surely. Similarly 𝔼(1_B1_C|𝒦)
    = 𝔼(1_B|𝒦) 
    𝔼(1_C|𝒦) 
    almost surely. Thus,

    𝔼(1_A 1_B1_C|𝒦) = 
        𝔼(1_A 1_B ∩ C|𝒦) 
           = 
        𝔼(1_A|𝒦) 𝔼(1_B ∩ C|𝒦) 
       = 𝔼(1_A|𝒦) 𝔼(1_B1_C|𝒦)
        =𝔼(1_A|𝒦) 𝔼(1_B|𝒦) 
        𝔼(1_C|𝒦)

    almost surely, which proves (i). 
    
    
    Consider now (ii). Define P = 𝔼(1_A1_B|ℋ) and Q = 𝔼(1_A|ℋ)𝔼(1_B|ℋ). 
   As P and Q are ℋ-measurable and C ∈ℋ is arbitrary, if 𝔼(P1_C) = 𝔼(Q1_C) then 
   we can use uniqueness of conditional expectation to conclude P = Q almost surely, 
   so that (ii) holds. Since ℋ⊆𝒢,  B ∩ C ∈𝒢.
    Moreover A ∈ℱ 
    and ℱ,𝒢 
    are independent, so
    
    𝔼(P1_C) = 𝔼(1_A 1_B1_C) =ℙ(A ∩ B ∩ C) = ℙ(A)ℙ(B ∩ C).

    The first equality in (<ref>) comes from definition of conditional expectation. Since A ∈ℱ 
    and ℱ,𝒢 are independent, 𝔼(1_A|𝒢) = 𝔼(1_A). 
    So by the tower property,
    
    𝔼(1_A|ℋ)
            = 𝔼(𝔼(1_A|𝒢)|ℋ) = 
            𝔼(𝔼(1_A)|ℋ)
            =𝔼(1_A)=
            ℙ(A).

    Moreover since 1_C is ℋ-measurable, 𝔼(1_B|ℋ)1_C = 𝔼(1_B1_C|ℋ). Thus,
    
    𝔼(Q1_C)    = 𝔼(𝔼(1_A|ℋ)
            𝔼(1_B|ℋ)1_C)
       =𝔼(ℙ(A)
            𝔼(1_B1_C|ℋ)) 
       = ℙ(A)𝔼(𝔼(1_B∩ C|ℋ))
       = ℙ(A)𝔼(1_B∩ C) = ℙ(A)ℙ(B ∩ C),

    with the last line using the tower property.
    Now (ii) follows from (<ref>)-(<ref>).



Adopt the notation of Assumption <ref>. 
Below let j,k,ℓ denote positive integers. It is easy to check that 
r_j = r_k if and only if k-j is an integer 
multiple of R, while m_k-ℓ R = m_k-ℓ when 
ℓ≤ m_k. Thus, 

    {m_j: j <k, r_j = r_k}   = 
    {m_k-ℓ R: ℓ≤  m_k}
       = {0,1,…,m_k-1},

with both sides empty if m_k = 0.
See Figure <ref>. 
Define 
    ζ^r = inf{t≥ 0:Y^r(t) ∉ U}.

Fix k ≥ 2 and note that, by definition of the fragments,

    {ζ^r_j>(m_j+1)Δ t,   ∀ j <k}⊆{T_k-1>Δ t,…,T_1>Δ t},

where E ⊆ E' indicates the event E' occurs whenever E occurs. By (<ref>), 

    {Y^r_j(t) ∈ U    ∀ t ∈ [m_jΔ t,(m_j+1)Δ t],j≤ k}
       ⊆{Y^r_k(t) ∈ U    ∀  t ∈ [0,(m_k+1)Δ t]}.

By definition of the fragments and (<ref>),

    {T_k-1>Δ t,…,T_1>Δ t}   =  {Y^r_j(t) ∈ U,    ∀ t ∈ [m_jΔ t,(m_j+1)Δ t],j<k}
       ⊆{Y^r_j(t) ∈ U,    ∀ t ∈ [0,(m_j+1)Δ t],j<k}
       = {ζ^r_j>(m_j+1)Δ t,   ∀ j <k}.

Combining (<ref>) and (<ref>), 

    {T_k-1>Δ t,…,T_1>Δ t} = {ζ^r_j>(m_j+1)Δ t,   ∀ j<k}.

Due to 
independence of Y^r(t)_t ≥ 0 over r and Lemma <ref>(ii),

    conditional on {ζ^r_j > (m_j + 1)Δ t,   ∀ j<k    s.t.    r_j = r_k},
       the event {ζ^r_j > (m_j + 1)Δ t ,   ∀ j<k   s.t.   r_j  r_k}
       is independent 
        of Y^r_k(m_kΔ t).

Again using (<ref>),

    {ζ^r_k > m_kΔ t} = {ξ^r_j>(m_j+1)Δ t,    ∀ j <k    s.t.   r_j = r_k}.

Combining (<ref>), (<ref>), and (<ref>), and using (<ref>),

    ℒ(X_k(0)|T_k-1>Δ t,…,T_1>Δ t) 
       = ℒ(Y^r_k(m_kΔ t)|ζ^r_j>(m_j+1)Δ t,   ∀ j<k)
       = ℒ(.Y^r_k(m_kΔ t)|ζ^r_j > (m_j + 1)Δ t,    ∀ j<k    s.t.   r_j = r_k)
       = ℒ(.Y^r_k(m_kΔ t)|ζ^r_k>m_kΔ t)= ρ.

As Y^1(t)_t ≥ 0
is a copy of X(t)_t ≥ 0 
with ℒ(X(0))= ρ, in particular 
ℒ(Y^1(0)) = ρ. Thus, 

    ℒ(X_1(0)) = ℒ(Y^r_1(m_1)) = ℒ(Y^1(0)) = ρ.

We have now established (<ref>) of 
Assumption <ref>. 

Consider now (<ref>). Let k ≥ 1.
Due to independence of Y^r(t)_t ≥ 0 over r and Lemma <ref>(ii), 
we see that 
conditional on X_k(0),  
(X_ℓ(t)_0 ≤ t ≤Δ t)_r_ℓ = r_k is independent 
of (X_ℓ(t)_0 ≤ t ≤Δ t)_r_ℓ r_k. 
For k≥ 2, the Markov property of 
Y^r_k(t)_t ≥ 0 
and (<ref>)
show that,
conditional on X_k(0), 
X_k(t)_0 ≤ t ≤Δ t is 
independent of (X_ℓ(t)_0 ≤ t ≤Δ t)_ℓ < k, r_ℓ = r_k. By Lemma <ref>(i) with ℱ = σ((X_ℓ(t)_0 ≤ t ≤Δ t)_ℓ<k, r_ℓ r_k), 
𝒢 =σ((X_ℓ(t)_0 ≤ t ≤Δ t)_ℓ<k, r_ℓ= r_k), 
ℋ = σ(X_k(t)_0 ≤ t ≤Δ t) and 𝒦 = σ(X_k(0)), we have,
for k≥ 2,

    conditional on X_k(0), 
     X_k(t)_0 ≤ t ≤Δ t is 
    independent of (X_ℓ(t)_0 ≤ t ≤Δ t)_ℓ < k.


Now define the fragments' irrelevant 
futures as follows. Let X_k(t)_t ≥Δ t be 
copies of X(t)_t ≥ 0 that evolve 
forward of time independently 
of everything else. That is, for each 
k ≥ 1, conditional on X_k(Δ t), 
X_k(t)_t ≥Δ t, X_k(t)_0 ≤ t ≤Δ t, and (X_ℓ(t)_t ≥ 0)_ℓ<k 
are mutually independent. From (<ref>) 
it is easy to see this is possible, 
as the irrelevant futures have no 
bearing on the definitions of the 
fragments. Now by construction of 
the irrelevant futures and (<ref>), 
it is easy to see that for k≥ 2,
conditional on X_k(0), 
X_k(t)_t ≥ 0 is 
independent of (X_ℓ(t)_t ≥ 0)_ℓ<k. 
This proves (<ref>) 
in Assumption <ref>.








Adopt the notation of Assumption <ref>. 
This proof will follow 
the same basic steps as 
the proof of Proposition <ref>, 
but the justifications will be 
different. Let ζ^r = inf{t≥ 0:Y^r(t) ∉ U} be as above. 

Fix k ≥ 2. We first 
claim that (<ref>) still 
holds. Let 
n ∈{0,1,…,m_k-1}. 
By the surjectivity assumption 
in (iii) there is 
j such that r_j = r_k 
and m_j = n. Since m_j = n < m_k 
and r_j = r_k, from (ii) 
we have t_wall^r_j(m_j) 
≤ t_wall^r_k(m_k). Since 
j  k, using monotonicity 
in (iii) we conclude 
j <k. 
Thus {0,1,…,m_k-1}⊆{m_j:j <k, r_j = r_k}. 
Now consider m_j such 
that j <k and 
r_j = r_k. By monotonicity 
in (iii) we must have 
t_wall^r_j(m_j)<t_wall^r_k(m_k). Then by (ii) we can conclude 
m_j < m_k. Thus 
{m_j:j <k, r_j = r_k} 
= {0,1,…,m_k-1}.

Next we establish (<ref>).
Equipped with (<ref>), 
we see that (<ref>) holds. 
Moreover, since (<ref>) 
agrees with (<ref>), 
the same steps as 
in the Proof of Proposition <ref> show 
that (<ref>) holds. 
On the other 
hand, (<ref>) 
holds because of (i), Lemma <ref>(ii),  
and independence
of Y^r(t)_t ≥ 0 over r. 
The sequence of equalities 
in (<ref>) 
then holds, with the 
last equality using 
(i) again. It 
remains to show that ℒ(X_1(0)) = 
ρ. Suppose m_1 > 0. By surjectivity in 
(iii) there is j>1 such 
that m_j = 0 and r_j = r_1. 
But then (ii) implies t_wall^r_j(m_j) ≤ t_wall^r_1(m_1), which 
contradicts monotonicity in 
(iii). Thus m_1 = 0, so we can apply (i) 
to conclude ℒ(X_1(0)) = ℒ(Y^r_1(m_1)) = ρ.  
Thus (<ref>) in in Assumption <ref> holds.

Consider now (<ref>). 
By (i) and independence of Y^r(t)_t ≥ 0 over r, 
conditional on X_k(0), 
(X_ℓ(t)_0 ≤ t ≤Δ t)_r_ℓ = r_k is independent 
of (X_ℓ(t)_0 ≤ t ≤Δ t)_r_ℓ r_k. 
Recall that (<ref>) still holds. Thus for k≥ 2, by the Markov property of 
Y^r_k(t)_t ≥ 0 and (<ref>), 
conditional on X_k(0), 
X_k(t)_0 ≤ t ≤Δ t is 
independent of (X_ℓ(t)_0 ≤ t ≤Δ t)_ℓ < k, r_ℓ = r_k. 
By Lemma <ref>(i) 
we conclude (<ref>) holds for k ≥ 2. 
Let the trajectory fragments' irrelevant futures 
be independent of everything else 
as in the proof of Proposition <ref>. 
Following the reasoning in that proof
we see that (<ref>) 
in Assumption <ref> holds.








The statements (i) 
and (ii) follow from Propositions <ref> 
and <ref>, 
respectively, with
((ξ_n^r,θ_n^r)_n ≥ 0)^r=1,…,R 
taking the place of 
(Y^r(t)_t ≥ 0)^r=1,…,R, 
and
with t_m ≡Δ t = 1 and g(ξ,θ) = ∫_0^θ f(ψ(t,ξ)) dt.



The statements (i) and (ii)  follow from Propositions <ref> 
and <ref>, 
respectively, with (Z^r(t)_t ≥ 0)^r=1,…,R 
taking the place of 
(Y^r(t)_t ≥ 0)^r=1,…,R, 
and with t_m ≡Δ t and g = f.







 §.§ Supplementary results


We first show that the decoupling 
of wall-clock times from 
the speed of computing 
X(t)_t ≥ 0 is a necessary 
condition for consistency. 
Below we
break assumption 
(i) in Proposition <ref> by assuming the 
wall-clock times to 
obtain the
initial QSD samples 
Y^r(0), r=1,…,R,
in the parallel step 
are correlated with 
the positions of those 
samples.

In Proposition <ref>, 
if (i) does 
not hold,
then the conclusions 
of Theorem <ref> 
and Theorem <ref> 
may not hold. 


Let X(t)_t ≥ 0 be a simple 
random walk on ℤ, 
meaning X(t+1)-X(t) = 1 or -1, 
each with probability 1/2. 
Let U = {0,1}. 
The QSD ρ of X(t)_t ≥ 0 in 
U is simply the 
uniform distribution on U. Assume 
X(t)_t ≥ 0 has initial 
distribution ℒ(X(0)) = ρ, 
and let (Y^r(t)_t ≥ 0)^r=1,…,R 
be independent copies of X(t)_t ≥ 0. Suppose 

    t_wall^r(0) < t_wall^s(0) 
            whenever 
        Y^r(0) = 0   
    and Y^s(0) = 1.

Notice that (<ref>) violates 
(i) of Proposition <ref>. 
Assume however that (ii)
and (iii) in Proposition <ref> hold. 
Then arguments similar 
to those in the proof 
of Proposition <ref> 
show that {Y^r_1(m_1) = 1} 
= 
{Y^r(0) = 1   ∀ r}. 
Adopt the notation of 
Algorithm <ref>. 
Then by the above and the definition (<ref>) of the fragments, 























    ℙ(T_par = 1,X_par = 2)    = ℙ(Y^r_1(m_1) = 1,Y^r_1(m_1+1) = 2) 
       = 1/2ℙ(Y^r_1(m_1) = 1)
       = 1/2ℙ(Y^r(0) = 1   ∀ r) = 1/2(1/2)^R.

Similarly,

    ℙ(T = 1,X(T)=2) = 
        ℙ(X(0) =1, X(1) = 2) 
        = 1/2ℙ(X(0)=1) =1/4.

Notice when R>1, (<ref>) and (<ref>) show  the conclusion of Theorem <ref> does not 
hold, as
ℙ(T_par = 1,X_par = 2) 
ℙ(T = 1,X(T)=2).
A similar construction shows 
the conclusion of Theorem <ref> can fail 
when (i) does not hold.



The next two results below are formal 
calculations related to claims 
made in the text above. 
These results could be 
made precise using results 
in <cit.>. 
However we stick to formal 
computations for brevity.





Suppose (<ref>) holds for all 
 x ∈ℝ^d-1 and i ∈ℐ. Then e^-V(x) 
is formally invariant for a PDMP generated by (<ref>).
 

Let L be defined as in (<ref>).
We will show that 
    ∑_i∈ℐ∫_ℝ^d-1 Lf(x,i) π(x,i) dx = 0
 
provided (<ref>) holds and
π(x,i)∝ e^-V(x), where
 ∝ 
indicates proportional to.
Write 
    λ_i(x,i) = -∑_j iλ_j(x,i).

With sufficient 
regularity we can integrate by parts to get

    ∑_i ∫ Lf(x,i) π(x,i) dx    ∝∑_i ∫(d_i ·∇ f(x,i) + ∑_j λ_j(x,i)f(x,j))e^- V(x) dx
       = -∑_i ∫ f(x,i) ∇·(d_i e^- V(x))dx 
             
    + ∑_i ∑_j ∫λ_i(x,j)f(x,i)e^- V(x) dx 
       = ∑_i ∫ f(x,i)(d_i ·∇ V(x)  +∑_j λ_i(x,j))e^-  V(x) dx.

Above, all the sums are over ℐ and integrals 
are over ℝ^d-1. If (<ref>) holds, 

    d_i ·∇ V(x)  +∑_j λ_i(x,j)    = 
      d_i ·∇ V(x)  +∑_j i (λ_i(x,j)-λ_j(x,i)) = 0.

Comparing with (<ref>) gives the result.

Note that the calculation in 
Remark <ref> 
shows (<ref>) is a necessary 
condition for (<ref>) to define a 
PDMP with an invariant distribution of the form π(x,i)∝ e^-V(x).





e^-β V(x) 
is formally invariant for a PDMP 
generated by (<ref>). 


Let L be defined as in (<ref>). We will show that 
    ∑_k=0^N-1∫_Ω Lf(x,k) π(x,k) dx = 0
 provided
π(x,k) ∝ e^-β V(x). Recall d_0,…,d_N-1∈ℝ^d-1 sum to 0 and we
consider the
indices of the d_k's as elements of ℤ_N, 
the integers modulo N. Write 

    F_k,ℓ(x) = β(d_k+…+d_k+ℓ)·∇ V(x).

With 
sufficient regularity we can integrate by parts to get

    ∑_k =0^N-1∫_Ω Lg(x,k)π(x,k) dx 
       ∝∑_k =0^N-1∫_Ω(d_k ·∇ g(x,k) + [g(x,k-1)-g(x,k)]max_0 ≤ℓ≤ N-1F_k,ℓ(x))e^-β V(x) dx 
       = -∑_k =0^N-1∫_Ω g(x,k)∇·(d_k e^-β V(x))dx
             + ∑_k =0^N-1∫_Ω g(x,k) (max_0 ≤ℓ≤ N-1 F_k+1,ℓ(x) - max_0 ≤ℓ≤ N-1F_k,ℓ(x))e^-β V(x) dx
       =∑_k =0^N-1∫_Ω  g(x,k)(β d_k ·∇ V(x)+max_0 ≤ℓ≤ N-1 F_k+1,ℓ(x) -  max_0 ≤ℓ≤ N-1F_k,ℓ(x) )e^-β V(x) dx,

where when Ω = ℝ^d-1 we assume V grows 
sufficiently fast at ∞ 
so that we can 
neglect the boundary term 
from the integration by 
parts.
Observe that, because ∑_ℓ=0^N-1 d_k+ℓ = 0 and 
d_k+N = d_k, we have

    {d_k+d_k+1,d_k+d_k+1+d_k+2,…,d_k+…+d_k+N-1,d_k+…+d_k+N}
        = {d_k+d_k+1,d_k+d_k+1+d_k+2,…,d_k+…+d_k+N-1,d_k}
       = {d_k,d_k+d_k+1,…,d_k+N-1}.

It follows that

    β d_k ·∇ V(x)+max_0 ≤ℓ≤ N-1 F_k+1,ℓ(x) -  max_0 ≤ℓ≤ N-1F_k,ℓ(x) 
       = βmax_0 ≤ℓ≤ N-1
        (d_k+…+d_k+ℓ+1)·∇ V(x) - βmax_0 ≤ℓ≤ N-1(d_k+…+d_k+ℓ)·∇ V(x)
       = 0.

This proves the desired result.




e^-β V is invariant for Z(nδ t)_n ≥ 0 
defined in Algorithm <ref>.



Write the acceptance 
probability in Algorithm <ref> as 

    A_k(x) = min_0 ≤ℓ≤ N-1exp(β V(x)-β V(x+d_kδ t+… + d_k+ℓδ t)).

Arguing similarly as in the formal proof of Remark <ref>, since 
∑_k=0^N-1 d_k = 0
we have 

    A_k(x)/A_k+1(x+d_kδ t) = exp(β V(x)-β V(x+d_kδ t)).

Now let π(x,k) ∝ e^-β V(x). Then the last display shows that 

    π(x+d_kδ t,k) = π(x,k)A_k(x)+π(x+d_kδ t,k+1)(1-A_k+1(x+d_kδ t)).

Inspecting Algorithm <ref>, 
we see that (<ref>) 
demonstrates the required result.






§ ACKNOWLEDGEMENTS


The author gratefully 
thanks Peter Christman for 
producing the numerical 
results leading to Figures <ref>, <ref> and <ref> in Section <ref>, 
as well as 
Petr Plecháč, Gideon Simpson, and Ting Wang for helpful conversations. The author 
also gratefully acknowledges support 
from the
National Science Foundation 
via the awards
NSF-DMS-1522398 and 
NSF-DMS-1818726.





siamplain



