Piecewise deterministic generative models

Authors: andrea  bertazzi
dario  shariatian
umut  simsekli
eric  moulines
alain  durmus
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 10:58:43

Abstract: We introduce a novel class of generative models based on piecewise deterministic Markov processes (PDMPs), a family of non-diffusive stochastic processes consisting of deterministic motion and random jumps at random times. Similarly to diffusions, such Markov processes admit time reversals that turn out to be PDMPs as well. We apply this observation to three PDMPs considered in the literature: the Zig-Zag process, Bouncy Particle Sampler, and Randomised Hamiltonian Monte Carlo. For these three particular instances, we show that the jump rates and kernels of the corresponding time reversals admit explicit expressions depending on some conditional densities of the PDMP under consideration before and after a jump. Based on these results, we propose efficient training procedures to learn these characteristics and consider methods to approximately simulate the reverse process. Finally, we provide bounds in the total variation distance between the data distribution and the resulting distribution of our model in the case where the base distribution is the standard $d$-dimensional Gaussian distribution. Promising numerical simulations support further investigations into this class of models.

Main paper: 
 introduction: 
 Diffusion-based generative models [Ho et al., 2020 , Song et al., 2021] have recently achieved stateof-the-art performance in various fields of application [Dhariwal and Nichol, 2021 , Croitoru et al., 2023 , Jeong et al., 2021 , Kong et al., 2021] . In their continuous time interpretation [Song et al., 2021] , these models leverage the idea that a diffusion process can bridge the data distribution µ ⋆ to a base distribution π, and its time reversal can transform samples from π into synthetic data from µ ⋆ . Anderson [1982] showed that the time reversal of a diffusion process, i.e., the backward process, is itself a diffusion with explicit drift and covariance functions that are related to the score functions of the time-marginal densities of the original, forward diffusion. Consequently, the key element of these generative models is learning these score functions using techniques such as (denoising) score-matching [Hyvärinen, 2005 , Vincent, 2011] .In this work we propose a new family of generative models which use piecewise deterministic Markov processes (PDMPs) as noising processes instead of diffusions. PDMPs were introduced around forty years ago [Davis, 1984 [Davis, , 1993] ] and since then have been successfully applied in various fields, including communication networks [Dumas et al., 2002] , biology [Berg and Brown, 1972, Cloez, Bertrand et al., 2017] , risk theory [Embrechts and Schmidli, 1994] , and the reliability of complex systems [Zhang et al., 2008] . More recently, PDMPs have been intensively studied in the context of Monte Carlo algorithms [Fearnhead et al., 2018] as alternatives to Langevin diffusionbased methods and Metropolis-Hastings mechanisms. This renewed interest in PDMPs has led to the development of novel processes, such as the Zig-Zag process (ZZP) [Bierkens et al., 2019a] , the Bouncy Particle Sampler (BPS) [Bouchard-Côté et al., 2018] , and the Randomised Hamiltonian Monte Carlo (RHMC) [Bou-Rabee and Sanz-Serna, 2017] . PDMPs offer several advantages compared to Langevin-based methods, such as better scalability and reduced computational complexity in high-dimensional settings [Bierkens et al., 2019a] . In the context of generative modelling PDMPs offer several potential advantages over diffusion processes. A key strength is their ability to effectively model data distributions supported on constrained or restricted domains. By adjusting their deterministic dynamics, PDMPs can easily incorporate boundary behaviour, making them straightforward to implement in such settings [Bierkens et al., 2018 , Davis, 1993] . Similarly, PDMPs can model data on Riemannian manifolds by employing flows that respect the manifold's geometry (see, e.g., Yang et al. [2022] for a PDMP on the sphere). Moreover, PDMPs are well-suited for modelling data distributions that combine a continuous density and a positive mass on a lower dimensional manifold [Bierkens et al., 2022] .Our contributions are the following:1) Leveraging the existing literature on time reversals of Markov jump processes [Conforti and Léonard, 2022] , we characterise the time reversal of any PDMP under appropriate conditions. It turns out that this time reversal is itself a PDMP with characteristics related to the original PDMP; see Proposition 1.2) We further specify the characteristics of the time-reversal processes associated with the three aforementioned PDMPs: ZZP, BPS, and RHMC. For these processes, Proposition 2 shows the corresponding time-reversals are PDMPs with simple reversed deterministic motion and with jump rates and kernels that depend on (ratios of) conditional densities of the velocity of the forward process before and after a jump. In contrast to common diffusion models, the emphasis is on distributions of the velocity, similar to the case of the underdamped Langevin diffusion [Dockhorn et al., 2022] , which includes an additional velocity vector akin to the PDMPs we consider. Moreover, the structure of the backward jump rates and kernels closely connects to the case of continuous time jump processes on discrete state spaces [Sun et al., 2023 , Lou et al., 2024] .3) We define our piecewise deterministic generative models employing either ZZP, BPS, or RHMC as forward process, transforming data points to a noise distribution of choice, and develop methodologies to estimate the backward rates and kernels. Then, we define the corresponding backward process based on approximations of the time reversed ZZP, BPS, and RHMC obtained with the estimated rates and kernels. In Section 4 we test our models on simple toy distributions. 4) We obtain a bound for the total variation distance between the data distribution and the distribution of our generative models taking into account two sources of error: first, the approximation of the characteristics of the backward PDMP, and second, its initialisation from the limiting distribution of the forward process; see Theorem 1. 
pdmp based generative models 2.1 piecewise deterministic markov processes: 
 Informally, a PDMP [Davis, 1984 [Davis, , 1993] ] on the measurable space (R D , B(R D )) is a stochastic process that follows deterministic dynamics between random times, while at these times the process can evolve stochastically on the basis of a Markov kernel. In order to define a PDMP precisely, we need three components, called characteristics of the PDMP: a vector field Φ : R + × R D → R D , which governs the deterministic motion, a jump rate λ : R + × R D → R + , which defines the law of random event times, and finally a jump kernel Q : R + × R D × B(R D ) → [0, 1], which is applied at event times and defines the new state of the process. Let us give an informal description of the evolution of a PDMP Z t , clarifying the role of the three characteristics. Suppose at time T ∈ R + the PDMP is at state z ∈ R D , that is Z T = z. The deterministic motion of the PDMP is described by the ODE dZ T+s = Φ(T + s, Z T+s )ds for s ⩾ 0, with initial condition Z T = z. We introduce the differential flow φ : (t, s, z) → φ t,t+s (z), which solves the ODE in the sense that dφ t,t+s (z) = Φ(T + s, φ t,t+s (z))ds for s ⩾ 0. The process evolves deterministically according to φ until the next event time T + τ , where τ is a random variable with law P(τ > s|Z T = z) = exp(-s 0 λ(φ T,T+u (z))du), i.e. the exponential distribution with non-homogeneous rate s → λ(φ T,T+s (z)). We can, at least in principle, simulate τ by solving τ = inf t > 0 :t 0 λ(T + u, φ T,T+u (z))du ⩾ E (1)where E ∼ Exp(1). The process is then defined on [T, T + τ ) by Z T+t = φ T,T+t (Z T ) for t ∈ [0, τ ). At time T + τ the process jumps to a new state that is drawn from the Markov kernel Q, hence we set Z T+τ ∼ Q(T + τ, φ T,T+τ (z), • ). A realisation of the path of a PDMP for a given time horizon can then be obtained following this procedure (see also Algorithm 1 in Appendix C.1 for a pseudo-code). The formal construction of a PDMP can be found in Appendix A.1.Typically a PDMP has several types of jumps, which can be represented by a family of jump rates and kernels (λ i , Q i ) i∈{1,...,ℓ} . A PDMP of such type can be obtained with the construction we have described by settingEQUATIONAn alternative, equivalent construction of a PDMP with λ, Q satisfying (2) is given in Appendix A.2. Finally, we say a PDMP is homogeneous (as opposed to the non-homogeneous case we have described) when the characteristics do not depend on time, that is Φ :EQUATIONQ : R D × B(R D ) → [0, 1].In all this work, we suppose that the PDMPs that we consider are non-explosive in the sense of Davis [1993] , that is they are such that the time of the n-th random event goes to +∞ as n → +∞, almost surely (see Durmus et al. [2021] for conditions ensuring this).We now introduce the three PDMPs we consider throughout the paper. All these PDMPs are timehomogeneous and live on a state space of the formE = R d × V, for V ⊂ R d , assuming V 0 ∈ V.Then, Z t can be decomposed as Z t = (X t , V t ), where X t ∈ R d is the component of interest and has the interpretation of the position of a particle, whereas V t ∈ V is an auxiliary vector playing the role of the particle's velocity. In the sequel, if there is no risk of confusion, we take the convention that any z ∈ R d × V, and we write z = (x, v) for x ∈ R d and v ∈ V. All the PDMPs below have a stationary distribution of the form π(dx) ⊗ ν(dv), where π has density proportional to x → e -ψ(x) , for ψ : R d → R a continuously differential potential, and ν is a simple distribution on V for the velocity vector. In our experiments we take π to be the standard normal distribution, while ν is the standard normal when V = R d or the uniform distribution when V is a compact set. Figure 1 shows sample paths for the position vector of the three PDMPs we introduce below.Figure 1 : Trace plots for ZZP (left), BPS (centre), RHMC (right). In all cases λ r = 1 and T f = 10.The Zig-Zag process The Zig-Zag process (ZZP) [Bierkens et al., 2019a ] is a PDMP with state spaceE Z = R d × {-1, 1} d .The deterministic motion is determined by the homogeneous vector field Φ Z (x, v) = (v, 0) T , i.e. the particle moves with constant velocity v. For i ∈ {1, . . . , d} we define the jump ratesλ Z i (x, v) := (v i ∂ i ψ(x)) + + λ r, where (a) + = max(0, a), ∂ i denotes the i-th partial derivative, and λ r ⩾ 0 is a user chosen refreshment rate. The corresponding (deterministic) jump kernels are given by Q , dw) , where δ z denotes the Dirac measure at z ∈ E. Here, R Z i is the operator that reverses the sign of the i-th component of the vector to which it is applied, i.e.Z i ((x, v), (dy, dw)) = δ (x,R Z i v) (dyR Z i v = (v 1 . . . , v i-1 , -v i , v i+1 , . . . , v d ).The ZZP falls within our definition of PDMP taking λ, Q as in (2). As shown in Bierkens et al. [2019a] , the ZZP has invariant distribution π ⊗ ν, where ν is the uniform distribution over {±1} d . Moreover, Bierkens et al. [2019b] shows that for any λ r ⩾ 0 the law of the ZZP converges exponentially fast to its invariant distribution e.g. when π is a standard normal distribution.The Bouncy Particle sampler The Bouncy Particle sampler (BPS) [Bouchard-Côté et al., 2018 ] is a PDMP with state space isE B = R d ×V B , where V B = R d or V B = S d-1 := {v ∈ R d : ∥v∥ = 1}.The deterministic motion is governed as ZZP by the homogeneous vector field defined for z = (x, v) ∈ E by Φ B (x, v) = (v, 0) T . Now we introduce two jump rates which correspond to two types of random events: reflections and refreshments. Reflections enforce that µ(x, v) = π(x)ν(v) is the invariant density of the process, where π(dx) ∝ exp(-ψ(x))Leb(dx) is a given distribution and ν is either a standard normal distribution when V B = R d or the uniform distribution on S d-1 when V B = S d-1 . Reflections are associated to the homogeneous jump rate(x, v) → λ B 1 (x, v) = ⟨v, ∇ψ(x)⟩ + , while refreshments are associated to (x, v) → λ B 2 (x, v) = λ r for λ r > 0. The corresponding jump kernels are Q B 1 ((x, v), (dy, dw)) = δ (x,R B x v) (dy, dw) , Q B 2 ((x, v), (dy, dw)) = δ x (dy)ν(dw), where R B x v = v -2( ⟨v,∇ψ(x)⟩ /|∇ψ(x)| 2 )∇ψ(x) . The operator R Bx reflects the velocity v off the hyperplane that is tangent to the contour line of ψ passing though point x. The norm of the velocity is unchanged by the application of R B , and this gives the interpretation that R B is an elastic collision of the particle off such hyperplane. As observed in Bouchard-Côté et al. [2018] , BPS requires a strictly positive λ r to avoid being reducible, that is to make sure the process can reach any area of the state space. Exponential convergence of the BPS to its invariant distribution was shown by Deligiannidis et al. [2019] , Durmus et al. [2020] .Randomised Hamiltonian Monte Carlo Randomised Hamiltonian Monte Carlo (RHMC) [Bou-Rabee and Sanz-Serna, 2017] refers to the PDMP with state space E H = R d × R d which is characterised by Hamiltonian deterministic flow and refreshments of the velocity vector from the standard normal distribution. The flow is governed by the homogeneous vector field defined by(x, v) → Φ H (x, v) = (v, -∇ψ(x)) T, where ψ is the potential of π. The jump rate coincides with the refreshment part of BPS, i.e., it is the constant function λ H : (x, v) → λ r > 0 and jump kernel Q H ((x, v), (dy, dw)) = δ x (dy)ν(dw). When the stationary distribution π is a standard Gaussian, the deterministic dynamics (x t , v t ) t⩾0 satisfy dx t = v t dt, dv t = -x t dt, which for t ⩾ 0 has solution x t = x 0 cos(t) + v 0 sin(t) and v t = -x 0 sin(t) + v 0 cos(t), where (x 0 , v 0 ) is the initial condition. It is well known that Hamiltonian dynamics preserve the density µ(x, v) = π(x)ν(v) [Neal, 2010] , where ν is the standard normal distribution, while velocity refreshments are necessary to ensure the process is irreducible. Exponential convergence of the law of this PDMP to µ was shown in Bou-Rabee and Sanz-Serna [2017] .Remark 1 (Noise schedule) Similarly to diffusion models we can introduce a noise schedule β(t) that regulates the amount of randomness injected at time t. This can be achieved using the time change of a given PDMP with characteristics (Φ, λ, Q) as forward process, resulting in the PDMP with characteristics (Φ β , λ β , Q) for Φ β (t, z) = β(t)Φ(t, z) and λ β (t, z) = β(t)λ(t, z).In this section we characterise the time reversal of a PDMP. This key result, stated in Proposition 1, is essential to be able to use PDMPs for generative modelling. The time reversal of a PDMP (Z t ) t∈[0,T f ] with initial distribution µ 0 is the process that at time t ∈ [0, T f ] has distribution µ 0 P T f -t , where µ 0 P s denotes the law of Z s . It follows that the law of the time reversal at time T f is µ 0 , which is the key observation in the context of generative modelling. Characterisations of the law of time reversed Markov processes with jumps were obtained in Conforti and Léonard [2022] and in the following statement we adapt their Theorem 5.7 to our setting, showing that the time reversal of a PDMP with characteristics (Φ, λ, Q) is a PDMP with reversed deterministic motion and jump rates and kernels satisfying (3).Proposition 1 Consider a non-explosive PDMP (Z t ) t⩾0 with characteristics (Φ, λ, Q) and initial distribution µ 0 on R D . In addition, let T f be a time horizon. Suppose that Φ is locally bounded, (t, z) → λ(t, z) is continuous in both its variables, andT f 0 E[λ(t, Z t )]dt < ∞.Assume the technical conditions H3, H4, postponed to the appendix. Then, the corresponding time reversal process is a PDMP with characteristics (← - Φ , ← - λ , ← - Q ), where ← - Φ (t, z) = -Φ(T f -t, z) and ←λ , ← -Q are the unique solutions to the following balance equation: for almost all t ∈ [0, T f ],EQUATION) where µ 0 P t stands for the distribution of Z t starting from µ 0 .The proof is postponed to Appendix A.4. The condition H3 is standard in the literature on PDMPs [Davis, 1993] and is verified for ZZP, BPS, and RHMC. H4 is a technical assumption on the do-main of the generator of the forward PDMP and has been shown to hold e.g. for the ZZP. In the next proposition we derive expressions for the backward jump rate and kernel satisfying (3) corresponding to a forward PDMP with characteristics with the same structure as those of ZZP, BPS, and RHMC. We state the result assuming the PDMP has only one jump type, but the generalisation to the case of ℓ > 1 jump mechanisms of the form (2) can be immediately obtained applying Proposition 2 to each pair (λ i , Q i ) for i ∈ {1, . . . , ℓ}. We refer to Appendix A.6 for the details.Proposition 2 Consider a non-explosive PDMP (X t , V t ) t⩾0 with characteristics (Φ, λ, Q) and initial distribution µ X 0 ⊗µ V 0 on R 2d . In addition, let T f be a time horizon. Suppose that Φ and λ satisfy the same conditions as Proposition 1, in particular the technical conditions H3, H4 postponed to the appendix. Suppose in addition that for any t ∈ (0, T f ], the conditional distribution of V t given X t has a transition density (x, v) → p t (v|x) with respect to some reference measure µ V ref on R d .(1) (Deterministic jumps). Suppose Q((y, w), (dx, dv)) = δ y (dx)δ Ryw (dv) where for anyy ∈ R d , R y : R d → R d is an involution which preserves µ V ref , i.e., R -1 y = R y and µ V ref (dR y w) = µ V ref (dw).Then for almost all t ∈ [0, T f ] and any (y, w) ∈ R 2d such that p T f -t (w|y) > 0 it holds that← - λ (t, (y, w)) = p T f -t (R y w|y) p T f -t (w|y) λ(T f -t, (y, R y w)) , ← - Q ((y, w), (dx, dv)) = δ y (dx)δ Ryw (dv) .(2) (Refreshments). Suppose Q((y, w), (dx, dv)) = δ y (dx)ν(dv|y), where ν is a transition kernel on R d × B(R d ), and λ(t, (y, w)) = λ(t, y). Suppose also for any y ∈ R d , ν(•|y) is absolutely continuous with respect to µ V ref .Then for almost all t ∈ [0, T f ] and any (y, w) ∈ R 2d such that p T f -t (w|y) > 0 it holds that← - λ (t, (y, w)) = ( dν /dµ V ref )(w|y) p T f -t (w|y) λ(T f -t, y), ← - Q (t, (y, w), (dx, dv)) = δ y (dx)p T f -t (v|x)µ V ref (dv).The proof is postponed to Appendix A.5. We remark that we consider that µ V 0 is a distribution on R d also when µ V 0 (V) = 1 for V ⊂ R d , in which case the reference measure can simply be chosen such that µ V ref (V) = 1. Applying Proposition 2 we are able to derive explicit expressions for the characteristics of the time reversals of ZZP, RHMC, and BPS. The rigorous statements and their proofs can be found in Appendix A.7. For ZZP and BPS we assume the following condition on π, the limiting distribution for the position vector of the forward process.H1 Recall π(x) ∝ e -ψ(x) . It holds that ψ ∈ C 2 (R d ) and sup x∈R d ∥∇ 2 ψ(x)∥ < +∞.This assumption is satisfied e.g. by any multivariate normal distribution. For BPS and RHMC we suppose that for any t ∈ (0, T f ], the conditional distribution of V t given X t has a transition density (x, v) → p t (v|x) with respect to the Lebesgue measure. Moreover, for all samplers we assume H4. Time reversal of ZZP In order to apply Proposition 2 we additionally assume that |∂ i ψ(x)|dµ ⋆ (x) < ∞ for all i = 1, . . . , d. We find that the deterministic motion is defined by ← -Φ Z (y, w) = (-w, 0) T for any (y, w) ∈ R 2d , while the backward rates and kernels are for i = 1, . . . , d and for all (y, w) ∈ R 2d such that p T f -t (w|y) > 0,← - λ Z i (t, (y, w)) = p T f -t (R Z i w|y) p T f -t (w|y) λ Z i (y, R Z i w) , ← - Q Z i ((y, w), (dx, v)) = δ (y,R Z i w) (dx, v) . (4)Time reversal of BPS Whereas in Appendix A.7 we consider the case where the velocity of BPS is initialised on S d-1 , we can formally apply Proposition 2 to the case of ν is the standard ddimensional Gaussian distribution assuming that |∇ψ(x)|dµ ⋆ (x) < ∞. The drift of the backward BPS is clearly the same as for the backward ZZP, while jump rates and kernels are for all t ∈ [0, T f ] and (y, w) ∈ R 2d such that p T f -t (w|y) > 0← - λ B 1 (t, (y, w)) = p T f -t (R B y w|y) p T f -t (w|y) λ B 1 (y, R B y w), ← - Q B 1 ((y, w), (dx, dv)) = δ (y,R B y w) (dx, dv) , ← - λ B 2 (t, (y, w)) = λ r ν(w) p T f -t (w|y) , ← - Q B 2 (t, (y, w), (dx, dv)) = p T f -t (v|y)δ y (dx)dv . (5)Time reversal of RHMC. The deterministic motion of the backward RHMC follows the system of ODEs ← -Φ H (x, v) = (-v, ∇ψ(x)) T , which, when the limiting distribution π is Gaussian, has solution x t = x 0 cos(t) -v 0 sin(t) and v t = x 0 sin(t) + v 0 cos(t). The backward refreshment rate and kernel coincide with those of BPS as given in (5).Remark 2 (Variance exploding PDMPs) Similarly to the case of diffusion models [Song et al., 2021] , we can define variance exploding PDMPs choosing ψ(x) = 0 for all x ∈ R d , that is when π(dx) is the Lebesgue measure. In this case, the deterministic motion of RHMC coincides with ZZP and BPS, and all three processes have only velocity refreshment events.In Section 2.2 we showed that the backward jump rates and kernels of ZZP, BPS, and RHMC, involve the conditional densities of the velocity vector of the forward process given its position vector at all times t ∈ [0, T f ]. These conditional densities are unavailable in analytic form, hence in this section we develop methods to learn the jump rates and kernels of our time reversed PDMPs. In Appendix D we give the pseudo codes and more detailed descriptions of the training procedure for our models, together with a comparison with diffusion models.Approximating the jump rates of the backward ZZP via ratio matching In the case of ZZP, we need to approximate for any i ∈ {1, . . . , d}, the rates in (4). Since the terms λ Z i (x, R Z i v) are known, it is sufficient to estimate the density ratios r Z i (x, v, t) := pt(R Z i v|x) /pt(v|x) for all states (x, v) such that p t (v|x) > 0. To this end, we introduce a class of functions {s θ :R d × {-1, 1} d × [0, T f ] → R d + :θ ∈ Θ} for some parameter set Θ ⊂ R d θ and aim to find a parameter θ ⋆ ∈ Θ such that for any i ∈ {1, . . . , d}, the i-th component of s θ⋆ , denoted by s θ⋆ i (•), is an approximation of r Z i . We then approximate the backward ZZP by using the rates λZi (t, (x, v)) = s θ⋆ i (x, v, T f -t) λ Z i (x, R Z i v).To address the problem of fitting θ, we consider different loss functions inspired by the ratio matching (RM) problem considered in Hyvärinen [2007] .From a discrete probability density p ± on {-1, 1} d , RM consists in learning the d ratios v → p±(Riv) /p±(v) for i ∈ {1, . . . , d}. This problem was motivated in Hyvärinen [2007] as a means to estimate p ± without requiring its normalising constant, similarly to score matching applied to estimate continuous probability densities [Hyvärinen, 2005] . In our context we are interested only in the ratios, hence as opposed to Hyvärinen [2007] we do not model the conditional distributions (x, v) → p t (v|x), but directly the ratios r Z i . Adapting the ideas of Hyvärinen [2007] to our context, we introduce the function G : r → (1 + r) -1 and define the Explicit Ratio Matching objective functionEQUATION)where ω : [0, T f ] → R * + is a probability density, and(X t , V t ) t⩾0 is a ZZP initialised from µ ⋆ ⊗ Unif({-1, 1} d ).This objective function considers simultaneously the square error in the estimation of both(x, v, t) → r i (x, v, t) and (x, v, t) → r i (x, R Z i v, t), where the function G improves numerical stability, particularly when one of the two ratios is very small. Clearly ℓ E (θ) = 0 if and only if s θ i (x, v, t) = r i (x, v, t) for almost all x, v, t and all i. Moreover, the choice of G allows us to optimise without knowledge of the true ratios, as shown in the following result.EQUATIONwhere(X t , V t ) t∈R+ is a ZZP starting from µ ⋆ ⊗ Unif({-1, 1} d ).Therefore we aim to solve the minimisation problem associated with ℓ I , which has for empirical counterpartθ → 1 N N n=1 d i=1 G 2 (s θ i (X n τ n , V n τ n , τ n )) + G 2 (s θ i (X n τ n , R Z i V n τ n , τ n )) -2G(s θ i (X n τ n , V n τ n , τ n ))where {τ n } N n=1 are i.i.d. samples from ω, independent of {(X n t , V n t ) t⩾0 } N n=1 , which are N i.i.d. realisations of the ZZP respectively starting at the n-th training data point with velocity V n 0 , where{V n 0 } N n=1 are i.i.d. observations of Unif({-1, 1} d ).Notice that the loss above has computational cost increasing linearly in d because d + 1 evaluations of the model are needed for each datum. This can be improved considering an estimate for the ratio which does not take as input the whole velocity vector (see Appendix D.1 for the details). This variation has computational cost that is constant in the dimension, but might have lower accuracy.Approximating the characteristics of BPS and RHMC For BPS and RHMC, Proposition 2 shows that if we aim to sample from the backward process, we have to estimate both ratios of the conditional density of the velocity of the forward PDMP given its position at any time t ∈ [0, T f ], and also to be able to sample from such densities as prescribed by the backward jump kernel (5). In order to address both requirements, we introduce a parametric family of conditional probability distributions {p θ : θ ∈ Θ} of the form (x, v, t) → p θ (v|x, t), where Θ ⊂ R d θ , which we model with the framework of normalising flows (NFs) [Papamakarios et al., 2021] . The advantage of NFs lays in their feature that, once the network is learned, it is possible both to obtain an estimate of the density at a given state and time, and also to generate samples which are approximately from (x, v, t) → p t (v|x). However, training conditional NFs can be challenging in high dimensions.Focusing on BPS, we now illustrate how we can use NFs to learn the backward jump rates and kernels. We aim to find a parameter θ B ⋆ such that p θ B ⋆ (v|x, t) is a good approximation of p t (v|x), the conditional density of the forward BPS with respect to the Lebesgue measure. We choose to optimise θ following the maximum likelihood approach, which gives the theoretical lossEQUATIONwhere ω : [0, T f ] → R * + is a probability density, and (X t , V t ) t⩾0 is a a BPS initialised from µ ⋆ ⊗ ν, with ν denoting the density of the d-dimensional standard normal distribution. The optimal parameter θ B ⋆ can then be found minimising the empirical counterpart of ℓ ML (θ):EQUATIONwhere {τ n } N n=1 are i.i.d. samples from ω, independent of {(X n t , V n t ) t⩾0 } N n=1 , which are N i.i.d. realisations of the ZZP respectively starting at the n-th training data point with velocity V n 0 , where {V n 0 } N n=1 are i.i.d. observations from the multivariate standard normal distribution. Once we have obtained the optimal parameter θ B ⋆ , we can define our approximation of the backward refreshment mechanism of BPS taking the rate λB2 (t, (x, v)) = λ r × ν(v) /p θ B ⋆ (v|x,T f -t) and the kernel QB 2 (t, (y, w), (dx, dv)) = p θ B ⋆ (v|y, T f -t)δ y (dx)dv.Similarly, we estimate the backward reflection ratio of BPS as λB1 (t, (x, v)) = λ B 1 (x, R B x v) × p θ B ⋆ (R B x v|x,T f -t) /p θ B ⋆ (v|x,T f -t).We now discuss how we can simulate the backward PDMP with exact backward flow map (t, x, v) → φ -t (x, v) and jump characteristics λ and Q that are approximations of the jump rates and kernels of the time reversed PDMPs obtained as discussed in Section 2.3. We recall that the backward rates have the general form λ(t,(x, v)) = s θ (x, v, T f -t)λ(x, Rv), where s θ is an estimate of a density ratio and R is a suitable involution. In principle such a PDMP can be simulated following the procedure described in Section 2.1, but the generation of the random jump times via (1) requires the integration of λ(t, φ -t (x, v)) with respect to t. This cannot be achieved since λ is defined through a neural network. A standard approach in the literature (see e.g. Bertazzi et al. [2022 Bertazzi et al. [ , 2023]] ) is to discretise time and (informally) approximate the integral in (1) with a finite sum.Here we focus on approximations based on splitting schemes discussed in Bertazzi et al. [2023] , adapting their ideas to the non-homogeneous case. Such splitting schemes approximate a PDMP with a Markov chain defined on the time grid {t n } n∈{0,...,N } , with t 0 = 0 and t N = T f . The key idea is that the deterministic motion and the jump part of the PDMP are simulated separately in a suitable order, obtaining second order accuracy under suitable conditions (see Theorem 2.6 in Bertazzi et al. [2023] ). Now, we give an informal description of the splitting scheme that we use for RHMC, that is based on splitting DJD in Bertazzi et al. [2023] , where D stands for deterministic motion and J for jumps. We define our Markov chain based on the step sizes {δ j } j∈{1,...,N } , where δ j = t j -t j-1 . Suppose we have defined the Markov chain on {t k } k∈{0,...,n} for n < N and that the state at time t n is (x tn , v tn ). The next state is obtained following three steps. First, the particle moves according to its deterministic motion for a half-step, that is we define an intermediate state (x tn , ṽtn ) = φ -δ n+1/2 (x tn , v tn ). Second, we turn our attention to the jump part of the process. In this phase, the particle is only allowed to move through jumps and there is no deterministic motion. This means that the rate is frozen to the value λ(t n + δn+1 /2, (x tn , ṽtn )) and thus the integral in (1) can be computed trivially. The proposal for the next event time is then given byτ n+1 ∼ Exp(λ(t n + δn+1 /2, (x tn , ṽtn ))). If τ n+1 ⩽ δ n+1 , we draw w ∼ Q(t n + δn+1 /2, (x tn , ṽtn ), •)and update ṽtn = w, else we leave ṽtn unchanged. Finally we conclude with an additional half-step of deterministic motion, letting (x tn+1 , v tn+1 ) = φ -δ n+1/2 (x tn , ṽtn ). We refer to Appendix C.2 for a detailed description of the schemes used for each process together with the pseudo-codes. 
error bound in total variation distance: 
 In this section, we give a bound on the total variation distance between the data distribution µ ⋆ and the law of the synthetic data generated by a PDMP with initial distribution π ⊗ ν and approximate characteristics obtained, e.g., with the methods described in Section 2.3. We obtain our result comparing the law of such PDMP to the law of the exact time reversal obtained in Section 2.2, that is the PDMP with the analytic characteristics of Proposition 2 and with initial distribution L(X T f , V T f ), i.e. the law of the forward PDMP at time T f when initialised from µ ⋆ ⊗ ν. In Theorem 1 below we then take into account two of the three sources of error in our models: (i) the error introduced initialising the backward PDMP from the limiting distribution of the forward, (ii) the error due to the approximation of the backward rates and kernels. For simplicity we neglect the discretisation error caused by the methods discussed in Section 2.4.We shall assume the following condition, which deals with the error introduced by initialising the backward PDMP from π ⊗ ν.The forward PDMP with semigroup (P t ) t⩾0 is such that there exist γ, C > 0 for which∥π ⊗ ν -µ ⋆ ⊗ νP t ∥ TV ⩽ Ce -γt .Informally, H2 is verified for some C < ∞ when π is a multivariate standard Gaussian distribution for ZZP and BPS if the tails of µ ⋆ are at least as light as those of π, while for RHMC it is enough if µ ⋆ has finite second moments. We refer to Appendix E.1 for a more detailed discussion on this aspect. We are now ready to state our result.Theorem 1 Consider a non-explosive PDMP (X t , V t ) t⩾0 with initial distribution µ ⋆ ⊗ν, stationary distribution π⊗ν, and characteristics (Φ, λ, Q). Let T f be a time horizon. Suppose the assumptions of Proposition 1 as well as H2 hold. LetEQUATION)where The proof is postponed to Appendix E.2. The first term in ( 10) is caused by initialising the process from π ⊗ ν, while the second term represents the error introduced by the approximate jump rate λ and kernel Q. For the sake of illustration we obtain a simple upper bound to (10) in the case of ZZP (for the details see Appendix E.3). We assume the conditions of Theorem 1 are satisfied and also that the expected error of the learned rates λZ i is bounded by a constant, in the sense thatEQUATIONE[|r Z i (X t , V t , T f -t) -s θ i (X t , V t , T f -t)|λ Z i (X t , R Z i V t )]⩽ M for all t ∈ [0, T f ] and i ∈ {1, . . . , d}. The latter condition is similar to the standard assumption asked on the approximation of the score in diffusion models [Chen et al., 2023] . Under these conditions we obtain the boundEQUATION)4 Numerical simulationsIn this section, we test our piecewise deterministic generative models on simple synthetic datasets.Design We compare the generative models based on ZZP, BPS, and RHMC with the improved denoising diffusion probabilistic model (i-DDPM) given in Nichol and Dhariwal [2021] . For all of our models, we choose the standard normal distribution as target distribution for the position vector, as well as for the velocity vector in the cases of BPS and RHMC. The accuracy of trained generative models is evaluated by the kernel maximum mean discrepancy (MMD). We refer to Appendix F for a detailed description of the parameters and networks choices.In Table 1 we report the MMD score for five, 2-dimensional toy distributions. We observe that the PDMP based generative models perform well compared to i-DDPM in all of these five datasets. In particular, ZZP and i-DDPM are implemented with the same neural network architecture, hence ZZP appears to compare favourably to i-DDPM with the same model expressivity.The results of Table 1 are supported by the plots of generated data shown in Figure 2 , illustrating how ZZP and BPS are able to generate more detailed edges compared to i-DDPM.In Figure 4 , we compare the output of RHMC and i-DDPM for a very small number of reverse steps. We observe how in this setting the data generated by RHMC are noticeably closer to the true data distribution compared to i-DDPM. This phenomenon is observed also for BPS as shown in Table 2 , and is intuitively caused by the refreshment kernel, which is able to generate velocities that correct wrong positions. Respecting this intuition, ZZP does not perform as well as BPS and RHMC for a small number of reverse steps since its velocities are constrained to {-1, 1}. Nonetheless, ZZP generates the most accurate results in our experiments given a large enough number of reverse steps.Table 2 and Figure 3 show that PDMP-based models require a smaller computational time to generate samples of a given quality compared to i-DDPM. This is the case because PDMP models require considerably less backward steps than i-DDPM, although each step is more expensive (see Table 2 ). Additional results can be found in Appendix F, including promising results applying ZZP to the MNIST dataset.  
discussion and conclusions: 
 We have introduced new generative models based on piecewise deterministic Markov processes, developing a theoretically sound framework with specific focus on three PDMPs from the sampling literature. While this work lays the foundations of this class of methods, it also opens several directions worth investigating in the future.Similarly to other generative models, our PDMP based algorithms are sensitive to the choice of the network architecture that is used to approximate the backward characteristics. Therefore, it is crucial to investigate which architectures are most suited for our algorithms in order to achieve state of the art performance in real world scenarios. For instance, in the case of BPS and RHMC it could be beneficial to separate the estimation of the density ratios and the generation of draws of the velocity conditioned on the position and time. For the case of ZZP, efficient techniques to learn the network in a high dimensional setting need to be investigated, while network architectures that resemble those used to approximate the score function appear to adapt well to the case of density ratios. Moreover, there are several alternative PDMPs that could be used as generative models and that we did not consider in detail in this paper, as for instance variance exploding alternatives.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: we provide all the necessary codes to reproduce our experiments. Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: All the design choices are specified in the paper or in the supplementary material.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We report the standard deviations to our tables, or report when the results are based on a single run.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We provide information on the type of compute workers that were used.Guidelines:• The answer NA means that the paper does not include experiments. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper poses no such risks. 
 
 Appendix: The Appendix is organised as follows. Appendix A includes the details and the proofs regarding Section 2.1 and Section 2.2. Appendix B contains details and proofs regarding the framework of density ratio matching. Appendix C gives details and pseudo-codes for the exact simulation of our forward processes (see Appendix C.1), and also for the splitting schemes that are used to approximate the backward processes (see Appendix C.2). Appendix D details the algorithms used to train our generative models and outlines the computational differences with the popular framework of denoising diffusion models. Appendix E contains the proof for Theorem 1 and some related details Finally, Appendix F contains the details on the numerical simulations, as well as additional results.A.1 Construction of a PDMPHere we describe the formal construction of a PDMP with the characteristics (Φ, λ, Q). To this end, consider the differential flow φ : (t, s, z) → φ t,t+s (z), which solves the ODE, dz t+s = Φ(t + s, z t+s )ds for s ⩾ 0, i.e. z t+s = φ t,t+s (z t ). We define by recursion on n ∈ N the process on (Z t ) t∈[0,Tn] on [0, T n ] and the increasing sequence of jump times (T n ) n∈N starting from an initial state Z 0 and setting T 0 = 0. Assume that (T i ) i∈{0,...,n} and (Z t ) t∈[0,Tn] are defined for some n ∈ N. We now define (Z t ) t∈ [Tn, Tn+1] . First, we defineEQUATIONwhere E n+1 ∼ Exp(1), and set the n + 1-th jump timeT n+1 = T n + τ n+1 . The process is then defined on [T n , T n+1 ) by Z Tn+t = φ Tn,Tn+t (Z Tn ) for t ∈ [0, τ n+1 ). Finally, we set Z Tn+1 ∼ Q(T n+1 , φ Tn,Tn+τn+1 (Z Tn ), •).The process (Z t ) t⩾0 is a Markov process by [Jacobsen, 2005, Theorem 7.3.1] .In this section we describe the formal construction of a non-homogeneous PDMP with the characteristics (Φ, λ, Q) where λ, Q are of the formEQUATIONRecall the differential flow φ : (t, s, z) → φ t,t+s (z), which solves the ODE, dz t+s = Φ(t + s, z t+s )ds for s ⩾ 0, i.e. z t+s = φ t,t+s (z t ). Similarly to the case of one type of jump only, we start the PDMP from an initial state Z 0 , assume it is defined as (Z t ) t∈[0,Tn] on [0, T n ] for some n ∈ N, and we now define define (Z t ) t∈ [Tn, Tn+1] . First, we define the proposals (τ i n+1 ) i∈{1,...,ℓ} for next event time asτ i n+1 = inf t > 0 : t 0 λ i (T n + u, φ Tn,Tn+u (Z Tn ))du ⩾ E i n+1where E i n+1 ∼ Exp(1) for i ∈ {1, . . . , ℓ}. Then define i * = arg min i∈{1,...,ℓ} τ i n+1 and set the next jump time toT n+1 = T n + τ i * n+1 . The process is then defined on [T n , T n+1 ) by Z Tn+t = φ Tn,Tn+t (Z Tn ) for t ∈ [0, τ n+1 ). Finally, we set Z Tn+1 ∼ Q i * (T n+1 , φ Tn,Tn+τn+1 (Z Tn ), •).In order to obtain the generator of a PDMP, Theorem 26.14 of Davis [1993] requires "standard conditions" on the characteristics (see conditions (24.8) in Davis [1993] ). We state these conditions for a non-homogeneous PDMP in the next assumption.The non-homogeneous characteristics (Φ, λ, Q) satisfy the following conditions:1. Φ is locally Lipschitz and the associated flow map φ has infinite explosion time; 2. λ is such that u → λ(φ t,t+u (x)) is integrable on [0, ε(x, t)) for some ε(x, t) > 0 and all (t, x) ∈ R + × E.3. Q is measurable and such that Q(t, x, {x}) = 0 for all (t, x) ∈ R + × E.4. Let (T n ) n∈{0,1,... } be the random sequence of event times of the PDMP and defineN t = ∞ k=0 1 t⩾T k . It holds that E x [N t ] < ∞ for all (t, x) ∈ R + × E.Notably, the PDMP is required to be non-explosive in the sense that the expected number of random events after any time t starting the PDMP from any state should be finite. These conditions are verified for all the three PDMPs we consider as forward processes. Assuming H3 we can apply Theorem 26.14 in Davis [1993] to the homogeneous PDMP obtained including the time variable, which gives that the extended generator of the non-homogeneous PDMP with characteristics (Φ, λ, Q) is given byEQUATIONfor all functions f ∈ dom(L t ), that is the space of measurable functions such thatM f t = f (Z t ) -f (Z 0 ) - t 0 L s f (Z s )dsis a local martingale. We also introduce the Carré du champ Γ t (f, g):= L t (f g) -f L t g -gL t f, with domain dom(Γ t ) := {f, g : f, g, f g ∈ dom(L t )} which in the case of a PDMP with generator (15) takes the formΓ t (f, g)(z) = λ(t, z) R d (f (y) -f (z))(g(y) -g(z))Q(t, z, dy) .In order to prove Proposition 1 we apply Conforti and Léonard [2022, Theorem 5.7 ] and hence in this section we verify the required assumptions. Before starting, we state the following technical condition which we omitted in Proposition 1 and is assumed in Conforti and Léonard [2022, Theorem 5.7] .H4 It holds C 2 c (R d ) ⊂ dom(L t )for any t ∈ R + . We now turn to verifying the remaining assumptions in Conforti and Léonard [2022, Theorem 5.7 ]. The "General Hypotheses" of Conforti and Léonard [2022] are satisfied since we assume the vector field Φ is locally bounded, the switching rate (t, z) → λ(t, z) is a continuous function, and the jump kernel Q is such that Q(t, x, •) is a probability distribution. In particular these assumptions imply thatsup t∈[0,T ],|z|⩽ρ R d (1 ∧ |z -y| 2 )λ(t, z)Q(t, z, dy) ⩽ sup t∈[0,T ],|z|⩽ρ λ(t, z) < ∞ for all ρ ⩾ 0.Then, Conforti and Léonard [2022, Theorem 5.7 ] requires a further integrability condition, which is satisfied when[0,T ]×R d ×R d (1 ∧ |z -y| 2 )µ 0 P t (dz)λ(t, z)Q(t, z, dy) < ∞.It is then sufficient to have thatEQUATIONFinally, Conforti and Léonard [2022, Theorem 5.7] requires some technical assumptions which we now discuss. Introduce the class of functions that are twice continuously differentiable and compactly supported, denoted by C 2 c (R d ), and for f ∈ C 2 c (R d ) consider the two following conditions:EQUATIONEQUATIONWe define F := {f ∈ C 2 c (E) : (17), (18) hold }. We need to verify that F ≡ C 2 c (E). Let us start by considering (17): we findT 0 µ 0 P t (dz)|L t f (z)|dt ⩽ T 0 µ 0 P t (dz) |⟨Φ(t, z), ∇f (z)⟩| + λ(t, z) |u(y) -u(z)|Q(t, z, dy) dt. Since f ∈ C 2 c (E)we have that |⟨Φ(t, z), ∇f (z)⟩| is compactly supported and hence integrable, while the second term is finite assumingT 0 E[λ(t, Z t )]dt < ∞.Under the latter assumption, (18) can be easily verified.Let us denote the initial condition of the forward PDMP by µ 0 = µ X 0 ⊗ µ V 0 . First of all, notice that, for a PDMP with position-velocity decomposition and homogeneous jump kernel, the flux equation (3) becomesµ 0 P t(dy, dw) ← - λ (t, (y, w)) ← - Q (t, (y, w), (dx, dv)) = µ 0 P t(dx, dv)λ( t, (x, v))Q((x, v), (dy, dw))where t = T f -t. Moreover, since the jump kernel leaves the position vector unchanged we obtain that this is equivalent toµ 0 P t(dw|y) ← - λ (t, (y, w)) ← - Q (t, (y, w), (dx, dv)) = µ 0 P t(dv|y)λ( t, (y, v))Q((x, v), (dy, dw)),where µ 0 P t (dw|y) is the conditional law of the velocity vector given the position vector at time t with initial distribution µ 0 .Suppose first that Q((y, w), (dx, dv)) = δ y (dx)δ Ryw (dv) for an involution R y . Then we find µ 0 P t(dw|y) ←λ (t, (y, w)) ← -Q (t, (y, w), (dx, dv)) = µ 0 P t(dR y w|y)λ( t, (y, R y w))δ y (dx)δ Ryw (dv)where we used that δ Ryw (dv) = δ Ryv (dw) since R y is an involution. Under our assumptions we have µ 0 P t(dR y w|y) = p t(R y w|y)µ λ( t, y), ← -Q (t, (y, w), (dx, dv)) = µ 0 P t(dv|y)δ y (dx).Proposition 2 considers PDMPs with one type of jump, while here we discuss the case of characteristics of the form ( 14), which is e.g. the case of ZZP and BPS. In this setting we can assume the backward jump rate and kernel have a similar structure, that is← - λ (t, z) = ℓ i=1 ← - λ i (t, z) , ← - Q (t, z, dz ′ ) = ℓ i=1 ← - λ i (t, z) ← - λ (t, z) ← - Q i (t, z, dz ′ ) ,in which case the balance condition (3) can be rewritten asµ 0 P T f -t (dy) ℓ i=1 ← - λ i (t, y) ← - Q i (t, y, dz) = µ 0 P T f -t (dz) ℓ i=1 λ i (T f -t, z)Q i (T f -t, z, dy) .It is then enough thatµ 0 P T f -t (dy) ← - λ i (t, y) ← - Q i (t, y, dz) = µ 0 P T f -t (dz)λ i (T f -t, z)Q i (T f -t, z, dy)holds for all i ∈ {1, . . . , ℓ}. It follows that it is sufficient to apply Proposition 2 to each pair(λ i , Q i ) to obtain ( ← - λ i , ← - Q i ) such that (3) holds.In this section we give rigorous statements regarding time reversals of ZZP, BPS, and RHMC. For all samplers we rely on Proposition 2 and hence we focus on verifying its assumptions. In the cases of ZZP and RHMC we assume the technical condition H4 since proving it rigorously is out of the scope of the present paper. We remark that this can be proved with techniques as in Durmus et al.[2021], which show H4 in the case of BPS.Proposition 4 (Time reversal of ZZP) Consider a ZZP (X t , V t ) t∈[0,T f ] with initial distribution µ ⋆ ⊗ ν, where ν = Unif({±1} d ) and invariant distribution π ⊗ ν, where π has potential ψ satisfying H1. Assume that H4 holds and that µ ⋆ (dx)|∂ i ψ(x)| < ∞ for all i = 1, . . . , d. Then the time reversal of the ZZP has vector field ← -Φ Z (x, v) = (-v, 0) T and jump rates and kernels are given for all (y, w) ∈ R 2d such that P T f -t (w|y) > 0 by← - λ Z i (t, (y, w)) = p T f -t (R Z i w|y) p T f -t (w|y) λ Z i (y, R Z i w), ← - Q Z i ((y, w), (dx, v)) = δ (y,R Z i w) (dx, v) for i = 1, . . . , d.Proof We verify the conditions of Proposition 2 corresponding to deterministic transitions and rely on Appendix A.6 to apply the proposition to each pair (λ Z i , Q Z i ). First notice the vector field Φ(x, v) = (v, 0) T is clearly locally bounded and (t, x) → λ(x, v) is continuous since ψ is continuously differentiable. Moreover, the ZZP can be shown to be non-explosive applying Durmus et al. [2021, Proposition 9] . Then, we need to verify (16). First, observe that E[λ i (X t , V t )] ⩽ E[|∂ i ψ(X t )|]. Then E[|∂ i ψ(X t )|] = E ∂ i ψ(X 0 ) + 1 0 ⟨X t -X 0 , ∇∂ i ψ(X 0 + s(X t -X 0 ))⟩ds ⩽ E[|∂ i ψ(X 0 )|] + E 1 0 |⟨X t -X 0 , ∇ 2 ψ(X 0 + s(X t -X 0 ))e i ⟩|dswhere e i is the i-th vector of the canonical basis. Notice that |X t -X 0 | ⩽ td. Thus we findE[|∂ i ψ(X t )|] ⩽ E[|∂ i ψ(X 0 )|] + t √ d sup x∈R d ∥∇ 2 ψ(x)∥and thereforeT f 0 E[λ(X t , V t )]dt ⩽ T f d i=1 E|∂ i ψ(X 0 )| + T f 2 √ d sup x∈R d ∥∇ 2 ψ(x)∥ .Since E µ⋆ |∂ i ψ(X)| < ∞ and because we are assuming H1, we obtain (16). Finally, notice that P t (dv|x) is absolutely continuous with respect to the counting measure on {1, -1} d , which is clearly invariant with respect to R Z i . □Proposition 5 (Time reversal of BPS) Consider a BPS (X t , V t ) t∈[0,T f ] with initial distribution µ ⋆ ⊗ ν, where ν = Unif(S d-1 ), and invariant distribution π ⊗ ν, where π has potential ψ satisfying H1. Assume that E µ⋆ [|∇ψ(X)|] < ∞. Then there exists a density p t (w|y) := d(µ0Pt)(dw|y) /ν(dw).Moreover, the time reversal of the BPS has vector field ← -Φ B (x, v) = (-v, 0) T ,while the jump rates and kernels are given for all t, y, wEQUATION)where t = T f -t.Remark 3 Under the assumption that ν is the uniform distribution on the sphere, it is natural to take µ V ref = ν, which gives that dν /dµ V ref = 1 and hence the backward refreshment rate is as in (19). When ν is the d-dimensional Gaussian distribution, the natural choice is to let µ V ref be the Lebesgue measure and hence we obtain a rate as given in (5).Proof We verify the general conditions of Proposition 2, then focusing on the deterministic jumps and the refreshments relying on Appendix A.6. The BPS was shown to be non-explosive for any initial distribution in Durmus et al. [2021, Proposition 10] . Since λ(t, (x, v)) = λ(x, v) = ⟨v, ∇ψ(x)⟩ + , with a similar reasoning of the proof of Proposition 4 we haveE[λ(X t , V t )] = E ⟨V t , ∇ψ(X 0 )⟩ + 1 0 ⟨V t , ∇ 2 ψ(X 0 + s(X t -X 0 ))(X t -X 0 )⟩ds + .Taking advantage of |V t | = 1 we have |X t -X 0 | ⩽ t and thus we findE[λ(X t , V t )] ⩽ E |∇ψ(X 0 )| + 1 0 |∇ 2 ψ(X 0 + s(X t -X 0 ))(X t -X 0 )|ds ⩽ E [|∇ψ(X 0 )|] + t sup x∈R d ∥∇ 2 ψ(x)∥.This is sufficient to obtain ( 16) since E[|∇ψ(X 0 )|] < ∞ and we assume H1. Moreover, H4 holds by Durmus et al. [2021, Proposition 23] . Finally notice that P t (dv|x) is absolutely continuous with respect to µV ref = Unif(S d-1 ), which satisfies µ B ref (R B (x)v) = µ B ref (v) for all x, v ∈ R d × S d-1. All the required assumptions in Proposition 2 are thus satisfied. □ Proposition 6 (Time reversal of RHMC) Consider a RHMC (X t , V t ) t∈[0,T f ] with initial distribution µ ⋆ ⊗ ν, where ν is the d-dimensional standard normal distribution, and invariant distribution π ⊗ ν, where π has potential ψ ∈ C 1 (R d ). Suppose that H4 holds and that for any y ∈ R d , P t (dw|y) is absolutely continuous with respect to Lebesgue measure, with density p t (w|y). Then the time reversal of the RHMC has vector field ← -Φ H (x, v) = (-v, ∇ψ(x)) T ,while the jump rates and kernels are given for all (y, w) ∈ R 2d such that p T f -t (w|y) > 0 by← - λ H 2 (t, (y, w)) = λ r ν(w) p T f -t (w|y) , ← - Q H 2 (t, (y, w), (dx, dv)) = p T f -t (v|y)δ y (dx)dv.Proof First of all, RHMC is non-explosive by Durmus et al. [2021, Proposition 8] . Then Φ is locally bounded and ( 16) is trivially satisfied. Finally, we can take µ V ref to be the Lebesgue measure. □We now describe a general approach to approximate ratios of densities based on the minimisation of Bregman divergences [Sugiyama et al., 2011] , which as we discuss is closely connected to the loss of Hyvärinen [2007] .For a differentiable, strictly convex function f we define the Bregman divergence B f (r, s) := f (r)f (s) -f ′ (s)(r -s). Given two time-dependent probability density functions on R 2d , p, q, we wish to approximate their ratio r(x, v, t) = pt(x,v) /qt(x,v) for t ∈ [0, T f ] with a parametric functions θ : R d × R d × [0, T f ] → R + by solving the minimisation problem min θ T f 0 ω(t) E B f (r(X t , V t , t), s θ (X t , V t , t)) dt,where the expectation is with respect to the joint density q t (x, v), that is (X t , V t ) ∼ q t , while ω is a probability density function for the time variable. Well studied choices of the function f include e.g. f (r) = r log r -r, that is related to a KL divergence, or f (r) = (r -1) 2 , related to the square loss, or f (r) = r log r -(1 + r) log(1 + r), which corresponding to solving a logistic regression task. Ignoring terms that do not depend on θ we can rewrite the minimisation asmin θ T f 0 ω(t) E pt f ′ (s θ (X t , V t , t))s θ (X t , V t , t) -f (s θ (X t , V t , t)) -E qt f ′ (s θ (X t , V t , t)) dt.Notably this is independent of the true density ratio and thus it is a formulation with similar spirit to implicit score matching. Naturally, in practice the loss can be approximated empirically with a Monte Carlo average.In the next statement, we show that the loss ℓ I defined in (6), or equivalently its explicit counterpart ℓ E (see Proposition 3), can be put in the framework of Bregman divergences.Corollary 1 Recall G(r) = (1 + r) -1 and let f (r) = (r-1) 2 /2. The task min θ ℓ E (θ) is equivalent to min θ d i=1 E pt B f (G(s θ i (X t , V t , t)), G(r i (X t , V t , t))) + B f (G(s θ i (X t , R Z i V t , t)), G(r i (X t , R Z i V t , t)))Proof The result follows by straightforward computations. □The proof follows the same lines as Hyvärinen [2007, Theorem 1] . We findℓ E (θ) = C + T f 0 ω(t) d i=1 E pt G 2 (s θ i (X t , V t , t)) + G 2 (s θ i (X t , R Z i V t , t)) -2G(r i (X t , V t , t))G(s θ i (X t , V t , t)) -2G(s θ i (X t , R Z i V t , t))G(r i (X t , R Z i V t , t)) dt,Algorithm 1: Pseudo-code for the simulation of a homogeneous PDMP Input : Time horizon T, initial condition (x, v).Set t = 0, (X 0 , V 0 ) = (x, v); while t < T do draw E ∼ Exp(1); compute the next event time, that is τ ∈ R + such that τ 0 λ(φ u (X t , V t ))du = E ; if t + τ > T then set Z T = φ T-t (Z t ); else simulate Z t+τ ∼ Q(φ s (Z t ), •); end set t = t + τ ; endwhere C is a constant independent of θ. Then plugging in the expression of G we can rewrite the last term asE pt G(s θ i (X t , R Z i V t , t))G(r i (X t , R Z i V t , t)) = v∈{±1} d p t (x, v)G(s θ i (x, R Z i v, t)) p t (x, R Z i v) p t (x, v) + p t (x, R Z i v) dx = E pt G(s θ i (X t , V t , t)) p t (X t , R Z i V t ) p t (X t , V t ) + p t (X t , R Z i V t ).Therefore we findℓ E (θ) = C + T f 0 ω(t) d i=1 E pt G 2 (s θ i (X t , V t , t)) + G 2 (s θ i (X t , R Z i V t , t)) - 2G(s θ i (X t , V t , t)) p t (X t , V t ) p t (X t , V t ) + p t (X t , R Z i V t ) - 2G(s θ i (X t , V t , t)) p t (X t , R Z i V t ) p t (X t , V t ) + p t (X t , R Z i V t ) dt = C + T f 0 ω(t) d i=1 E pt G 2 (s θ i (X t , V t , t)) + G 2 (s θ i (X t , R Z i V t , t)) -2G(s θ i (X t , V t , t)) dt.The forward PDMPs that we consider can all be simulated exactly by solving the integral (13), at least when the limiting distribution is the multivariate standard normal. This is possible because for each process we can easily simulate the random event times as well as their deterministic dynamics.The general procedure for the simulation of a time-homogeneous PDMP up to a fixed time horizon T is given in Algorithm 1. In the remainder of the section we give additional details on the simulation of each process.The case of RHMC is trivial, since the random events have the exponential distribution with constant parameter λ r and the deterministic dynamics are given by x t = x 0 cos(t) + v 0 sin(t) and v t = -x 0 sin(t) + v 0 cos(t), where (x 0 , v 0 ) is the initial condition. Hence all the steps in Algorithm 1 can be performed and the state (X T , V T ) can be easily obtained.ZZP: Notice the event rates of ZZP are of the form λ i (x, v) = (v i x i ) + when the stationary distribution is the standard normal. In this case we find that each coordinate of the ZZP is independent, that is the evolution of ((X t ) i , (V t ) i ) is not affected by ((X t ) j , (V t ) j ) for i, j = 1, . . . , d with i ̸ = j. Therefore we can simulate each coordinate of the ZZP in parallel following the procedure in Algorithm 1. Let us illustrate how one can obtain the next event time τ of the i-th coordinate when the process is at (x i , v i ). We have that τ solves 0 (v i x i + u) + du -E = 0 for E ∼ Exp(1). This gives the following quadratic equation for τ :τ 2 2 + v i x i τ -v i x i (-v i x i ) + - 1 2 (-v i x i ) 2 + -E = 0, which has solution τ = -v i x i + (v i x i ) 2 + + 2E.BPS: In the case of BPS one has to simulate a proposal for both the next reflection event, τ b and for the next refreshment event, τ r , then accepting the smallest of the two as event time. Obtaining the proposal for the next refreshment time is straightforward since the rate is constant. The proposal for the following reflection time can be obtained similarly to the case of ZZP, but noticing that in this case the event rate is λ(x, v) = ⟨v, x⟩ + . Then τ b solves 1). Noticing that it must be τ b > ⟨v,x⟩ /|v| 2 , this gives the following quadratic equation forτ b 0 (⟨v, x⟩ + |v| 2 u) + du -E = 0 for E ∼ Exp(τ b : |v| 2 2 τ 2 b + ⟨v, x⟩τ b + 1 2 (-⟨v, x⟩) 2 + -E = 0,which has solutionτ b = -⟨v, x⟩ + ⟨v, x⟩ 2 + + 2|v| 2 E |v| 2 .Here we discuss the splitting schemes we use to discretise the backward PDMPs. For further details on this class of approximations we refer the reader to Bertazzi et al. [2023] .We have already discussed the splitting scheme DJD for RHMC in Section 2.4, and we give the pseudo-code in Algorithm 2.Algorithm 2: Splitting scheme DJD for the time reversed RHMC Initialise either from(X 0 , V 0 ) ∼ π ⊗ ν or (X 0 , V 0 ) ∼ π(dx)p θ * (dv|x, T f ); for n = 0, . . . , N -1 do ( X, V ) = φ H -δ n+1/2 (X tn , V tn ) ; t = T f -t n -δn+1 2 ; Estimate ratio: s = ν( V ) /p θ * ( V | X, t ) ; Draw proposal τ n+1 ∼ Exp(sλ r ) ; if τ n+1 ⩽ δ n+1 then Draw V ∼ p θ * ( • | X, t ); end (X tn+1 , V tn+1 ) = φ H -δ n+1/2 ( X, V ); end ZZP: For ZZP we apply the splitting scheme DJD discussed in Section 2.4, with the only difference that we allow multiple velocity flips during the jump step similarly to Bertazzi et al. [2023] . Algorithm 3 gives a pseudo-code.In the case of BPS, we follow the recommendations of Bertazzi et al. [2023] and adapt their splitting scheme RDBDR, where R stands for refreshments, D for deterministic motion, and B for bounces. We give a pseudo-code in Algorithm 4. We remark that an alternative is to use the scheme DJD for BPS, simulating reflections and refreshments in the J part of the splitting. This choice has the advantage of reducing the number of model evaluations.Backward Kernels We learn s θ by minimising the empirical counterpart of the implicit ratio matching loss ℓ I given in (7). We define such empirical loss lI as the function lI :ZZP for i = 1, . . . , d for i = 1, . . , d ← - λ Z i (t, (y, w)) = p Tf -t (R Z i w | y) p Tf -t (w | y) λ Z i y, R Z i w ← - Q Z i (y, w), (dx, dv) = δ y,R Z i w (dx, dv) RHMC ← - λ H (t, (y, w)) = λ r ν(w) p Tf -t (w | y) ← - Q H t, (y, w), (dx, dv) = p Tf -t (v | y) δ y (dx) dv BPS ← - λ B 1 (t, (y, w)) = p Tf -t (R B y w | y) p Tf -t (w | y) λ B 1 y, R B y w ← - Q B 1 (y, w), (dx, dv) = δ y,R B y w (dx, dv) ← - Q B 2 t, (y, w), (dx, dv) = p Tf -t (v | y) δ y (dx) dv ← - λ B 2 (t, (y, w)) = λ r ν(w) p Tf -t (w | y)EQUATIONwhere {τ n } N n=1 are i.i.d. samples from ω, independent of {(X n t , V n t ) t⩾0 } N n=1 , which are N i.i.d. realisations of the ZZP respectively starting at the n-th training data point X n 0 with velocity V n 0 , where{V n 0 } N n=1 are i.i.d. observations of Unif({-1, 1} d ).Algorithm 5 shows the pseudo code for our training algorithm. We remark that the simulation of the forward ZZP follows the guidelines explained in Appendix C.1.Input: Time distribution ω on [0, T f ], model s θ while θ has not converged do Get random data batch{X n 0 } B n=1 ; Sample {τ n } B n=1 ∼ ω ⊗B ; Sample {V n 0 } B n=1 ∼ Unif({±1} d ) ⊗B ; for n = 1 to B do Simulate (X n τ n , V n τ n ) by running a ZZP from (X n 0 , V n 0 ); Compute loss lI (θ) ← LI (s θ , (X n τ n , V n τ n , τ n ) B n=1) ; Perform optimization step on lI (θ);Output: trained model s θ * A computationally efficient model for the ZZP. The computation of LI in each iteration of Algorithm 5 requires evaluating d + 1 times the model s θ for each data-point. Indeed, the loss function (7) requires evaluating s θ for each component flip of the velocity vector. Hence the computational cost of the algorithm scales linearly in the data dimension d.In order to overcome this computational burden we build an alternative model leveraging the following observation:EQUATIONwhere v -i denotes the vector containing all components of v other than the i-th. The approximation in ( 21) is motivated by the fact that we expect the components of the velocity to be nearly independent conditional on the position vector. This is because the position of the process holds most of the information regarding the velocity of each coordinate. Under (21) we find that a good approximation for the ratio r Z i (x, v, t) is given byr Z i (x, v i , t) := p t (-v i |x) p t (v i |x) .Hence it is reasonable to estimate the simplified ratios r Z i rather than the true r Z i (x, v, t). In particular this can be achieved with a model s θ that requires the current position and time, and only the i-th velocity component rather than the full vector v.Following the reasoning above we build an alternative model which takes as input the position vector and the time variable and outputs a 2d-dimensional vector, that is {s θ : R d × [0, T f ] → R 2d + : θ ∈ Θ}. We introduce the following notation for the output of the neural network s θ :EQUATIONwhere s θ + and s θ -are both vectors in R d and denote the two, d-dimensional blocks in the output of s θ . The model will be trained in such a way that s θ +,i (x, t), that is the i-th component of the vector s θ + (x, t), approximates r Z i (x, +1, t), i.e. in the case v i = +1. Similarly, we estimate r Z i (x, -1, t) with s θ -,i (x, t). Let us now describe how we can train this model in such a way that the computational cost remains constant in the dimensionality of the data. For any w ∈ {1, -1} d , we introduce the projection operator Π i w defined for anyEQUATIONwhere we have sign(w i ) = + when w i = +1 and sign(w i ) =when w i = -1. In words, Π i w s θ (y, t) selects either s θ + or s θ -in ( 22) based on sign(w i ) and returns the estimate of the i-th ratio at (y, t) corresponding to the velocity w i . Using the operator Π i w we can re-formulate the loss (7) aslI (θ) = T f 0 dt ω(t) d i=1 E G 2 (Π i Vt s θ (X t , t)) + G 2 (Π i -Vt s θ (X t , t)) -2G(Π i Vt s θ (X t , t)) .The associated empirical loss is lI :EQUATIONwhere {τ n } N n=1 and {(X n t , V n t ) t⩾0 } N n=1 are obtained as in ( 20). This simplified model can then be trained using the same procedure shown in Algorithm 5, but where the loss above is used instead of the loss (20). The computational cost for the training of this model is clearly constant in the data dimension d, since only a single evaluation of the model s θ is needed for each data-point.In the case of RHMC and BPS both the backward rate and backward jump kernel are characterised by p t (v|x). We introduce a parametric family of conditional probability distributions {p θ : θ ∈ Θ} of the form (x, v, t) → p θ (v|x, t), where Θ ⊂ R d θ , which we model with normalising flows (NFs) [Papamakarios et al., 2021] , as it permits both obtaining an estimate of the density, at a given state and time, and also to generate samples according to this distribution. To train our model, we rely on the maximum likelihood population loss method, as derived in (8), and use its empirical counterpart (9). The final BPS and RHMC training algorithm is given in Algorithm 6. The simulation of the forward BPS and RHMC follows the methods listed in Appendix C.1.We provide a short comparison of the computational complexity of our PDMP generative methods with diffusion models, as each generative method admits the same design: a fixed forward process {X t } 0⩽t⩽T f ran from time 0 to T f , with X T f approximately distributed as a standard Gaussian N (0, I d ) (using the Variance-Preserving process in the case of diffusion models [Song et al., 2021] ), a corresponding backward process { ← -X t } 0⩽t⩽T f , and a generative process { ← -X θ t } 0⩽t⩽T f being an approximation to the true backward process, initialized from ← -X θ 0 ∼ N (0, I d ).Computational differences in the generative processes Since we are using the splitting scheme for our PDMPs, the backward processes admits a computational complexity growing linearly with the chosen number of backward steps N , alike diffusion models [Song et al., 2021] . The latter models require only a single network inference per backward step. However, as measured in Table 2 , each backward step is costlier in the case of our PDMP samplers. We provide an explanation for each sampler:• ZZP Using the simplified model introduced in Appendix D.1 and Algorithm 3, only one network inference is required per backward step, to approximate the backward rate. However, this approach requires a slight modification to the neural network architecture, involving a doubled channel output and adjustments for the element selection mechanisms (projection operator ( 23)), resulting in a higher cost per step than diffusion models.• RHMC Using Algorithm 2, one likelihood computation from the normalizing flow is needed, to approximate the backward rate. Additionally, at most one random variable must be drawn from the normalizing flow, which approximates the backward kernel.• BPS Using Algorithm 4, three likelihood computations are required from the normalizing flow, to approximate backward rates. Moreover, at most two random variables need to be sampled from the normalizing flow, which approximates the backward kernel. The final cost per step depends on the chosen normalizing flow architecture, but is higher for BPS than for all other samplers.It should be noted that, in our experiments, each generative PDMP model requires less backward steps than the diffusion model, leading to a smaller overall computational time for equal generation quality, as showed in Figure 3 .E Discussion and proof for Theorem 1 E.1 Discussion on H2In this section we discuss H2 in the case of ZZP, BPS, and RHMC. For all three of these samplers, existing theory shows convergence of the formEQUATIONwhere V : R 2d → [1, ∞) is a positive function and ∥µ∥ V := sup |g|⩽V |µ(g)| is the V -norm. When the initial condition of the process is µ ⋆ ⊗ ν, we obtain the bound∥µ ⋆ ⊗ νP t -π ⊗ ν∥ V ⩽ C ′ e -γt µ ⋆ ⊗ ν(V ),which translates to a bound in TV distance, since we assume V ⩾ 1. Conditions on π ensuring (26) can be found for ZZP in Bierkens et al. [2019b] , for BPS in Deligiannidis et al. [2019] , Durmus et al. [2020] , and for RHMC in Bou-Rabee and Sanz-Serna [2017] . Observe that we can set the constantC in H2 to C = C ′ µ ⋆ ⊗ ν(V ). Clearly, C is finite whenever µ ⋆ ⊗ ν(V ) < ∞. Since V is such that lim |z|→∞ V (z) = +∞,showing C is finite requires suitable tail conditions on the initial distribution µ ⋆ ⊗ ν. Inspecting the results of the papers mentioned above, one can verify that µ ⋆ ⊗ ν(V ) < ∞ when π is a multivariate standard Gaussian distribution as long as: (i) the tails of µ ⋆ are at least as light as those of π for ZZP and BPS, (ii) µ ⋆ has finite second moments for RHMC.First notice that ∥µEQUATIONhence we focus on bounding the right hand side. Under our assumptions, the forward PDMP(X t , V t ) t∈[0,T f ] admits a time reversal that is a PDMP ( ← - X t , ← - V t ) t∈[0,T f ] with characteristics ( ← - Φ , ← - λ , ← - Q ) satisfying the conditions in Proposition 1. Therefore, it holds µ ⋆ ⊗ ν = L( ← - X T f , ← - V T f) and so ( 27) can be written as∥µ ⋆ -L(X T f )∥ TV ⩽ ∥L( ← - X T f , ← - V T f ) -L(X T f , V T f )∥ TV , We find ∥ ← - Q Z (t, (x, v), •) -QZ (t, (x, v), •)∥ TV = sup A d i=1 1 (x,R Z i v)∈A - λ Z i (t, (x, v)) ← - λ Z (t, (x, v)) - λZ i (t, (x, v)) λZ (t, (x, v)) ⩽ sup A d i=1 1 (x,R Z i v)∈A ← - λ Z i (t, (x, v)) -λZ i (t, (x, v)) ← - λ Z (t, (x, v)) + λZ i (t, (x, v)) ← - λ Z (t, (x, v)) - λZ i (t, (x, v)) λZ (t, (x, v)) ⩽ d i=1 | ← - λ Z i (t, (x, v)) -λZ i (t, (x, v))| ← - λ Z (t, (x, v)) + | λZ (t, (x, v)) - ← - λ Z (t, (x, v))| ← - λ Z (t, (x, v))In the last inequality we used that λZ is non-negative. Therefore we findg 1 t (x, v) ⩽ 1 2 d i=1 | ← - λ Z i (t, (x, v)) -λZ i (t, (x, v))| + 1 2 | λZ (t, (x, v)) - ← - λ Z (t, (x, v))| ⩽ d i=1 | ← - λ Z i (t, (x, v)) -λZ i (t, (x, v))|.Noticing that| ← - λ Z i (t, (x, v)) -λZ i (t, (x, v))| = |r Z i (x, v, t) -s θ i (x, v, t)| λ Z i ((x, R Z i v).we findg t (x, v) ⩽ 2 d i=1 |r Z i (x, v, t) -s θ i (x, v, t)| λ Z i ((x, R Z i v)Finally, we use the inequality 1 -e -z ⩽ z, which holds for z ⩾ 0, to conclude thatE 1 -exp - T f 0 g T f -t (X t , V t )dt ⩽ 2 d i=1 E T f 0 |r Z i (X t , V t , T f -t) -s θ i (X t , V t , T f -t)| λ Z i (X t , R Z i V t )dt .We run our experiments on 50 Cascade Lake Intel Xeon 5218 16 cores, 2.4GHz. Each experiment is ran on a single CPU and takes between 1 and 5 hours to complete, depending on the dataset and the sampler at hand.For each forward PDMP, we take a time horizon T f equal to 5, and set the refreshment rate λ r to 1. For training, we choose the uniform distribution Uniform([0, T f ]) as the time distribution ω. For the simulation of backward PDMPs with splitting schemes, we use a quadratic schedule for the time steps, that is (δ n ) n∈{1,...,N } given by δn = T f × (( n /N) 2 -( n-1 /N) 2 ).For i-DDPM, we follow the design choices introduced in Nichol and Dhariwal [2021] and in particular we use the variance preserving (VP) process, the cosine noise schedule, and linear time steps. Detailed setup For ZZP and i-DDPM we use a neural network consisting of eight timeconditioned multi-layer perceptron (MLP) blocks with skip connections, each of which consisting of two fully connected layers of width 256. The time variable t passes through two fully connected layers of size 1 × 32 and 32 × 32, and is fed to each time conditioned block, where it passes through an additional 32 × 64 fully connected layer before being added element-wise to the middle layer. The model size is 6.5 million parameters. For ZZP, we apply the softplus activation function x → 1 /β log(1 + exp(βx)) to the output of the network, with β = 1, to constrain it to be positive and stabilise behaviour for outputs close to 1.In the case of RHMC and BPS, we use neural spline flows [Durkan et al., 2019] to model the conditional densities of the forward processes, as it shows good performance among available architectures. We leverage the implementation from the zuko package [Rozet et al., 2022] . We set the number of transforms to 8, the hidden depth of the network to 8 and the hidden width to 256. To condition on x, t, we feed them to three fully connected layers of size d × 8, 8 × 8 and 8 × 8, where d is either the dimension of X t , or d = 1 in the case of the time variable. The resulting vectors are then concatenated and fed to the conditioning mechanism of zuko. The resulting model has 3.8 million parameters.We take advantage of the approaches described in Section 2.3 to learn the characteristics of the backward processes.All experiments are conducted using PyTorch [Paszke et al., 2019] . The optimiser is Adam [Kingma and Ba, 2015] with learning rate 5e-4 for all neural networks.In Table 4 we show the accuracy in terms of the refreshment rate, while in Table 5 we show different choices of the time horizon. In both cases, we consider the Gaussian mixture data and we use the 2-Wasserstein metric to characterise the quality of the generated data. Figure 5 shows the generated data by the best model for each process.We consider the task of generating handwritten digits training the ZZP on the MNIST dataset. We use the simplified loss function given in Equation ( 24) in Appendix D.1 and use the simplified model and loss. In Figure 6 we show promising results obtained with the design choices described previously in Appendix F.1, apart from the differences that follow. The optimiser is Adam [Kingma and Ba, 2015] with learning rate 2e-4. We use a U-Net following the implementation of Nichol and Dhariwal [2021] , where the hidden layers are set to [128, 256, 256, 256] , where we fix the number of residual blocks to 2 at each level, and add self-attention block at resolution 16 × 16, with 4 heads. We duplicate the channel of the penultimate layer, and make each copy go through separate MLPs to obtain the two vectors (s 22)). We use an exponential moving average with a rate of 0.99. At every layer, we use the silu activation function, while we apply the softplus to the output of the network, with β = 0.2. We train the model for 40000 steps with batch size 128.θ + (•, •), s θ -(•, •)) ∈ R 2d + (as in (• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 


Decision: The reviewers, some of whom are experts in topics intimately related to PDMPs such as MCMC methods, and others of whom are experts in generative modeling, unanimously recommend acceptance. I concur.

Reviews: review: paper_summary: paper_summary: The paper proposes using Piecewise Deterministic Markov Processes (PDMPs) for generative modelling applications, by using the property that PDMPs also admit time reversals that themselves are PDMPs. There are three major contributions in my understanding - 

By characterizing certain families of PDMPs, i.e. Zig-Zag processes (ZZP, the Bouncy Particle Sampler (BPS), and the Randomised Hamiltonian Monte Carlo (RHMC), in terms of their jump rates and kernels, this paper shows how to obtain tractable closed-form and approximations for the jump rates and kernels of the time reversed PDMPs.
Theoretically, this paper then proposes a total variation bound between the “learnt” data distribution and the true data distribution when the base distribution is a Gaussian. This is a useful property, quite similar to bounds that have been proposed before in the literature (for example, for the Ornstein-Uhlenbeck process in [1, Theorem 5.2.]).
Finally, the paper proposes two empirical techniques to learn the time reversals akin to score-matching. First, for the ZZP process, inspired by score-matching techniques, the authors propose a ratio-matching technique. Secondly, for the BPS and RHMC processes, the authors learn normalising flows for the time reversal. They then show promising results in low-dimensional and MNIST generative modelling applications as a proof-of-concept.

[1] Dominique Bakry, Ivan Gentil, Michel Ledoux, et al. Analysis and geometry of Markov diffusion operators, volume 103. Springer, 2014.
strength_weakness: strength_weakness: Fundamentally, I think the paper lacks a convincing argument about why generative modelling with PDMPs would fundamentally be more useful than traditional generative modelling. I understand that there were some arguments made in the introduction of the paper, namely Lines 35-36 (“such as better scalability and reduced computational complexity in high-dimensional settings”). However, it is really unclear to me how this argument actually translates to the empirical score-matching (or normalising flow training) objectives that the authors formulate, vs an approach like DDPM. The experimental section is quite lacking in details and comparisons about how the PDMP approach improves along any number of axes, beyond the qualitative plots. For example, I can think of many axes of improvement that could be discussed - 
sample efficiency (how many training datapoints are needed to learn the time-reversal given that the process is partly deterministic), 
mixing rates towards the Gaussian for their time reversal. Usually, SDEs such as the Ornstein-Uhlenbeck process are quite quick at mixing towards a Gaussian, making them quite nice to use when reversing a Gaussian distribution as the base. For partly deterministic processes, is this easier or harder to do?
Are there any comparisons to regular Markov process methods that can show that having an irreversible Markov process is beneficial here? I believe that this is a big factor in why PDMPs are alluring, and their irreversibility makes them mix faster and use less data [3]. Any experiments showing sample efficiency and mixing rates would be really beneficial here.

I am worried that there are many subtleties in the training and sampling procedures of diffusion models, and indeed there are many papers focusing solely on the empirical training tricks that can improve generative modelling, and comparing to a vanilla DDPM model doesn’t properly ablate the technique. I would be hesitant to rely on these empirical results as a surefire sign of improved modeling, which is frustrating, as theoretically, the paper does seem to point to this being the case, and I really do want to believe. 

Furthermore, I think the paper could benefit from being a lot clearer about the specific advantages of PDMPs vs other stochastic processes for generative modeling. This is barely mentioned, but does form the crux of the empirical results. This made it difficult for me to read through the theoretical developments, proofs and theorems without knowing the reason why we would really want to do this in the first place.

I also think the paper can also benefit from being more explicit in how their developed score matching and normalizing flow training differs from traditional methods (maybe an algorithm block), as this would be something really interesting to practitioners looking to adopt existing codebases to using PDMPs instead.

[3] 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, 1288–1320. DOI: 10.1214/18-AOS1715
questions: questions: I’ve summarised a few of my questions in the Weaknesses section, but I list a few more questions here I would love some clarity on - 

You mention that [2] originally suggested time reversal theorems for jump Markov processes. I found it really difficult to figure out what the exact distinctions are between the theoretical formulations of [2], and the specific novelty of this paper. I’d love clarifications on what the additional mathematical frameworks introduced are, that adopt this framework to PDMPs.
I would really like a better discussion on Equation 9, talking about what the different components of the inequality correspond to, and if there are any insights to gain about why PDMPs perform better or differently to something traditional like an OU process.
In Line 238-239, you mention that it’s easy to subsample the loss across dimensions. Is this still an unbiased, consistent estimator for the continuous time loss?
You mention that PDMPs typically have several types of jump rates and kernels, and provide Equation 2 as a parametric family of possible jumps and kernels. Is this the only possibility for constructing jumps and kernels?
At the end of each section describing one of the 3 PDMPs, the authors mention that there has been work showing exponential convergence to invariant distributions for all the PDMPs. How does this exponential convergence compare to the OU process? Is there actually faster convergence due to the inherent irreversibility of PDMPs?
limitations: limitations: Yes, they have addressed any potential limitations.
paper_summary: paper_summary: This paper considers the development of generative models based on piecewise deterministic Markov processes. The key idea proposed in the paper is to use piecewise deterministic Markov processes instead of diffusions as the "noising process" of the generative model. This relies on the fact that time reversals of PDMPs are themselves PDMPs. Three specific instances of PDMPs are considered. The authors also derive a bound (in total variation distance) between the data distribution and the distribution of the generative model. The methodology is illustrated with some simple experiments.
strength_weakness: strength_weakness: To me, the descriptions of approximating the process characteristics with normalizing flows are unclear. This part should be written more with more details, perhaps in the supplement, as this is core to being able to reproduce the results. I do appreciate that the authors provided a description of the experiment in E.1, but it is not enough to put things together. I would be interested in replicating at least the simple experiment, but I don't think I can do it as the paper stands.
questions: questions: Are there particular data characteristics that these processes would be most suitable for? It seems that images with sharp boundary transitions are better modeled by the authors proposed approach. Would you say that is the case?

Do these methods have issues like being only able to generate from one mode if the data distribution is multimodal?

Why does RHMC do so well for such a small number of steps?
limitations: limitations: Yes.
paper_summary: paper_summary: This interesting paper on the popular topic of generative models introduce a new family of generative models which builds on the so-called piecewise deterministic Markov process (Zig-Zag process, Bouncy Particle Sampler, Randomised Hamiltonian Monte Carlo). In contrast to many of the existing models this family is not based on diffusion models. The paper includes a through analysis of the construction 
and it propose training procedures and methods for approximate simulation of the reverse process.
strength_weakness: strength_weakness: * Missing real-world examples
* The is a big jump in  the style of writing between Section 1 and 2. Do not get me wrong here, the technical developments are most interesting, but many readers would be helped by a more gentle transition between these sections. Space for this can be created by moving more of the technical details into the supplemental material.
questions: questions: * You provide convincing examples on simple synthetic datasets. Is you family of methods able to generate data also in more challenging real-world examples? 
* When you compare the results of the various methods you show results for "steps". Is the a fair way of comparing the methods? Would it not be better to have some kind of computational cost on the x-axis? My reasoning here is that a "step" is not a well defined concept and it could mean very different things to different methods. In any case it would be really interesting to see the performance with the computational cost on the x-axis.
limitations: limitations: Real-world examples missing.
paper_summary: paper_summary: The paper explores generative models utilizing PDMPs, a type of stochastic process characterized by deterministic motion interspersed with random jumps at random times. These models offer an alternative approach to diffusion-based generative models, which have become very popular in the AI community in recent years. The authors focus on three specific PDMPs: the Zig-Zag process (ZZP), the Bouncy Particle Sampler (BPS), and Randomised Hamiltonian Monte Carlo (RHMC). The authors leverage the existing literature on the time reversals of Markov jump processes, characterizing the time reversal of any PDMP under appropriate conditions. They show that the time reversal of a PDMP is itself a PDMP with characteristics related to the original process. The authors also specifically outline the characteristics of the time-reversal processes for ZZP, BPS, and RHMC. The jump rates and kernels of these time reversals admit explicit expressions based on the conditional densities of the PDMPs before and after jumps.

The authors provide bounds on the total variation distance between the data distribution and the distribution of their generative models, considering errors from the approximation of the backward PDMP’s characteristics and its initialization from the forward process’s limiting distribution. Some initial but promising numerical simulations on simple toy distributions are presented, showcasing the potential of PDMP-based generative models. The results support further investigation into this class of models on more challenging data structures.
strength_weakness: strength_weakness: Originality: While the introduction of PDMP-based generative models is novel, the paper primarily focuses on theoretical aspects and toy datasets. There is a limited exploration of how these models can be applied to more complex, real-world problems, which could showcase their true originality and practical utility. Although PDMPs offer a new approach, the paper does not extensively compare these models with a wide variety of existing generative models, even on the simple toy examples that were considered here. This makes it difficult to fully appreciate the originality and benefits of PDMPs over other state-of-the-art methods.

Quality: The paper relies on several technical assumptions and conditions (e.g., H3, H4). While these are necessary for the theoretical results, their practical applicability might be limited and the authors do not discuss whether or not these assumptions hold for the toy examples which they consider. The paper could be strengthened by discussing the feasibility of these assumptions in real-world scenarios. The empirical validation is primarily limited to simple toy distributions. While these are useful for initial validation, the lack of experiments on more complex datasets (e.g., image or text data) reduces the overall impact and persuasiveness of the empirical results. It would be beneficial to include comparisons on more challenging benchmarks.

Clarity: The paper uses dense mathematical notation and detailed proofs, which may be challenging for readers who are not specialists in stochastic processes or PDMPs. This is going to be challenging for the authors as the high-level of technical detail provided in the paper does lead to a very robust paper. However, perhaps more intuitive explanations or visual aids could help make the content more accessible, if not in the main paper then in the appendix. The practical implementation details, particularly regarding the training procedures and simulation methods, are somewhat sparse (even though there are more details in the appendix). Providing a step-by-step guide or pseudocode could help practitioners better understand how to apply the proposed methods. Given the space constraints, this would have to be added to the appendix. This is covered in the case of splitting schemes, but could perhaps be modified to be more user-friendly to people new to this area of research. 

Significance: The paper’s significance is somewhat limited by the focus on theoretical and synthetic examples. Without demonstrating the effectiveness of PDMP-based models on real-world data, it is challenging to gauge their practical significance and potential impact in applied settings. While the paper claims that PDMPs offer better scalability and reduced computational complexity, there is limited empirical evidence to support these claims. Benchmarking the computational performance against existing generative models would provide a clearer picture of the advantages and limitations in terms of scalability and efficiency.
questions: questions: - I've outlined quite a few points in the weaknesses section that the authors may wish to address as questions.
- Could the authors include computational time with their simulations?
- Under H2, the authors state that the condition is satisfied when "the tails of $\mu_*$ are sufficiently light." Could this be clarified as to what "sufficiently light" means? I'm interpreting this as lighter than Gaussian. 
- A small query, in Proposition 2, on Line 166, it states that $\mu_0^X \otimes \mu_0^V$ on $\mathbb{R}^2d$. But how does this work for ZZP when $V=[-1,+ 1]^d$ (sorry - I can't get the curly brackets to work)?
limitations: limitations: NA
paper_summary: paper_summary: This paper proposes a new generative model on continuous random variables similar to diffusion-based generative models, but this method uses a different family of perturbation schemes. In particular, for the forward processes, which transform the clean data to a stationary distribution, the authors propose using piecewise deterministic Markov processes (PDMPs), which consist of a combination of discrete diffusions and ODEs. Specifically, the initial value, i.e., data, is transformed via simple ODEs, whose velocity term is often chosen to be a constant. However, the velocity term (also called vector fields) is updated at any random time with a constant rate $\lambda$. The occurrence of a new event is independent, and the time interval of the new event since the last one is exponentially distributed. Moreover, the new velocity will be assigned by user-chosen schemes; for example, one can flip the sign of the velocity, or another one can sample from the standard normal distribution. Thus, while the forward process transforms the data with a simple ODE, it randomly changes the ODEs, hence, the title of the paper.

In addition, the paper shows that under some mild assumptions, the time-reversed processes of the PDMPs are also PDMPs, and their transition probability distribution and update rate can be written in terms of the forward processes. Therefore, we can consider modeling such terms with parametric models similar to popular diffusion-based generative models.
The paper proposes the explicit ratio matching objective function, as in Equation 6, to train such terms.

Moreover, the authors suggest using three popular PDMPs for the forward process: Zig-Zag process (ZZP), bouncy particle sampler (BPS), and randomized Hamiltonian Monte Carlo (RHMC). Consequently, the paper shows the time-reversed PDMPs and potential parameterizations that fit each case.

Finally, the paper demonstrates the efficacy of the proposed method through a few toy experiments.

----------------------------------
Update the rating from 5 to 7 after the authors' rebuttal.
strength_weakness: strength_weakness: The proposed method is novel and interesting. However, the paper didn't characterize some potential drawbacks well. For example, due to the nature of PDMPs, the training of the proposed method inherits the problem of discrete diffusion models. In particular, the ratio matching should be done for each dimension independently, and this costs a lot for high-dimensional data compared to the denoising score matching and others for previous methods. This problem could be the reason that the current submission didn't include the experiments on popular high-dimensional datasets. In this regard, the paper suggests that Monte Carlo estimates can reduce the computational cost, but it trades off the variance of the learning signals, which would be critical for large-scale experiments.

I consider that a new method doesn't always need to be better than previous methods. Nevertheless, proper information about the proposed method in practice would be important for potential readers to evaluate the significance of the paper. Thus, some additional discussions related to this problem need to be added.

Moreover, given that many alternative choices exist other than the proposed method, in-depth discussions on the proposed method's merits would be appreciated for evaluating its importance.

In addition, I find that the presentation should be improved. For example, many statements contain multiple prepositional phrases, which are difficult to parse and understand.
questions: questions: N/A
limitations: limitations: N/A

Comments: title: Response to authors' rebuttal
comment: Thank you for your responses, they've definitely cleared up a few doubts I had regarding the paper. I think the paper is in a great spot with the addition of your rebuttal figures (1 and 2 in the pdf), and I am quite confident with my initial score, and will be happy to advocate for a paper acceptance, pending other reviewers' responses. Overall, I think the rebuttal responses were well-explained, and I thank the authors for the additional effort, and I'm confident this will help make the paper stronger.
Thank you to the authors for providing a comprehensive rebuttal. I am satisfied that the authors have addressed my questions and I thank them for the incorporating my feedback into their revised paper.
Thank you for these clarifications. I am satisfied with the responses.