A piecewise deterministic Monte Carlo method for diffusion bridges

We introduce the use of the Zig-Zag sampler to the problem of sampling conditional diffusion processes (diffusion bridges). The Zig-Zag sampler is a rejection-free sampling scheme based on a non-reversible continuous piecewise deterministic Markov process. Similar to the L\'evy-Ciesielski construction of a Brownian motion, we expand the diffusion path in a truncated Faber-Schauder basis. The coefficients within the basis are sampled using a Zig-Zag sampler. A key innovation is the use of the fully local Algorithm for the Zig-Zag sampler that allows to exploit the sparsity structure implied by the dependency graph of the coefficients and by the subsampling technique to reduce the complexity of the algorithm. We illustrate the performance of the proposed methods in a number of examples.


Introduction
Diffusion processes are an important class of continuous time probability models which find applications in many fields such as finance, physics and engineering.They naturally arise by adding Gaussian random perturbations (white noise) to deterministic systems.We consider diffusions described by a one-dimensional stochastic differential equation of the form dX t = b(X t )dt + dW t , X 0 = u, where (W t ) t≥0 is a driving scalar Wiener process defined in some probability space and b is the drift of the process.The solution of equation (1), assuming it exists, is an instance of
If ξ is standard normal and {ξ i,j } is a sequence of independent standard normal random variables (independent of ξ), then converges almost surely on [0, T ] (uniformly in t) to a Brownian motion as N → ∞ (see e.g.Section 1.2 of McKean, 1969).The basis formed by φ and {φ i,j } is known as the Faber-Schauder basis (see Figure 1).The larger i, the smaller the support of φ i,j , reflecting that higher order coefficients represent the fine details of the process.A Brownian bridge starting in u and ending in v can be obtained by fixing ξ = v/ √ T and adding the function φ(t)u = (1 − t/T )u to (2).By sampling ξ N := (ξ 0,0 , ξ 1,0 , ..., ξ N,2 N −1 ) (which in this case are standard normal), approximate realisations of a Brownian bridge can be obtained.

Zig-Zag sampler for diffusion bridges
Let Q u denote the Wiener measure on C[0, T ] with initial value X 0 = u (cf.section 2.4 of Karatzas and Shreve, 1991) and let P u denote the law on C[0, T ] of the diffusion in (1).Under mild conditions on b, the two measures are absolutely continuous and their Radon-Nikodym derivative dP u dQ u is given by the Girsanov formula.Denote by P u,v T and Q u,v T the measures of the diffusion bridge and the Wiener bridge respectively, both starting at u and conditioned to hit a point v at time T .Applying the Bayes' law for conditional expectations (Klebaner, 2005, Chapter 10) we obtain: dP u,v T dQ u,v T (X) = q(0, u, T, v) p(0, u, T, v) where p and q are the transition densities of X under P, Q respectively so that for s < t, p(s, x, t, y)dy = P (X t ∈ dy | X s = x).As p is intractable, the Radon-Nikodym derivative for the diffusion bridge is only known up to proportionality constant.The main idea now consists of rewriting the Radon-Nikodym derivative in (3), evaluating it in X N and running the Zig-Zag sampler for ξ N targeting this density.Technicalities to actually get this to work are detailed in Section 3. A novelty is the introduction of a local version of the Zig-Zag sampler, analogously to the local bouncy particle sampler (Bouchard-Côté, Vollmer, and Doucet, 2015).This allows for exploiting the sparsity in the dependence structure of the coefficients of the Faber-Schauder expansion efficiently, resulting in a reduction of the complexity of the algorithm.The methodology we propose is derived for one dimensional diffusion processes with unit diffusivity.However, diffusions with state-dependent diffusivity can be transformed to this setting using the Lamperti transform (an example is given in Subsection 5.3).In Subsection 6.1 we generalize the method to multivariate diffusion processes with unit diffusivity, assuming the drift to be a conservative vector field.

Contributions of the paper
The Faber-Schauder basis offers a number of attractive properties: (a) The coefficients of a diffusions have a structural conditional independence property (see Section 4 and Appendix A) which can be exploited in numerical algorithms to improve their efficiency.
(b) A diffusion bridge is obtained from the unconditioned process by simply fixing the coefficient ξ.
(c) It will be shown (see for example Figure 8) that the non-linear component of the diffusion process is typically captured by coefficients ξ ij in equation ( 2) for which i is small.This allows for a low dimensional representation of the process and yet a good approximation.Therefore, the approximation error caused by leaving out fine details is equally divided over [0, T ], contrary to approaches where a proxy for the diffusion bridge is simulated by Euler discretisation of an SDE governing its dynamics.In the latter case, the discretisation error accumulates over the interval on which the bridge is simulated.
(d) It is very convenient from a computational point of view as each function is piecewise linear with compact support.
We adopt the Zig-Zag sampler (Bierkens, Fearnhead, and Roberts, 2019) which is a sampler based on the theory of piecewise deterministic Markov processes (see Fearnhead et al., 2018, Bouchard-Côté, Vollmer, and Doucet, 2015, Andrieu and Livingstone, 2019, Andrieu et al., 2018).The main reasons motivating this choice are: (a) The partial derivatives of the log-likelihood of a diffusion bridge measure usually appear as a path integral that has to be computed numerically (introducing consequently computational burden derived by this step and its bias).The Zig-Zag sampler allows us to replace the gradient of the log-likelihood with an unbiased estimate of it without introducing bias in the target measure.This is done in Subsection 4.4 with the subsampling technique which was presented in Bierkens, Fearnhead, and Roberts (2019) for applications for which the evaluation of the log-likelihood is expensive due to the size of the dataset.(c) The method is a rejection-free sampler, differing from most of the methodologies available for simulating diffusion bridges.
(d) The Zig-Zag sampler is defined and implemented in continuous time, eliminating the choice of tuning parameters appearing for example in the proposal density of the Metropolis-Hastings algorithm.This advantage comes at the cost of a more complicated method which relies upon bounding from above rates which are model specific and often difficult to derive (see Section 5 for our specific applications).
(e) The process is non-reversible: as shown, for example, in Diaconis, Holmes, and Neal (2000), non-reversibility generally enhances the speed of convergence to the invariant measure and mixing properties of the sampler.For an advanced analysis on convergences results for this class of non-reversible processes, we refer to the articles Andrieu and Livingstone (2019) and Andrieu et al. (2018).
The local Zig-Zag sampler relies on the conditional independence structure of the coefficients only.This translates to other settings than diffusion bridge sampling, or other choices of basis functions.For this reason, Section 4 describes the algorithms of the sampler in their full generality, without referring to our particular application.A documented implementation of the algorithms used in this manuscript can be found in Schauer and Grazzi (2020).

Outline
In Section 2 we set some notation and recap the Zig-Zag sampler.In Section 3 we expand a diffusion process in the Faber-Schauder basis and prove the aforementioned conditional dependence.The simulation of the coefficients ξ N presents some challenges as it is high dimensional and its density is expressed by an integral over the path.We give two variants of the Zig-Zag algorithm which enables sampling in a high dimensional setting.In particular, in Section 4 we present the local and fully local Zig-Zag algorithms which exploit a factorization of the joint density (Appendix A) and a subsampling technique which, in this setting, is used to avoid the evaluation of the path integral appearing in the density (which otherwise would severely complicate the implementation of the sampler).In Section 5 we illustrate our methodology using a variety of examples, validate our approach and compare the Zig-Zag sampler with other benchmark MCMC algorithms.We conclude by sketching the extension of our method to multi-dimensional diffusion bridges, carrying out an informal scaling analysis and providing several remarks for future research (Section 6 and Section 7).

Preliminaries
Throughout, we denote by ∂ i the partial derivative with respect to the coefficient ξ i , the positive part of a function f by (f ) + , the ith element and the Euclidean norm of a vector x respectively by [x] i and x .The cardinality of a countable set A is denoted by |A|.

Notation for the Faber-Schauder basis
To graphically illustrate the Faber-Schauder basis, a construction of a Brownian motion with the representation of the basis functions is given in Figure 1.The Faber-Schauder functions are piecewise linear with compact support.The length of the support and the height of the function is determined by the first index while the second index determines the location.All basis functions with first index i are referred to as level i basis functions.For convenience, we often swap between double and single indexing of Faber-Schauder functions.
Denote the double indexing with (i, j) and the single indexing with n.We go from one to the other through the transformations where • denotes the floor function.The basis with truncation level N has M := 2 N +1 − 1 coefficients.Let ξ N denote the vector of coefficients up to level N , i.e.
and let X ξ N := X N when we want to stress the dependencies of X N on the coefficients ξ N .Using double indexing, we denote by S i,j = supp φ i,j .Coefficients q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.4 0.8 1.2 0.00 0.25 0.50 0.75 1.00

Approximated process
Figure 1: Lévy-Ciesielski construction of a Brownian motion on (0, 1).On the left the Faber-Schauder basis functions up to level N = 3, on the top-right the values of the corresponding coefficients located at the peak of their relative FS basis function and on the bottomright the resulting approximated Brownian path X N (black line) compared with a finer approximation (red line).The truncated sum defines the process in 2 N +1 + 1 finite dyadic points (black dots) with linear interpolation in between points.A finer approximation corresponds to Brownian fill-in noise between any two neighboring dyadic points.

The Zig-Zag sampler
A piecewise deterministic Markov process (Davis, 1993) is a continuous-time process with behaviour governed by random jumps at points in time, but deterministic evolution governed by an ordinary differential equation in between those times (yielding piecewise-continuous realizations).If the differential equation can be solved in closed form and the random event times can be sampled exactly, then the process can be simulated in continuous time without introducing any discretization error (up to floating number precision) making it attractive from a computational point of view.
By a careful choice of the event times and deterministic evolution, it is possible to create and simulate an ergodic and non-reversible process with a desired unique invariant distribution (Fearnhead et al., 2018).The Zig-Zag sampler (Bierkens, Fearnhead, and Roberts, 2019) is a successful construction of such a processes.We now recap the intuition and the main steps behind the Zig-Zag sampler.
The one-dimensional Zig-Zag sampler is defined in the augmented space (ξ, θ) ∈ R × {+1, −1}, where the first coordinate is viewed as the position of a moving particle and the second coordinate as its velocity.The dynamics of the process t → (ξ(t), θ(t)) (not to be confused with the time indexing the diffusion process) are as follows: starting from (ξ(0), θ(0)), (a) its flow is deterministic and linear in its first component with direction θ(0) and constant in its second component until an event at time τ occurs.That is, (ξ(t), θ(t)) = (ξ(0) + tθ(0), θ(0)), 0 ≤ t ≤ τ .
(b) At an event time τ , the process changes the sign of its velocity, i.e. (ξ(τ ), θ(τ The event times are simulated from an inhomogeneous Poisson process with specified rate The d-dimensional Zig-Zag sampler is conceived as the combination of d one-dimensional Zig-Zag samplers with rates λ i (ξ, θ), i = 1, ..., d, where the rates create a coupling of the independent coordinate processes.The following result provides a sufficient condition for the d-dimensional Zig-Zag sampler to have a particular d-dimensional target density π as invariant distribution.Assume that the target d-dimensional distribution has strictly positive density with respect to the Lebesgue measure i.e.
has π as invariant density.Condition (5) is derived in the supplementary material of Bierkens, Fearnhead, and Roberts (2019).Condition (5) is equivalent to for some γ i (ξ) ≥ 0. Throughout, we set γ i (ξ) = 0 because generally the algorithm is more efficient for lower Poisson event intensity (see for example Andrieu and Livingstone, 2019, Subsection 5.4).Assume the target density is π(ξ) = cπ(ξ).The process targets the specific distribution function through the Poisson rate λ which is a function of the gradient of ξ → ψ(ξ) = − log(π(ξ)), so that any proportionality factor of the density disappears.Throughout we refer to the function ψ as the energy function.As opposed to standard Markov chain Monte Carlo methods, the process is not reversible and it is defined in continuous time.
Example 2.1.Consider a d-dimensional Gaussian random variable with mean µ ∈ R d and positive definite covariance matrix Σ ∈ R d×d .Then Notice that if Σ is diagonal, then λ k (ξ, θ) = 0 whenever the process is directed towards the mean so that no jump occurs in the kth component when one of the following conditions is satisfied: In Figure 2 we simulate a realization of the Zig-Zag sampler targeting a univariate standard normal random distribution.
Algorithm 1 shows the standard implementation of the Zig-Zag sampler.Given a fixed time t ≥ 0 and a position (ξ(t), θ(t)), the first event time τ * after t is determined by taking the minimum of event times τ 1 , τ 2 , . . ., τ d simulated according to the Poisson rates λ i , i = 1, 2, ..., d.At event time τ * , the velocity vector becomes θ(τ * ) = F i * (θ(t)), with i * = arg min(τ 1 , . . ., τ d ).The algorithm iterates this step moving forward each time until the next simulated event time exceeds the final clock τ final .
Although we consider the velocities for each dimension of a d-dimensional Zig-Zag process to be either 1 or −1, these can be taken to be any non-zero values (θ i , −θ i ) for i = 1, ..., d.A finetuning of θ 1 , ..., θ N can improve the performance of the sampler.Note that the only challenge in implementing Algorithm 1 lies on the simulation of the waiting times which correspond to the simulation of the first event time of d inhomogeneous Poisson processes (IPPs) with rates λ 1 , λ 2 , ..., λ d which are functions of the state space (ξ, θ) of the process.Since the flow of the process is linear and deterministic, the Poisson rates are known at each time and are equal to To lighten the notation, we write λ i (t) := λ i (t; ξ, θ) when ξ, θ are fixed.Given an initial position ξ and velocity θ, the waiting times τ 1 , ..., τ d are computed by finding the roots for x of the equations where (u i ) i=1,2,...,d are independent realisations from the uniform distribution on (0, 1).When it is not possible to find roots of equation (7) efficiently, for example in closed form, it suffices to find upper bounds for the rate functions for which this is possible; Subsection 4.4 treats this problem for our particular setting.The linear evolution of the process and the jumps of the velocities are always trivially computed and implemented.Algorithm 1 returns a skeleton of values corresponding to the position of the process at the event times.From these values, it is straightforward to reconstruct the continuous path of the Zig-Zag sampler.Given a sample path of the Zig-Zag sampler from 0 to τ final , we can obtain a sample from the target distribution in the following way: • Denote by ξ(τ ) the value of the vector ξ at the Zig-Zag clock τ < τ final .Fixing a sample frequency ∆τ , we can produce a sample from the density π by taking the values of the random vector ξ at time τ burn-in + ∆τ, τ burn-in + 2∆τ, ...., τ final where τ burn-in is the initial burn-in time taken to ensure that the process has reached its stationary regime.Throughout the paper, we create samples using this approach.

Zig-Zag sampler for Brownian bridges
The previous subsections contain all ingredients necessary to run the Zig-Zag sampler in a finite dimensional projection of the Brownian bridge measure Q 0,v on the interval [0, T ].We fix ξ to v and run the Zig-Zag sampler for ξ N as defined in (4) targeting a multivariate normal distribution.Figure 3 shows 100 samples obtained from one sample run of the Zig-Zag sampler where the coefficients are mapped to samples paths using (2).The final clock of the Zig-Zag is set to τ final = 500 with initial burning τ burn-in = 10.
Both Brownian motion and the Brownian bridge are special in that all coefficients in the Faber-Schauder basis are independent.Of course, these processes can directly be simulated without need of a more advanced method like the Zig-Zag sampler.However, for a diffusion process with nonzero drift this property is lost.Nevertheless, we will see that when the process is expanded in the Faber-Schauder basis, many coefficients are still conditionally independent.This implies that the dependency graph of the joint density of the coefficients is sparse.We will show in Section 4 how this property can be exploited efficiently using the Zig-Zag sampler in its local version.

Faber-Schauder expansion of diffusion processes
We extend the results of Section 2 to one-dimensional diffusions governed by the SDE in (1).Although the density is defined in infinite dimensional space, in this section we justify both intuitively and formally that the diffusion can be approximated to arbitrary precision by considering a finite dimensional projection of it.
The intuition behind using the Faber-Schauder basis is that, under mild assumptions on the drift function b, any diffusion process behaves locally as a Brownian motion.Expanding the diffusion process with the Faber-Schauder functions, this notion translates to the existence of a level N such that the random coefficients at higher levels which are associated to the Faber-Schauder basis are approximately independent standard normal and independent from ξ N under the measure P.
Define the function where the first integral is understood in the Itô sense and X ≡ (X s , s ∈ [0, T ]).
For sufficient conditions for verifying that this assumption applies, we refer to Remark 3.6, Remark 3.9 and Liptser, Aries, and Shiryaev (2013), Chapter 6.
Moreover, a weak solution of the stochastic differential equation exists which is unique in law.
As we consider diffusions on [0, T ] with T fixed, we denote Z(X) := Z T (X).Due to the appearance of the stochastic Itô integral in Z(X), we cannot substitute for X its truncated expansion in the Faber-Schauder basis.Clearly, whereas the approximation has finite quadratic variation, X has not.Assuming that b is differentiable and applying Itô's lemma to the function B(x) = x 0 b(s)ds, the stochastic integral can be replaced and equation ( 8) is rewritten as where b is the derivative of b.
Definition 3.3.Let X be a diffusion governed by (1).Let X N be the process derived from X by setting to zero all coefficients of level exceeding N in its Faber-Schauder expansion (see equation ( 2)).Set We define the approximating measure P N by the change of measure where Note that the measure P u N associated to the approximated stochastic process is still on an infinite dimensional space and such that the joint measure of random coefficients ξ N is different from the one under Q u while the remaining coefficients stay independent standard normal and independent from ξ N .This is equivalent to approximating the diffusion process at finite dyadic points with Brownian noise fill-in in between every two points.We now fix the final point v T by setting ξ = v T .Define the approximated stochastic bridge with measure P u,v T N in an analogous way of equation (11), so that where The following is the main assumption made.Proof.In the following we lighten the notation by omitting the initial point u from the notation, which will be assumed fixed to u = x 0 .We wish to show that P v T N converges weakly to P v T as N → ∞.This is equivalent to showing that f dP v T N → f dP v T for all bounded and continuous functions f .Write c v T ∞ = p(0, x 0 , T, v T )/q(0, x 0 , T, v T ).By equation ( 3) and ( 9), and we have that where we used Assumption 3.1 for applying the change of measure between the conditional measures.Notice that Z N (X) = Z(X N ).The mapping X → Z(X), as a function acting on C(0, T ) with uniform norm, is continuous, since B, b, and b are continuous.Therefore, it follows from the Lévy-Ciesielski construction of Brownian motion (see Section 1.1.1)and the continuous mapping theorem that Now notice that, under conditional measures Q v T and P v T , the term B(X T )−B(X 0 ) is fixed.By the assumptions on b and b , Z is a bounded function and by dominated convergence we get that lim giving convergence to zero of the first term in (13).This implies that also the constant ∞ so that all the terms in (13) converge to 0.
We now list some technical conditions for the process to satisfy Assumptions 3.1 and 3.4.Proof.Assumption 3.4 is trivially satisfied.By Remark 3.6, also Assumption 3.1 is satisfied.
In Subsection 5.3 we will present an example where the drift b is not globally Lipschitz, yet Assumption 3.4 is satisfied.

Assumption 3.8. There exists a non-decreasing function
The above integrability condition is for example satisfied if h(|x|) = c(1 + |x|) for some c > 0.
By Lemma 3.10, below Then for a sequence of stopping times (τ k ) diverging to infinity such that (Z τ k t ) 0≤t≤T is a martingale for all k, we have Proof.The maximum M T = max 0≤t≤T X t of a Brownian motion is distributed as the absolute value of a Brownian motion and thus has density function 2 √ 2πT exp(−x 2 /(2T )), see Karatzas and Shreve (1991), Subsection 2.8.We have P(N T ≥ y) ≤ 2P(M T ≥ y) from which the result follows.
Finally we mention that Theorem 3.5 can be generalized in the following way to diffusions without a fixed end point.Proposition 3.11.If Assumption 3.4 is satisfied and B is bounded, then P N converges weakly to P.
The proof follows the same steps of the one of Theorem 3.5.In this case we need to pay attention on B, as for unconditioned process, the final point is not fixed.If B is bounded, then Assumption 3.8 is satisfied.By Remark 3.9 also Assumption 3.1 is satisfied so that we can apply Theorem 3.2 for the change of measure.Finally, by the assumptions on b and B, the function Z is bounded and by dominated convergence we get that lim

A local Zig-Zag algorithm with subsampling for high-dimensional structured target densities
In Subsection 4.4 we will show that the task of sampling diffusion bridges boils down to the task of sampling a high-dimensional vector ξ N ∈ R M under the measure P u,v T N .Define by P ξ N the distribution of the vector ξ N .Under the target measure, We take the density π to be the M -dimensional invariant density (target density) for the Zig-Zag sampler.An efficient implementation of piecewise deterministic Monte Carlo methods, including the local and fully local Zig-Zag sampler can be found in Schauer and Grazzi (2020).

Subsampling technique
In our setting, the integral appearing in the Girsanov formula (10) poses difficulties when finding the root of equation (7) and would require numerical evaluation of the integral, hence also introducing a bias.By adapting the subsampling technique presented in Bierkens, Fearnhead, and Roberts (2019) (Section 4) we avoid this problem altogether (see Subsection 4.4).In general this technique requires (a) unbiased estimators for ∂ i ψ i.e. random functions ∂ i ψi (ξ, U i ) such that for all i and ξ.These random functions create new (random) Poisson rates given by λi (t; ξ, θ; whose evaluation becomes feasible and computationally more efficient compared to the original Poisson rates given by equation ( 6).
(b) upper bounds λi : (R + × R d × {−1, +1} d ) → R + for all i = 1, ..., d such that for any point (ξ, θ) and t ≥ 0 we have As we show in Algorithm 2 and in Section 5, these upper bounds are used for finding the roots of equation (7).
Algorithm 2 gives the algorithm for the Zig-Zag sampler with subsampling.It can be proved (see Bierkens, Fearnhead, and Roberts, 2019) that the Zig-Zag sampler with subsampling has the same invariant distribution as its original and therefore does not introduce any bias.Note that we slightly modified the algorithm from Bierkens, Fearnhead, and Roberts (2019) in order to reduce its complexity.In particular it is sufficient to draw new waiting times and to save the coordinates only when the if condition at the subsampling step of Algorithm 2 is true.
Algorithm 2 d-dimensional Zig-Zag sampler with subsampling procedure ZigZag ws(τ final , ξ, θ) Recall that by the definition of λ i (see equation ( 6)), the ith partial derivative of the negative loglikelihood determines the sets N i .Now let us suppose that the first event time τ is triggered by the coordinate i so that at event time, the velocity θ i is flipped.For all λ k which are not function of this coordinate (k ∈ N i ), we have which implies that the waiting times drawn before τ , are still valid after switching the velocity i.This allows us to rescale the previous waiting time and reduce the number of computations at each step.The sets N 1 , ..., N d are connected to the factorisation of the target distribution and define its conditional dependence structure.Indeed, take a d-dimensional target distribution with the following decomposition where ξ (i) := {ξ j : j ∈ Γ i } and Γ i ⊂ {1, 2, ..., N } defines a subset of indices.We have that where the ith term in the sum is equal to 0 if k / ∈ Γ i .Since the Poisson rates ( 6) are defined through the partial derivatives, the factorisation defines the sets N 1 , ..., N d of Assumption 4.1.
Algorithm 3 shows the implementation of the local sampler which exploits any conditional independence structure so that the complexity of the algorithm scales well with the number of dimensions.
The local Zig-Zag sampler simplifies to independent one-dimensional Zig-Zag processes if the coefficients are pairwise independent coefficients, as it was the case in the example of sampling a Brownian motion or Brownian bridge (see Subsection 2.3).On the other hand, it defaults to Algorithm 1 when the dependency graph is fully connected, that is if N i = {1, . . ., d}, ∀i.

Fully local Zig-Zag sampler
Combining the subsampling technique and the local ZZ can lead to a further reduction of the complexity of the algorithm.Indeed the bounds for the Poisson rates might induce sparsity Algorithm 4 Implementation of the d-dimensional fully local Zig-Zag sampler Input: The bounds λi depend only on ξ k , θ k , for k ∈ Ni and the random Poisson rates λi (eq.( 14)) depends only on U i (the randomizing argument of ∂i ψ) and ξ k , θ k for k ∈ Ñi (U i ) procedure ZigZag fully local(τ final , ξ, θ) end while return reflection tuples ((i (l) , s (l) , ξ (l) )) l=1,...,k end procedure as λi can be function of few coordinates (see for example Subsection 5.2).This means that, after flipping θ i , λold j (τ + t) = λnew j (t) for almost all j = i making the if statement in the local step of Algorithm 3 almost always satisfied and improving the efficiency of the algorithm.This means that, after flipping θ i , we have that λold j (τ + t) = λnew j (t) for almost all j = i or, in other words, the cardinality of the set N i in the local step of Algorithm 3 is small.Furthermore, the evaluation of λi (t, ξ, θ) and λi (t, ξ, θ) for i = 1, 2, ..., d does not necessarily require to access the location of all the coordinates ξ j so that, by assigning an independent time for each coordinate and updating only the coordinates needed for the evaluation of λi and λi , the algorithm can be made more efficient.This is shown in the fully local ZZ sampler (Algorithm 4) where Ni , Ñi (U i ) define respectively the subset and the random subset of the coordinates required for the evaluation of λi (•; ξ, θ) and λi (•; ξ, θ; U i ).

Sampling diffusion bridges
In order to employ the Zig-Zag sampler to simulate from the bridge measure we choose the truncation level N in equation (2).Then, under P This is a straightforward consequence of the change of measure in (12) and the Lévy-Ciesielski construction.
We need to make one further assumption: Assumption 4.2.The drift b of the diffusion process is twice differentiable.
Assumption 4.2 is necessary in order to compute the ξ k -partial derivative of the energy function, which becomes where As the index k in the Faber-Schauder basis function gets larger, both the magnitude of φ k and the size of its support decrease so that typically h k (s; ξ N )ds gets smaller and ∂ k ψ(ξ) ≈ ξ k which corresponds to the partial derivative of the energy function of a standardized Gaussian random variable with independent components.This justifies one more time the intuition that for high levels i, the random variables ξ ij , j = 1, ..., 2 i − 1 are approximately normally distributed and almost independent from the other random coefficients.
In order to avoid the evaluation of the integral appearing in (16) and the difficulty of drawing a Poisson time from its corresponding rate (6), we employ the subsampling technique.Considering ξ N nonrandom, we take as an unbiased estimator for where U k ∼ Unif(S k ) and as the bounding intensity rate where Φk = max s (φ k (s)) and The subsampling technique avoids the numerical computation of the time integral ( 16), thus avoiding a numerical bias and reducing the computational effort from O(T ) (for fixed discretization size) to O(1).The variance of this unbiased estimator can be reduced by averaging over multiple independent uniform draws or similar strategies (see for example Section 5.4), albeit at the cost of additional computations.In Section 5 we show specifically for each numerical experiment how we derived the Poisson upper bounds λi .
The compact support of the Faber-Schauder functions induce a sparse dependency structure on the target measure π.Indeed, X t only depends on those values of ξ l,k for which t ∈ S l,k .See Figure 4 for an illustration.It is easy to see that this implies that ∂ψ(ξ N ) ∂ξ (i,j) depends only on those ξ (k,l) for which the interior of S i,j ∩ S k,l is non-empty.In particular, 0 T S 3,0 S 3,1 S 3,2 S 3,3 S 3,4 S 3,5 S 3,6 S 3,7 Figure 4: Support of the Faber-Schauder functions (φ i,j : i ∈ {0, 1, . . ., N }, j = {0, 1, . . ., 2 i − 1} with N = 3.The coefficient ξ i,j is independent of the coefficient ξ k,l conditionally on the set of common ancestors (ξ m,n : define the relation ξ i,j ξ k,l to hold if S k,l ⊂ S i,j .If this happens, then we refer to ξ i,j as the ancestor of ξ k,l (and conversely ξ k,l as the descendant).Then the sets in Assumption 4.1 (using double indexing) can be chosen as where N is the truncation level.Formally, N i,j are the neighborhoods of the interval graph induced by ((S i,j : i ∈ {1, 2, . . ., N }, j ∈ {0, 1, . . ., 2 i −1})) with vertices {(i, j) : i ∈ {1, 2, . . ., N }, j ∈ {0, 1, . . ., 2 i − 1}}, where there is an edge between (i, j) and (l, k) if the interior of S i,j ∩S k,l is non-empty (see Figure 11).The factorization of the partial derivatives leads to a specific dependency structure of the coefficients under the target diffusion bridge measure: the coefficient ξ i,j is conditionally independent of the coefficient ξ k,l if S i,j ∩S k,l = ∅ conditionally on the set of common ancestors (ξ m,n : ξ m,n ξ i,j ∧ξ m,n ξ k,l ).This argument is made more formal by decomposing the likelihood function in Appendix A.

Numerical results
We show numerical results for three representative examples.In general, when applying our method, we start from a model (1), devise a representation of the approximate diffusion bridge (12), that we sample using generic implementations of algorithms 1-4 from our package, which are easily adapted to the task of sampling the coefficients of the Faber-Schauder expansion.To this end, we provide the k-th partial derivative of the energy function (16) or an upper bound to the Poisson rate (18) as argument for the sampler, as well as the sets N i,j as given in Section 4.4.The reader is referred to the file faberschauder.jl in the public repository https://github.com/SebaGraz/ZZDiffusionBridge/srcfor the implementation of the expansion and for the generic implementation of the different variants of the Zig-Zag sampler to our package (Schauer and Grazzi, 2020).
The first class of diffusion processes considered are diffusions with linear drift function (Subsection 5.1).This is a special case, where our method does not require the subsampling technique described in Subsection 4.1 and only Algorithm 3 has been employed.Notice that for this class, the transition kernel of the conditioned process is known.In Subsection 5.2, we apply our method for diffusions which substantially differ from Brownian motions, being highly non-linear and multimodal and therefore creating challenging bridge distributions for standard MCMC.Here we use the the fully local algorithm (Algorithm 4).In the specific example considered, the implementation of the Zig-Zag sampler is facilitated by the drift function and its derivatives being bounded and therefore a bounded Poisson rate for the subsampling technique is available.In view of this, we choose for the third numerical experiment a diffusion with unbounded drift (Subsection 5.3).For all the models, Assumptions 3.1, 3.4 and 4.2 are immediate to verify and Assumption 4.1 is satisfied.For each experiment, the burn-in τ burn-in and final clock τ final are manually tuned by inspecting the trace of ξ N and ensuring that the process reached stationarity before τ burn-in and fully explore the state space before the final clock τ final .The computations are performed with a conventional laptop with a 1.8GHz intel core i7-8550U processor and 8GB DDR4 RAM.We wrote the program in Julia 1.4.2 which allows profiling and optimizing the code for high performance.The program is publicly available on GitHub at https://github.com/SebaGraz/ZZDiffusionBridgewhere the reader can follow the documentation to reproduce the results.

Linear diffusions
A linear stochastic differential equation conditioned to hit a final point v T has the form for some (α, β) ∈ R 2 .Assumptions 3.1, 3.4 and 4.2 can be easily verified.In this case the energy function of the target distribution is for some constant C 1 , C 2 .Note that ψ is a quadratic function of ξ, which means that the target density is still Gaussian under P u,v T N .It follows that Interchanging the integral and the sum, this becomes where Φ k = φ k dt, Φ jk = φ k φ j dt, Φk = φφ k dt and Φ = φφ k dt.This is a linear function of ξ N and, for each i, the event times with rates λ i , see (6), can be directly simulated without upper bounds.Figure 5 shows samples from the resulting diffusion bridge measure with α = −5, β = −1 obtained with this method running the Zig-Zag sampler for τ final = 1000, with a burn-in time of τ burn-in = 10.The closed form of the expansion of linear processes, or more generally, reciprocal linear processes, with the Faber-Schauder basis was also found and used in Meulen, Schauer, and Waaij (2018) for the problem of nonparametric drift estimation of diffusion processes.The results are validated by computing analytically the density of the random variable X T /2 (which, for the linear case, is known in close form) and comparing this with its empirical density obtained from one sample of the Zig-Zag process (see Figure 7, left panel).

Non-linear multi-modal diffusions
The stochastic differential equation considered here has the form for some α ≥ 0. When α = 0 the process is a standard Brownian motion while for positive α, the process is attracted to its stable points (2k − 1)π, k ∈ N. Assumption 3.1, 3.4, 4.2 follow from drift, its primitive and its derivative being globally bounded.Fixing N , the energy function is given by Using trigonometric identities, we obtain that To avoid the need to find the roots of equation (7) we apply the subsampling technique described in Subsection 4.1.Since the drift and its derivatives are bounded, we can easily find the following upper bound for (14): with a 1 = Φk S k (α 2 + α)/2, Φk = max(φ k ) and ξ k (t) = ξ k + θ k t.In this case, the upper bound λi is a function only of the coefficient ξ i .Figure 6 shows the results obtained with this method setting α = 0.7.For this diffusion, the non-linearity and multiple modes make the mixing of the Zig-Zag sampler slower so we set τ final = 10000 and burn-in τ burn-in = 10.Analysing the goodness of the empirical diffusion bridge distribution obtained is a difficult task since the true conditional distribution is not known in a tractable form.We start by checking if some geometrical properties of the diffusion bridge distributions are preserved in the simulations.For example, in Figure 6, it can be noticed that the diffusion is attracted to the stable points ±π, ±3π, ..., and symmetric (geometrically speaking, after rotation) around the vertical axes t = T /2.We furthermore validate our method by simulating forward diffusion processes, using Euler discretization in a fine grid, and retaining only the paths which end in a -ball of a certain point at time T ( -ball forward simulation).If the final point is such that the probability of ending in this -ball is high enough, we can create in this way a sample from the approximated bridge and compare it to the samples obtained from the Zig-Zag.The right panel of Figure 7 shows the joint empirical distribution with the two methods of the first quarter and third quarter random variables.Finally, Figure 8 illustrates that the marginal distribution of the coefficients in higher levels is approximately Gaussian and the non-linearity of the process is absorbed by the coefficients in low levels.

Diffusions with unbounded drift
Here we consider stochastic exponential logistic models.For this class, the process grows exponentially with rate r until it reaches its saturation point K. Its dynamics are perturbed by noise which grows as the population grows.The resulting stochastic differential equation takes the form We can transform the process in order to get a new process with unitary diffusivity σ = 1 (Lamperti transform with X t = − log(Y t )/β).The transformed differential equation becomes with c 1 = β/2 − r/β and c 2 = r/(βK).Note that the drift function b of the transformed process is not global Lipschitz continuous.Nevertheless Assumptions 3.4 and 4.2 are satisfied and by Remark 3.9, also Assumption 3.1 is verified.In this case, the partial derivative of the energy function is given by where a 1 = 2r 2 /(βK), a 2 = a 1 /K.As before, it is not possible to simulate directly the first event time using the Poisson rates given by equation (6).The subsampling technique requires an upper bound for the unbiased estimator (14).Define the following quantities b with λ (1) k e 2β k t and z (1) k exp(−βb k , φi = max s φ i (s).
Using the superposition theorem (see for example Grimmett and Stirzaker, 2001), we can simulate a waiting time with Poisson rate ( 23) by means of simulating three waiting times according to the Poisson rates λ k and then take the minimum of the three realizations.Since at any time t > 0, either λ  k (t) is 0, we just need to evaluate two waiting times.Figure 9 shows the results obtained with our method for this process.The final clock of the Zig-Zag sampler is set to T = 1000 and initial burn-in time τ burn-in = 10.

Numerical comparisons
In this section we benchmark the fully local Zig-Zag sampler against the Metropolisadjusted Langevin algorithm (MALA) (Roberts and Tweedie, 1996), Hamiltonian Monte Carlo (HMC) (Duane et al., 1987) and another well known PDMP, the Bouncy particle sampler (Bouchard-Côté, Vollmer, and Doucet, 2015).The Bouncy Particle sampler can use the exact subsampling technique in a very similar way as explained in Subsection 4.1.According to the scaling limit results obtained in Bierkens, Kamatani, and Roberts (2018), the Zig-Zag is more efficient compared to the Bouncy Particle sampler in a high dimensional setting when the conditional dependency graph corresponding to the target measure exhibits sparsity (which clearly is the case here).The MALA sampler is a well known discrete-time Markov chain Monte Carlo method which performs informed updates through the gradient of the target distribution.HMC is considered a state-of-the-art algorithm.In contrast to PDMPs, for HMC and MALA the gradient needs to be fully evaluated and no subsampling methods can be exploited.Thus, the integral in (16) needs to be computed numerically, introducing bias.Furthermore, contrary to PDMPs, the resulting Markov chain is reversible.We study the performance of the samplers for the stochastic differential equation (20) with u, v = 0 and the time horizon T = 100 and we let α vary.As α increases, the target distribution on the coefficients presents higher peaks and valleys and is therefore a challenging distribution for general Markov chain Monte Carlo methods.We fix the refreshment rate of the Bouncy Particle sampler to 1 to avoid a degenerate behaviour and implement the MALA algorithm with adaptive step size over 250,000 iterations.We used the automatically tuned dynamic integration time HMC Algorithm (Betancourt, 2018) with 3, 000 iterations and with diagonal mass matrix and integrator step size both adaptively tuned in a warm-up phase of 2, 000 iterations, with the latter adapted using a dual-averaging algorithm (Hoffman and Gelman, 2014) with target acceptance statistic of 0.8.The algorithm is provided in the package AdvancedHMC.jl(Ge, Xu, and Ghahramani, 2018).The integral appearing in the gradient of the energy function is computed for the MALA sampler and for the HMC sampler numerically with a simple Euler integration scheme over 2 N +1 points, where N is the truncation level which is fixed to 6 for all the experiments.The final clock for the PDMPs is T = 25, 000.We also include the numerical results of two variants of the Zig-Zag sampler: (ZZv1) where the partial derivative in ( 16) is estimated by averaging over multiple independent realizations of (17), with the number of realizations is proportional to the length of the range of the integral in ( 16); (ZZv2) where the partial derivative in ( 16) is estimated by decomposing the range of the integral into N subintervals (with N proportional to the length of the range of the integral) and evaluating the integrand at a random point drawn inside each subinterval.
These variants of the Zig-Zag have been proposed after noticing that the coefficients at low levels are the ones deviating the most from normality and the partial derivative with respect to those coefficients have larger support.This suggests that refining the estimates of the partial derivative of the energy function only with respect to those coefficients can be beneficial and improve the performances of the PDMPs.Figure 10 shows the results obtained.The fully local Zig-Zag and its variants always outperform the Bouncy Particle sampler, the MALA and the HMC with respect to the statistics considered, namely the mean, median and minimum of the effective sample size computed for each coefficient of the Faber-Schauder expansion and the effective sample size of the coefficient ξ 0,0 , which gives the middle point X T /2 and, as shown in Figure 10, is one of the most difficult coefficients to sample.

Extensions
In this section we briefly sketch the extension of the approach presented in Section 3 to a class of multi-dimensional diffusion bridges.Then we study the scaling properties of the algorithm with respect to three quantities: the time horizon of the diffusion bridge T , the truncation level N and the dimensionality of the diffusion bridge d.The asymptotic variance estimate used for computing the ESS is obtained using batch means.Notice that, while the subsampling technique adopted for the piecewise deterministic Monte Carlo methods does not introduce bias on the target distribution, the numerical integration adopted for the MALA and HMC samplers introduces bias on the target distribution.

Multivariate diffusion bridge
Consider a d-dimensional diffusion bridge given the stochastic differential equation where (W t ) t≥0 is a d-dimensional Wiener process and ∇B : R d → R d is a conservative vector field, i.e. the gradient of some scalar-valued function B. Denote its law by P u,v T .Similarly to equation (10), under mild assumption on ∇B(X t ), we can write the change of measure between P u,v T and the standard d dimensional Wiener bridge measure Q where b = ∇B, ∆B is the Laplacian of B and C is a normalization constant which depends on u, v T and T .It is straightforward to derive an equivalent approximated measure as done in equation ( 12) and prove Theorem 3.11 in the multi-dimensional setting.In this case the d dimensional diffusion bridge measure is approximated by the same truncated expansion of equation ( 2) with coefficients ξ i,j , i = 0, ..., N ; j = 0, ..., 2 N which now are d dimensional random vectors.The total dimensionality of the target density for diffusion bridges becomes d(2 N +1 − 1) .Similarly to the one dimensional case, Proposition A.1 holds (the proof follows in a similar fashion of the proof of Proposition A.1 and is omitted for brevity).The Poisson rates λ k i,j (where, k ∈ {1, ..., d} defines the coordinate of the d dimensional process) are functions of the sets N k i,j which have maximum admissible size 1) so that Assumption 4.1 holds.

Scaling for large T, N, d
The following scaling analysis serves as preliminary work for future explorations.The expected run time of the fully local Zig-Zag sampler (Algorithm 4) is intimately related with the number of Poisson event times for a fixed final clock τ final and the conditional independence structure appearing in the target measure.The former is determined by the size of the Poisson bounding rates λ1 , ..., λM while the latter is defined by the sets N 1 , ..., N M and determine the complexity of the local step of Algorithm 3.
Proof.For every i = 0, 1, ..., N ; j = 0, 1, ..., 2 i − 1, the time horizon T and scaling index i enter in the bounding rates of (18) through the terms S i,j and φi,j .The first term is of O(T 2 −i ) and the second one is of O( √ T 2 −i/2 ).Proposition 6.1 helps understanding how the complexity of the algorithm scales as T grows and as the truncation level N grows.As T grows, the Poisson rates increase with order T 3/2 so that the total number of Poisson events for a fixed Zig-Zag clock increases with the same order.
Furthermore, as the truncation level N grows, the change of measure affects less and less the coefficients in high levels and the partial derivative of the energy function goes to zero with rate 2 −3N/2 ) implying that the for large N , λN,j ≈ C 2 = (ξ N,j θ N,j ) + (which is the Poisson rate for the Brownian bridge).As a consequence, the Poisson processes of the coefficients in high levels (i large) will be approximately independent with all the other coefficients and not function of the level i so that the complexity of Algorithm 4 scales approximately linearly with the number of mesh points.This is opposed to the standard Zig-Zag algorithm (Algorithm 1) which does not take advantage of the approximate independence of the coefficients in high levels so that the 2 N +1 − 1 waiting times have to be renovated at every reflection of each coefficient.
The scaling result under mesh refinement (when N grows) is unsatisfactory as the algorithm deteriorates when the resolution of the path increases.A partial solution can be obtained by letting the absolute value of the marginal velocities |θ N,j | decrease as N increases.This would enhance the scaling property of the algorithm under mesh refinement at the cost of a slow mixing of high level components.An alternative solution is considered in Bierkens et al. (2020) where the authors enhance the scaling property of the algorithm by replacing the Zig-Zag algorithm with the Factorised Boomerang sampler.The Factorised Boomerang sampler differs form the Zig-Zag by having curved trajectories which are invariant to a prescribed Gaussian measure.This allows the process to sample from the Gaussian measure (Brownian bridge measure) at barely no cost.However, the main drawback of the factorized Boomerang sampler is the current limiting techniques for simulating Poisson times given the curved trajectories which lead to Poisson upper bounds which are not tight.
Finally, when the dimensionality of the diffusion bridge is d 1, both the dimensionality of the target density of the Zig-Zag sampler and the sets N k i,j for i = 0, ..., N ; j = 0, ..., 2 i − 1; k = 1, ..., d grow linearly with d so that, in general, we expect the computational time to grow with rate d 2 .When the drift of the multidimensional bridge presents a sparse structure, i.e. not all coordinates of the differential equation interact directly with each other, as common in the high dimensional case arising from discretized stochastic partial differential equations (e.g.Mider, Schauer, and Meulen, 2020, Section 6), the size of those sets reduces considerably until the extreme case of d independent diffusion bridges where the sets N k i,j are not anymore a function of d and clearly the complexity grows linearly with the dimensionality d.

Conclusions
In this paper we have introduced a new method for the simulation of diffusion bridges which substantially differs from existing methods by using the Zig-Zag sampler and the basis of representation adopted.We motivated both choices and presented the method and its implementation.The resulting simulated bridge measures are shown to be close to the real measures, even for low dimensional approximations and bridges which are highly non-linear.We took advantage of the subsampling technique and a local version of the Zig-Zag to sample high dimensional approximation to conditional measures of diffusions with intractable transition densities.The subsampling technique is a key property in favour of using piecewise deterministic Monte Carlo methods for diffusion bridges (and whenever the target measure is expressed as an integral that requires numerical evaluation).However, the main limitation found for these methods is that they rely on upper bounds of the Poisson rates which are model-specific.Upper bounds for PDMC are easily derived in situations where the log-likelihood has a bounded Hessian.In our setting this means that we wish for the function b 2 (x) − b (x) to have bounded second derivative.In other cases, tailor made bounds need to be derived which can be substantially more complicated.Furthermore, the performance of these samplers can be affected if the upper bounds are too large.
In conclusion, this is the first time (to our knowledge) the Zig-Zag has been employed in a high dimensional practical setting.We claim that the promising results will open research toward applications of the Zig-Zag for high dimensional problems.We mention below some possible extensions of the methodology proposed which are left for future research: (a) The hierarchical structure of the Faber-Schauder basis suggests that the Zig-Zag should explore the space at different velocities to achieve optimal performance.Un-fortunately, it is not immediately clear how to tune the velocity vector; (b) In Section 6 we anticipated the possibility to simulate multidimensional diffusion bridges.In order to generalize the results presented in this paper, we assumed the drift being a conservative vector field.In order to relax this limiting assumption, new convergence results have to be derived which deal explicitly with the stochastic integral appearing in equation ( 8).
(c) The driving motivation for proposing this methodology is to perform parameter estimation of a discretely observed diffusion model.For this purpose, the Zig-Zag sampler runs jointly on the augmented path space given by the coefficients ξ and the parameter space Θ.
(b) In the same spirit as the local Bouncy Particle Sampler of Bouchard-Côté, Vollmer, and Doucet (2015) and Peters and With (2012), the local and the fully local Zig-Zag sampler introduced in Section 4 reduces the complexity of the algorithm improving its efficiency with respect to the standard Zig-Zag Algorithm as the dimensionality of the target distribution increases (see Subsection 6.2).This opens the way to high dimensional applications of the Zig-Zag sampler when the dependency graph of the target distribution is not fully connected and when using subsampling.The factorization of the log-likelihood and the local method we proposed is reminiscent of other work such as e.g.Faulkner et al. (2018), Michel, Tan, and Deng (2019) and Monmarché et al. (2020).

Figure 3 :
Figure 3: 100 samples from the Brownian bridge measure starting at 0 and hitting 0 at time 1 obtained by one run of the Zig-Zag sampler targeting the coefficients relative to the measure expanded with the Faber-Schauder basis.The resolution level is fixed to N = 6 and the Zig-Zag clock to τ final = 500 and initial burn in τ burn-in = 10.

Assumption 3 . 4 .Theorem 3 . 5 .
The drift b is continuously differentiable and b 2 + b is bounded from below.If Assumptions 3.1 and 3.4 are satisfied, then P u,v T N converges weakly to P u,v T as N → ∞.
Remark 3.9.If Assumptions 3.4 and 3.8 hold, then Assumption 3.1 is satisfied.Proof.By Subsection 3.5 in Karatzas and Shreve (1991), (Z t ) is a local martingale.Say b (x) + b 2 (x) ≥ −2C, where C ≥ 0. Using the assumptions, we have 2,...,k−1 end procedure 4.2 Local Zig-Zag sampler Subsection 3.1 of Bouchard-Côté, Vollmer, and Doucet (2015) proposes a local algorithm for the Bouncy Particle Sampler which is a process belonging to the class of piecewise deterministic Markov processes.Similar ideas apply to our setting.Assumption 4.1.The Poisson rate λ i for a d-dimensional target distribution is a function of the coordinates N i ⊂ {1, . . ., d},

Figure 6 :
Figure 6: Simulation of the diffusion bridge measure (200 samples) given by equation (20) with α = 0.7 starting at −π at time 0 and hitting 3π at T = 50.Truncation level N = 6, final clock τ final = 10000 and burn-in 10.The straight horizontal lines are the attraction points of the process.

Figure 7 :
Figure 7: On the left panel: comparison between empirical distribution (blue line, computed with a kernel estimator) and the exact distribution (red line) of the mid-point random variable X 5 for the linear diffusion (equation 19) with a = −5 and b = −1.The empirical distribution has been extracted from the same experiment shown in Figure 5. On the right panel: comparison between the joint distribution of the variables X T /4 and X 3T /4 of the process given in equation (20) starting at −π and hitting π at T = 50.The scatter plot with red dots are obtained with -ball Euler simulation with = 0.1 and discretization ∆t = 0.0005 while the blue continuous path is the Zig-Zag path.

Figure 8 :
Figure 8: Q-Q (quantile-quantile) plot against standard normal distributions of the sample path of 7 coefficients respectively at level 0, 1, 2, 3, 4, 5, 6 targeting the conditional bridge measure given by equation (20) with α = 0.7 and initial point u = 0 and final point v = 0 at T = 100.On the bottom right panel, the heatmap of the absolute value of the sample correlation between the coefficients at different levels.The blue straight lines correspond to the marginal measures of the coefficients relatively to a Brownian bridge.

Figure 9 :
Figure 9: Simulation of the diffusion bridge measure (100 samples) given by the logistic growth model (equation (22)) with parameters K = 2000, r = 0.08, β = 0.1, starting at the value 50 and hitting 1000 at time 200.Truncation level N = 6, final clock τ final = 1000 and burn-in τ burn-in = 10 .The blue smooth line is the solution of the deterministic logistic model without final condition.

Figure 10 :
Figure 10: Performance comparison of the fully local Zig-ZaZ (ZZ), its variants (ZZv1 and ZZv2), the Bouncy Particle sampler with Subsampling (refreshment rate set to 1), MALA and HMC sampler.The performance measure considered here are respectively the effective sample size (ESS) of the middle point X T /2 , the median and the minimum of the ESS over the dimension of the coefficients of the expansion.The target diffusion bridge with drift b(x) = α sin(x) with u, v = 0 and T = 50 and truncation level N = 6.The final clock for the PDMPs is set to T = 25000, the number of iterations for the MALA is set to be 250000 with adaptive time step targeting the acceptance rate 0.6 (Roberts and Rosenthal, 1998) and the number of iteration for the HMC is 3000 with the algorithm finetuned by the package AdvancedHMC.jl.All the quantities are normalized by the runtime of execution.The asymptotic variance estimate used for computing the ESS is obtained using batch means.Notice that, while the subsampling technique adopted for the piecewise deterministic Monte Carlo methods does not introduce bias on the target distribution, the numerical integration adopted for the MALA and HMC samplers introduces bias on the target distribution.