{'fraction_non_alphanumeric': 0.04595953978670028, 'fraction_numerical': 0.01427353279205131, 'mean_word_length': 3.8717639975918123, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 51, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adaptive schemes for piecewise deterministic Monte Carlo algorithms

The Bouncy Particle sampler (BPS) and the Zig-Zag sampler (ZZS) are continuous time, non-reversible Monte Carlo methods based on piecewise deterministic Markov processes. Experiments show that the speed of convergence of these samplers can be affected by the shape of the target distribution, as for instance in the case of anisotropic targets. We propose an adaptive scheme that iteratively learns all or part of the covariance matrix of the target and takes advantage of the obtained information to modify the underlying process with the aim of increasing the speed of convergence. Moreover, we define an adaptive scheme that automatically tunes the refreshment rate of the BPS or ZZS. We prove ergodicity and a law of large numbers for all the proposed adaptive algorithms. Finally, we show the benefits of the adaptive samplers with several numerical simulations.


Introduction
Piecewise Deterministic Markov processes (PDMP) have been recently used for Monte Carlo sampling as a continuous time alternative to Markov chain Monte Carlo (MCMC) methods. The Bouncy Particle sampler (BPS) and the Zig-Zag sampler (ZZS), introduced in [13] and [8] respectively, are primary examples of this new class of methods, with notable predecessors [31,28,29]. Both samplers are based on continuous time dynamics defined by piecewise linear, deterministic trajectories and changes in the velocity of the process at random times. The simplicity of these dynamics is such that the underlying processes can be simulated exactly in many relevant scenarios. A key aspect of these samplers is the non-reversibility of the underlying Markov process. This property has been observed to result in a lower asymptotic variance of Monte Carlo estimates for moments of the target density [2,7]. Moreover, these algorithms can be naturally modified to achieve asymptotically exact subsampling. In the Bayesian setting, this means that Piecewise Deterministic Monte Carlo (PDMC) algorithms need to access only a small portion of the data set at every iteration without introducing a bias. The reader is referred to [20,39] for an in depth description of the general methodology of PDMC algorithms.
Several papers have studied the convergence properties of the BPS and the ZZS in recent years. It was first observed in [13] that the BPS can fail to be ergodic unless a refreshment of the velocity vector is performed at random times. Exponential ergodicity of the BPS was proved in [17,18] for different distributions from which refreshments of the velocity vector can be drawn. Similarly, in [12] the ZZS was shown to converge exponentially to its invariant distribution under reasonable assumptions and without the need of velocity refreshments. Exponential convergence in the L 2 sense was established for both samplers in [1] using the hypocoercivity framework, and recently, using Poincaré inequalities in space-time, in [26]. A study of the scaling limits was conducted in [10], giving also a criterion to choose the refreshment rate of BPS. The asymptotic variance of these processes has been studied in e.g. [2,7]. It is also possible to design PDMC algorithms with non-linear trajectories, see e.g. [39,9].
Although PDMC sampling methods offer some important benefits as mentioned above, computation remains expensive, which requires us to investigate possible performance improvements. In particular, a strong performance degradation is observed when the target distribution π is anisotropic. Figure 1 illustrates this phenomenon in the case of Gaussian targets as a function of the correlation between all components. The performance drop occurs due to a combination of decreasing accuracy of the estimates and increasing computational complexity of the algorithms, which is implied by the growing number of velocity change events. Our idea to improve this issue is to let the process learn (part of) the covariance matrix Σ π and take advantage of it to enhance the mixing properties. The covariance estimate is used to linearly transform the target in such a way that it becomes more isotropic, i.e. with unitary covariance matrix. The standard samplers are then run targeting the transformed version of π, and the obtained sample is finally re-transformed to be approximately from π. The procedure is applied iteratively, and once a new estimate of Σ π is computed, it is used to define the linear transformation of π. The estimate will eventually be close to the true covariance matrix and the process targets an isotropic version of π. This scheme can also be interpreted as an application of a linear transformation directly to the standard ZZS and BPS. The natural applications of this procedure are then targets with elliptical level curves, although performance improvements can be observed also for distributions that deviate from this class.
Furthermore, we address the problem of automatically choosing the rate of velocity refreshments. In [10] the authors consider a BPS with a specific target and derive that in the limit it is optimal to have a ratio of number of refreshments over number of total events of 0.7812. In this paper we use this criterion as a basis to define an adaptive algorithm that iteratively adjusts the refreshment rate to obtain the right ratio. This same adaptive scheme can be applied to the ZZS. Indeed it could be the case that adding velocity refreshments to the ZZ process leads to a faster convergence to the invariant measure, although it comes with a larger asymptotic variance. For an analysis of these two results we refer to [40] and [7] respectively.
Both the schemes we discussed take advantage of what the process has learned up until the current time to tune a parameter or to improve the performance. This idea is at the core of adaptive Markov chain Monte Carlo algorithms. For an introduction to this area we refer to [3,35], while standard results on convergence of these methods can be found among others in [5,21,22,34]. It is well known that adaptive MCMC algorithms can lose the right invariant measure if not applied with care (see for instance [34] and the examples therein). Therefore we study in depth the convergence properties of the proposed algorithms. To fit into the existing adaptive MCMC literature we let the adaptation happen at fixed points in time. The main challenge consists of establishing a simultaneous geometric drift condition for a family of BPS's (see Lemma 4.5) and a simultaneous small set condition for a family of ZZ processes (see Lemma 4.1). The former result is obtained taking advantage of the Lyapunov function found in [18], while the latter is proved by extending the one-dimensional case considered in Lemma B.2. The ergodicity and a law of large numbers for the proposed adaptive PDMC algorithms are then established in Theorem 3.8 and Theorem 3. 10.
In Section 2 we introduce the adaptive schemes, while in Section 3 the theoretical aspects of the algorithms are studied. The proofs of the two main theorems can be found in Section 4. In Section 5 the adaptive BPS and ZZS are tested empirically on various Gaussian targets, on a Bayesian logistic regression problem with correlated data, and on a mixture of two Gaussian distributions. Appendix A gives details on the implementation of adaptive PDMC (with and without subsampling), as well as an alternative adaptive scheme for the refreshment rate. Appendix B contains the proofs of other theorems stated in Section 3, together with more technical results.

The adaptive schemes
We are interested in building adaptive strategies to make the ZZS and the BPS choose the refreshment rate themselves and/or converge faster to the target density. We begin with an introduction of the standard versions of both samplers, followed by a characterisation of the preconditioned processes in Section 2.2 and a discussion on the choice of the transformation matrix in Section 2.3. Finally, the adaptive algorithms are defined in Section 2.4.

The standard ZZS and BPS
Let the target density π be defined on X ⊂ R d as where U (ξ) is called potential or energy function, and ξ ∈ X . Let us now define the standard ZZS with invariant measure π. Throughout the paper the position and velocity vectors at time t of the standard processes, both ZZS and BPS, are denoted respectively as Ξ(t) and Θ(t). We distinguish between the Zig-Zag Sampler (i.e. a PDMC algorithm) and the Zig-Zag Process (i.e. the Markov process on which the algorithm is based). The Zig-Zag (ZZ) process is a Markov process (Ξ(t), Θ(t)) t≥0 with state space E = R d × {−1, +1} d that follows a linear trajectory in the position space with velocity Θ until one of its d-inhomogeneous Poisson clocks rings. When the i-th clock rings the velocity of the i-th coordinate switches sign. In mathematical terms, this means that the new velocity vector is F i Θ, where F i : {−1, +1} d → {−1, +1} d is the flip operator defined as The rates of the Poisson clocks are defined for i = 1, . . . , d as where (θ i ∂ i U (ξ)) + = max(0, θ i ∂ i U (ξ)) and γ i (ξ, θ) : E → R + is called excess switching rate and is such that γ i (ξ, θ) = γ i (ξ, F i θ). It was shown in [8,Theorem 2.2] that this choice of switching rates ensures that µ = π × Unif({−1, +1} d ) is the stationary distribution of the process. Moreover, the ZZ process is characterised by its infinitesimal generator where function f should be in the domain of the generator D(L). The linear trajectories are represented in the first term, while the second term represents the event of a velocity flip. Similarly, we denote the BPS as (Ξ(t), Θ(t)) t≥0 , but in this case the state space is E = R d × R d and so Θ(t) ∈ R d . In contrast to the ZZS, the BPS has two Poisson clocks. The first one depends on the gradient of the energy function and has inhomogeneous rate λ(ξ, θ) = θ, ∇U (ξ) + = max(0, θ, ∇U (ξ) ). At event time the particle is reflected on the level curve of the potential U following an elastic bounce, and thus preserving the norm of the velocity vector. After a bounce the new velocity vector is given by the inverse transformation. For this reason it is important that the matrix M is invertible, and thus that we can go from one state space to the other. This procedure is illustrated in Figure 2. An equivalent option is to simulate directly the process (X t , Θ t ) t≥0 that results from the scheme in Figure 2. We will conveniently alternate between these two formulations when studying the ergodic properties of the samplers, while we will use the latter formulation for our experiments. The dynamics of process (X t , Θ t ) t≥0 are studied in the remainder of this section. Let us first focus on the case in which (Ξ t , Θ t ) t≥0 is a standard ZZS with excess switching rate γ. In this case the switching rates areλ M,i (ξ, θ) = (θ i ∂ iŨM (ξ)) + + γ i (ξ, θ) for i = 1, . . . , d. Note that, unless stated otherwise, we will always use a tilde to indicate quantities related to the standard PDMC samplers with transformed target. In the next proposition we find the generator of the preconditioned ZZS. For a characterisation of the domain of the extended generator we refer to [16,Theorem 26.14].
Proposition 2.1. Let M ∈ R d×d be an invertible matrix. Let (Ξ(t), Θ(t)) t≥0 be a standard ZZS with targetπ M and excess switching rates γ : E → R d + . The process (X(t), Θ(t)) t≥0 = (M Ξ(t), Θ(t)) t≥0 has extended generator (L M , D(L M )) where for any h ∈ D(L M ) in which for i = 1, . . . , d where M i denotes the i-th column of M .
Proof. All the proofs of this section are postponed to Appendix B.1.
Proposition 2.1 shows that the transformed process is again a PDMP with linear trajectories between jumps. The transformation affects the velocity of the process, which is now v = M θ, and the switching intensities, which are as defined in (3). In particular, the available velocities for a fixed transformation M are in the following set: When a switch of the i-th velocity of the underlying standard ZZ process happens, the velocities of the transformed process change according to the Therefore, all components of the velocity are possibly affected by any single event. If M is diagonal the behaviour is more similar to the standard ZZ process and in particular we have thatF i ≡ F i . In the proposition below we check that (X(t), Θ(t)) t≥0 targets the correct density function.
Let us now apply the same transformation scheme shown in Figure 2 to the BPS. In this case the switching rate of the standard BPS with targetπ M isλ M (ξ, θ) = θ, ∇Ũ M (ξ) + , while reflections on the level curves ofŨ M are obtained by applying operatorR M . The following result, analogous to Proposition 2.1, studies the transformed process.
where we defined and Once again the true velocity of the process is v = M θ. When a velocity refreshment takes place the Observe also that the reflection rate is λ M (x, θ) = v, ∇U (x) + and thus preserves the same structure as in the standard BPS. It follows that the complexity of the simulation of event times remains unchanged. Finally consider the reflection operator in (6). This corresponds to a reflection in the opposite direction to the gradient in the transformed space, i.e. ∇ ξŨ (ξ) = M T ∇U (x). After the bounce the process moves with velocity Proposition 2.4. Consider the same setting of Proposition 2.3. Then, for any invertible M ∈ R d×d , the transformed BPS (X(t), Θ(t)) t≥0 has invariant distribution µ = π × ψ.

Choosing the transformation matrix
As explained above, we wish to transform the target to mitigate its anisotropies. To this end, some alternative choices of the transformation matrix M are the following: a) M = Cov π (X): this transformation is such that the targetπ M has unitary covariance matrix.
The downside of this choice is the additional O(d 3 ) computations that are introduced by the calculation of the square root of the covariance; b) M is a rotation matrix (det M = 1) such that the transformed density has a certain angle. Although an interesting case, it is not clear whether there is an optimal angle that speeds up the convergence; c) M is the diagonal matrix with M ii = Var π (X i ) for any 1 ≤ i ≤ d. This choice introduces O(d) computations due to the square root of the variances, which is a negligible additional computational burden. However, correlations in the target are not picked up and only a rescaling of the axes is performed. The main advantage of this choice is that the scenario in which some components are explored quickly and others slowly is avoided.
Both the first and the third option can potentially change the expected number of switching events, and this could be an inconvenience in certain settings. It is not difficult to modify these transformations in such a way that the expected switching rate is enforced to be close to that of the original standard PDMC algorithm. For example, consider the transformed BPS with generator as in Proposition 2.3. Then in stationarity we have Since the standard case corresponds to M = I d , we can normalise any M by dividing it by its Frobenius norm. Then the upper bound is the same for all such choices of M and the expected switching intensity will be close to the standard case. This does not make a difference from a computational point of view as it just amounts to reparametrization of the time parameter, but prevents unpredictable behaviour of the algorithm. Naturally the options above are not available in practice as the covariance matrix is unknown. It is the goal of the next section to propose an adaptive scheme that overcomes this issue.

Adaptive PDMC algorithms
In the previous sections we defined the transformation scheme and we discussed the effect it has on the underlying process, together with different choices of the preconditioning matrix. We now describe how this idea can be applied in practice by designing an adaptive PDMC algorithm. Our general strategy is to simulate the process in continuous time and store the states of the process at discrete times. Then at predefined times the stored states are used to update the adaptation parameters. In addition to the adaptive preconditioner, we incorporate an adaptation of the refreshment rate, which makes its choice automatic.
Let us then define a family of Markov semigroups by P := {(P t Γ ) t≥0 : Γ ∈ Y}, in which Γ is the adaptation parameter, Y is a compact space, and (P t Γ ) t≥0 is the semigroup of a modified ZZS or BPS. The modification is given by the adaptation parameter, which is then Γ = (M, λ r ) for BPS and Γ = (M, γ) for ZZS. Thus Y is a suitable compact space of preconditioners and refreshment rates/excess switching rates. Naturally, it is also possible to choose Γ = M or Γ = λ r only. Now that we have defined a family of Markov processes, we define a rule that establishes how to choose a P Γ ∈ P at every iteration. Let us begin by introducing a discretisation step ∆t, which defines a discretisation of Algorithm 1 Adaptive PDMC sampler 1: Input: family of kernels P = {P Γ : Γ ∈ Y} 2: Input: initial condition (x, θ) ∈ E, Γ 0 ∈ Y, set B, ∆t, {p n } n≥0 , number of steps N 3: Output: sequence of discrete samples {X n , Θ n } N n=0 4: Initialise n = 0, (X 0 , Θ 0 ) = (x, θ), Q 0 = Γ 0 5: while n ≤ N do 6: With probability p n , set Γ n+1 = Q n+1 10: With probability 1 − p n , set Γ n+1 = Γ n 11: end if 12: n = n + 1 13: end while the time variable. At each time step n ∈ N, which corresponds to continuous time t = n∆t, the adaptive scheme can update the parameter Γ n based on the new information available, that is the new state of the process (X n , Θ n ). This defines a random sequence {Γ n } n≥0 . Once Γ n is computed, the next state of the process is given by (X n+1 , Θ n+1 ) ∼ P ∆t Γn ((X n , Θ n ), ·). Then one updates the parameter, obtains the next state of the process, and so on.
The definition above defines the core ideas, which are written in pseudo-code form in Algorithm 1. A few issues remain to be clarified. A first question is how to simulate the PDMP semigroup of either ZZS or BPS. Details on how the processes can be simulated in the case of a target with dominated Hessian can be found in Appendix A.1. In a big data setting (large number of observations, moderate dimensionality of the problem) it can be beneficial to take advantage of subsampling techniques that can be implemented with PDMC algorithms. In Appendix A.2 details can be found on how to make use of subsampling in the context of the adaptive schemes here discussed. For further information on the general implementation of ZZS and BPS we refer to [8,13]. A second aspect of Algorithm 1 we focus on is the introduction of set B, and thus of the auxiliary sequence of random variables (Q n ) n≥0 . The idea is to update the adaptation parameter Γ n only if (X n , Θ n ) ∈ B. This is useful from a theoretical point of view as it ensures that the process remains bounded in probability. Note that set B is defined by the user and can be chosen large. The auxiliary variable Q n is updated even if the process is outside of B and then, as soon as the process enters B, or if it was already in B, we let Γ n = Q n . A third characteristic of Algorithm 1 is the sequence {p n } n≥0 . This is a sequence for which p n ∈ [0, 1] for all n ∈ N and p n → 0 as n → ∞. The meaning is that at time step n we update the parameter Γ n with probability p n (assuming (X n , Θ n ) ∈ B), and with remaining probability (1 − p n ) we set Γ n = Γ n−1 . This choice is helpful when proving ergodicity of the adaptive scheme, as it enforces that the quantity of adaptation diminishes and eventually vanishes.
The function update(Q n , (X n+1 , Θ n+1 )) outputs the updated parameter given the new observation (X n+1 , Θ n+1 ). As suggested in [3], the estimation of the covariance matrix can be done sequentially, or online, by applyingμ Here {r n } n≥0 is a positive, decreasing sequence such that r n → 0 as n → ∞. In our simulations we choose r n = 1/n. Equation (7) is then used to define M n+1 such thatΣ n+1 = M T n+1 M n+1 . The same principle can be used if one is not interested in estimating the full covariance matrix, but only the diagonal, or more generally only a subset of it. More advanced estimation techniques can be employed to preserve any existing conditional independence structure in the target, as discussed in [41]. We remark thatΣ n+1 needs to be positive definite in order for M n+1 to be invertible as required. This property is achieved by choosingΣ 0 = 1 d×d , i.e. the identity matrix, and then observing that the second equation in (7) can be reformulated aŝ Indeed (X n+1 −μ n )(X n+1 −μ n ) T is non-negative definite and by inductionΣ n is positive definite, and thereforeΣ n+1 is itself positive definite. Moreover, in Section 3 we will see that to show ergodicity of the adaptive algorithms it is required that M n lies in a compact space of positive definite matrices.
Observe that positive definiteness follows from the fact that M n is the square root of a positive definite matrix, while it is sufficient to set bounds on the norm of M n in order to force it to be in a compact space. In particular we can impose that M n is not updated if the norm of the new estimate is outside of a user chosen interval [M min , M max ]. As M min and M max can be chosen arbitrarily small and large respectively, this condition is not restrictive in practice, although the choice of the cut-off value may influence convergence properties of the algorithm. Then refreshment rate of the BPS is assumed to be constant and is updated iteratively as follows. At time step n, n refl reflections took place and thus we estimate the average reflection rate as λ refl (n) = n refl /(n∆t). Therefore, using the optimality criterion in [10] we have An alternative adaptive strategy for the refreshment can be found in Appendix A.3. The scheme above can be applied to the excess switching rate of the ZZS. Although the analysis in [7] suggests that the best choice in terms of asymptotic variance is γ ≡ 0, adding some diffusivity could speed up the convergence to the target measure. In practice the user can select the wanted ratio of velocity switches over total events and proceed as above. However, a criterion to choose this ratio is currently unavailable for ZZS, and thus in this paper we limit ourselves to a theoretical study of this option. Finally, we remark that in practice it is not reasonable to update the parameters at every iteration. The main reason for this is the computational cost of such an operation. In the most general case, the task of learning all components of Σ takes O(d 2 ) operations, while the computation of its square root, which is needed to obtain the transformation matrix M , is an O(d 3 ) operation. Therefore it is rather inconvenient to perform this at every time step. To avoid this issue it is sufficient to define the adaptive scheme such that adaptations happen every n adap time steps, where n adap is a user-defined integer. A possible choice is for instance n adap = 1000. This modification is beneficial also because it allows the process to explore the target distribution before updating the parameters. Similarly, it is reasonable to update the refreshment rate based on the previous n adap time steps, as in the long term this allows to stabilise around the wanted ratio. The covariance matrix can be updated as in (7) by simply processing the entire batch of n adap data points one at a time.

Theoretical results
In the context of adaptive MCMC algorithms, convergence to the target density is usually proved with simultaneous drift conditions and small set conditions. In Section 3.1 we introduce the notation and the main existing theorems we make use of, and we extend these results to more general conditions in Theorem 3.5. In Section 3.2 we state Theorem 3.6, which shows ergodicity for an adaptive MCMC algorithm based on a continuous time process. In this result, the assumptions are formulated directly in continuous time. Finally, Theorems 3.8 and 3.10 in Section 3.3 show that the adaptive ZZS and the adaptive BPS discussed in Section 2.4 are ergodic and satisfy a weak law of large numbers under reasonable growth conditions on the potential.

Theory of adaptive MCMC
We denote the parameter that specifies the kernel as Γ ∈ Y. At time step n a Y-valued random variable Γ n determines which transition kernel will be used to move to the next step. From here on each Markov transition kernel P Γ is assumed to define a Markov chain that has µ as stationary measure, and moreover it is aperiodic and irreducible. An adaptive MCMC algorithm is then said to be ergodic if where ||·|| TV is the total variation distance, i.e. µ− ν TV = sup A⊂E |µ(A)− ν(A)|. A crucial quantity turns out to be the ε-convergence time function M ε : E × Y → N, defined as The next theorem, proved in [34], is arguably the most important result for establishing ergodicity of adaptive MCMC methods. (a) (Containment condition) For all z 0 ∈ E, Γ 0 ∈ Y, ε > 0 the sequence {M ε (Z n , Γ n )} ∞ n=0 is bounded in probability given z 0 , Γ 0 ; (b) (Diminishing adaptation) The following limit holds in probability: The boundedness of {M ε (Z n , Γ n )} ∞ n=0 can be rephrased as for all z 0 ∈ E, Γ 0 ∈ Y, δ > 0, there exists N ∈ N such that P (M ε (Z n , Γ n ) ≤ N |z 0 , Γ 0 ) ≥ 1 − δ, for all n ∈ N.
We are then interested in sufficient conditions that imply containment. A first case is the following, and was studied in [5].
(a) C is a uniform (ν Γ , δ, n 0 )-small set, i.e. for each Γ, there exists a probability measure ν Γ (·) on C such that P n 0 Γ (z, ·) ≥ δν Γ (·) for all z ∈ C; Then [5, Theorem 3] establishes that an SGE family satisfies the containment condition. In Section 3.3 we use this result to show that containment holds for the adaptive ZZS when the class of preconditioners is restricted to diagonal matrices.
In practice it is often hard to show that the family of Markov kernels is SGE, as it is not trivial to find a Lyapunov function that satisfies the simultaneous geometric drift condition. In [15] the authors introduced a way around this problem, although in a different context, and in [14] this was applied to adaptive MCMC. The fundamental idea is that it is possible to weaken the simultaneous drift condition by allowing adaptations only at time steps n at which the process Z n is inside of a compact set B. This means that, defining an auxiliary random process {Q n } n≥1 that contains the current adaptation parameter independently of the position of Z n , Γ n is updated as This modification avoids unbounded detours of the process by sticking to the same ergodic kernel once the process exits a fixed compact set. The compact set can be chosen arbitrarily large, and therefore in most applications the process will not exit from it. With this is mind we introduce the following sets of assumptions, which we show in Theorem 3.5 to be sufficient to enforce the containment condition.

c) the adaptation parameter is allowed to be updated only if the process is inside of a compact set
B, as defined in (11). There exist α, λ ∈ (0, 1), (c) the adaptation parameter is allowed to be updated only if the process is inside of a compact set B, as defined in (11). The proof can be found in Appendix B.2.

Remark.
A weak law of large numbers (WLLN) for bounded and measurable functions follows immediately from containment and diminishing adaptation by Theorem 3.4 in [32]. Therefore under the conditions of Theorem 3.5 a WLLN holds.

Convergence properties of adaptive MCMC algorithms based on continuous time Markov processes
It could be the case, as it is in the present work, that one is interested in defining an adaptive scheme based on a family of continuous time Markov processes in continuous time. In this case a grid for the time variable needs to be introduced in order to indicate the times at which the adaptation occurs. In fact the adaptive chain only sees the process at times m∆t, where ∆t > 0 is the step size and m ∈ N.
Although the resulting chain is in discrete time, it is in most cases easier to work directly with the continuous time process. The following result, which is analogous to Theorem 3.5, is helpful in this sense.
Theorem 3.6. Consider a family of Markov processes with generators {L Γ : Γ ∈ Y}, each being irreducible, aperiodic, and having µ as invariant measure. Consider a grid for the time variable with step ∆t. Consider an adaptive scheme that at times m∆t, with m ∈ N, chooses a process from the aforementioned family. Furthermore, suppose that the adaptive scheme satisfies the diminishing adaptation condition (10) for P := P ∆t . Finally assume one of the following two sets of conditions holds: (c) it holds that ∆t = t 0 .
If the adaptation parameter is allowed to be updated only if the process is inside of a compact set B, as defined in (11), then the adaptive algorithm satisfies the containment condition and is thus ergodic.
The proof of this theorem can be found in Appendix B.3.
Remark. The restrictions on the step size can be milder than as stated in Theorem 3.6. For instance, if the minorisation condition (1a) of Theorem 3.6 holds for all t ≥ t 0 , then one is free to choose any step size ∆t > 0. Furthermore, in both cases the assumption that the parameter can be updated only if the process is inside of a compact set at the adaptation time can be dropped when a simultaneous geometric drift condition holds (Assumption 3.2(b)).

Convergence properties of adaptive PDMC algorithms
Relying on Theorem 3.6, in this section we turn our attention to the ergodicity of adaptive PDMC algorithms. The proofs of the two theorems are postponed to Section 4. First, let us consider the adaptive ZZS. We assume the following conditions on the potential.
Let us now state the ergodicity result for the adaptive ZZS.
Let P = {P M,γ : M ∈ M, γ ∈ Λ} be a family of preconditioned Zig-Zag processes with generators defined by Equation (2). Suppose Assumption 3.7 holds and assume either of the following conditions holds:

has no additional restrictions but adaptations are allowed only inside of a compact set B,
and let ∆t be the discretisation step. Then the containment condition holds. Moreover, if the adaptive strategy is as described in Section 2.4 and is such that p n → 0 as n → ∞, then the diminishing adaptation condition holds and thus for µ = π × Unif(Θ): Finally, for any bounded and measurable f : Remark. The time discretisation step ∆t can be chosen freely and is not subject to constraints. Moreover, we remark that under condition (a) on M the adaptive algorithm is SGE, i.e. satisfies Assumption 3.2, while under condition (b) it satisfies the first set of conditions in Theorem 3.6. Thus if one is interested in learning only the diagonal elements of the covariance, then it is possible to take B = R d and allow adaptations independently of the state of the process.
Remark. It was shown in [12] that the ZZS is geometrically ergodic under Assumption 3.7 also in the case γ = 0, whereas in Theorem 3.8 we require γ(x, θ) ≥ γ min > 0. This extra assumption is convenient when proving a simultaneous small set condition (see Lemma 4.1). Based on similar arguments as in [12], we expect the statement of Theorem 3.8 to remain valid even in the case γ min = 0. In practice, one is free to choose γ min very small and thus this assumption does not represent a severe limitation.
Below we introduce a set of assumptions that is used to show ergodicity of the adaptive BPS. Here we limit our attention to the case of ψ = N (0, ½ d ).
Assumption 3.9 (Assumptions A1, A2, and A7 in [18]). Let U :  (4), where λ r is the refreshment rate. Suppose Assumption 3.9 holds and let ∆t be the discretisation step. Assume that ψ = N (0, ½ d ). If adaptations are allowed only inside of a compact set as explained in Equation (11), then the containment condition holds. Furthermore, for the strategy discussed in Section 2.4 the diminishing adaptation holds as long as p n → 0 as n → ∞.

family of preconditioned BPS's as defined in Equation
Thus the ABPS is ergodic in the sense of Equation (12) and satisfies a WLLN of the form (13) for any bounded and measurable f : Remark. Proving ergodicity of the adaptive BPS with refreshments from S d−1 , i.e. the unit sphere centred at the origin, is more challenging due to the more involved drift condition proved in [17].
In particular, it is not straightforward to convert it into a simultaneous drift condition as required in assumption 2(b) of Theorem 3.6.
4. Proofs of the main theorems 4.1. Proof of Theorem 3.8 In order to prove the theorem we show that, for suitable families of preconditioners, either Assumption 3.2 holds for the family of discretised ZZ processes, or condition (1) in Theorem 3.6 holds for the family of continuous time ZZ processes. In the next two sections we show auxiliary results, while in Section 4.1.3 we assemble them to show the theorem.

Minorisation condition for the ZZS
The following lemma shows that a simultaneous small set condition holds for the family of ZZ processes. The strategy of the proof is to reduce the d-dimensional minorisation condition to 1-dimensional conditions for every component of the process. Then we can take advantage of Lemma B.2 in Appendix B.4.1, which establishes that a simultaneous minorisation condition holds for a 1-dimensional ZZ process as long as lower and upper bounds for the switching rates are available.
Then for any set of the form In particular, t 0 and δ do not depend neither on M nor on γ.
Proof. We first show that all compact sets are uniformly small for the family of standard Zig-Zag processes with targets in Let M ∈ M and γ ∈ Λ and denote as (Ξ M,γ (t), Θ M,γ (t)) t≥0 the ZZ process that targetsπ M (ξ) and has switching rateγ M (ξ, θ) = γ(M ξ, θ). LetP M,γ be the corresponding semigroup. Then observe that for any rectangle Then observe that, independently of the values at any time s < t of the other components of the process, where we have defined Here λ max is well defined because U ∈ C 1 . In particular neither of the bounds in (15) depend on the specificπ M and thus hold for any M ∈ M. Therefore, we can apply Lemma B.2 with t 0 = 2R + ε with ε > 0 to each component of the product (14). It then follows that for any t ≥ t 0 there exists δ 1 > 0 that satisfies where Most importantly, δ depends only on the bounds on the switching rates, which are uniform for all γ ∈ Λ. Therefore Equation (16) holds for the same δ and ν d for all M ∈ M, γ ∈ Λ, and any rectangle B ⊂ C.
Since the set of rectangles is a π-system and generates the Borel σ-algebra, by the monotone class theorem (Theorem 6.2 in [24]) it follows that the lower bound of Equation (16) holds for any Borel set B. Since any compact set can be included in a large enough hypercube, for all sets C = D × V , with D compact and V ⊆ Θ, there are ν d and t 0 > 0 such that for all t ≥ t 0 there exists δ d > 0 for which C is uniformly (t, δ d , ν d )-small for the family of standard ZZ processes defined above. Now we wish to translate this result to the family of linearly transformed Zig-Zag processes with excess switching rates in Λ and target π. Denote as (X M,γ (t), Θ M,γ (t)) t≥0 the process with generator L M,γ and observe that for any t > 0 and any set : (x, θ) ∈ C, M ∈ M} which is itself a compact set and depends on M, but not on the specific M ∈ M. Then, for all t ≥ t 0 , with t 0 large enough, it holds that and δ is chosen such that C M is a uniform (t, δ, ν)-small set for the family of standard ZZ processes with targetsπ M and switching ratesγ M . Therefore there is no dependence neither on M nor on γ in t 0 and δ and the proof is concluded.

Drift conditions for the Zig-Zag process
If we restrict our attention to the class of diagonal matrices with positive, bounded entries, then the Lyapunov function in [12,Lemma 11] satisfies also a simultaneous drift condition. This is shown in the following lemma.
. . , d, and where Λ is a set of excess switching rates γ : E → R d + such that for all γ ∈ Λ it holds that Let δ > 0 and α > 0 be such that 0 < (δγ max )/V min < α < 1 and define φ(s) = 1 2 sign(s) ln (1 + δ|s|). Then the function is a simultaneous Lyapunov function for the family of ZZ processes, that is there exist Proof. We follow the proof of Lemma 11 in [12]. Let M ∈ M and γ ∈ Λ. Applying the generator to V , which was defined in Equation (19), we obtain For the remaining term the best we can do is to derive the following bound where we have used that 0 ≤ φ ′ (s) ≤ δ/2. Finally we obtain for any M ∈ M which is independent of the specific M and γ and can be made arbitrarily small outside of a sufficiently large compact set C by our assumptions on U .
If we wish to consider a more general class of positive-definite matrices, then we have to settle for the following, weaker result. matrices and Λ is a space of excess switching rates γ : E → R + such that (18) is satisfied for some γ max . Let δ > 0 and α > 0 be such that 0 < δγ max < α < 1. Define for each M ∈ M the function where M i denotes the i-th column of M and φ : R → R was defined in Lemma 4.2. Then there are A 1 > 0, A 2 > 0, and a compact set C ⊂ E such that for all M ∈ M the following simultaneous drift condition holds: for all (x, θ) ∈ E.
In particular A 1 , A 2 , C do not depend neither on M nor on γ (but depend on M and Λ).
Then transforming the function in (21) we obtain a Lyapunov function for the standard ZZ process with transformed target:Ṽ  θ) . In particular the maximum over γ ∈ Λ does not cause problem as the γ's are uniformly bounded. Hence A 2 is independent of the specific M and γ.

Finalising the proof of Theorem 3.8
Let us first consider the case of diagonal preconditioners. Let ∆t > 0 be the discretisation step. Then by Lemma 4.1 for any set of the form C = D × V , with D compact and V ⊆ Θ, there exist δ > 0, n 0 := inf{n ∈ N : n ≥ t 0 ∆t }, ν M (·) such that C is a uniform (ν M , n 0 , δ)-small set for the family of discretised processes. Observe that no conditions on ∆t are required. Moreover, a simultaneous drift condition holds by Lemma 4.2 combined with Lemma B.3 for any ∆t. The condition µ(V ) < ∞ is satisfied by definition of the Lyapunov function V , and in fact it also holds that sup (x,θ)∈C V (x, θ) < ∞ by continuity of V in x. Therefore all the conditions of Assumption 3.2 are verified, which means the family is simultaneously geometrically ergodic and by Theorem 3 in [5] the containment condition is satisfied.
In the case of a non-diagonal transformation matrix, parts (a)-(c) of condition (1)  Proposition 4.6 implies that, under the assumption that p n → 0 as n → ∞, the adaptive strategy satisfies diminishing adaptation. Therefore ergodicity follows from diminishing adaptation and containment.

Proof of Theorem 3.10
In the next two sections we show respectively that condition (2) in Theorem 3.6 is verified for the family of preconditioned BPS and/or of BPS with refreshment rates in a compact set. In Section 4.2.3 we use these auxiliary results to show the theorem.

Simultaneous coupling inequality
The next lemma states that a simultaneous coupling inequality is satisfied for the BPS with adaptive preconditioner and/or adaptive refreshment rate. The proof is based on the proof of Lemma 12 in [18], which shows a coupling inequality result for the standard BPS.
In particular α is independent of M and λ r .
Our strategy is thus to apply Lemma 12 in [18] to each semigroupP t M,λr and then take advantage of Equation (22). Observe that it is possible to apply the lemma because part (a) of Assumption 3.9 implies that for all M ∈ M Hence the integrability condition holds for eachP t M,λr with respect to the corresponding targetπ M . Applying the lemma and Equation (22) it follows that for any compact set Observe that there is a dependence on M due to the fact that the set of initial conditions K is trans- (22). This translates to a dependence on M in the coefficient α, as a consequence of the dependence of α on set K M . This inconvenience can be avoided if we choose α such that the coupling inequality is satisfied for all ini- The last thing to do is showing that there exists a constant α * > 0 independent of M and λ r such that α(λ r , M, K) ≥ α * for all M ∈ M and λ r ∈ Λ r . This is indeed the case for the following reasons: • The dependence on M appears then in the following term, which in the statement of Lemma 12 in [18] appears as a factor in α(λ r , M, K): HereÑ = N + (1 + E 1 /λ r )R M for some N > 0, and E 1 , E 3 are independent exponential random variables with parameter 1. Because g(r) is a factor in α(λ r , M, K), we wish to bound it from below, and hence we should bound the supremum of ∇Ũ M (ξ) from below. Denoting ζ(λ r ) = (1 + E 1 /λ r )R M + (r/λ r )Ñ , this can be done as follows: where we used thatŨ M (ξ) = M T ∇U (x), that M −1 x ≥ x / M , and the compactness of M. This is sufficient to eliminate any dependence on M in α. • Depending on the specific factors in the statement of the lemma, the switching rate λ r can be conveniently bounded either above by λ max r or below by λ min r in order to bound α(λ r , M, K) from below. This can be done in every term in which λ r appears and it follows that the dependence on it can be easily eliminated.

Drift condition for the BPS
The second condition we need is uniformity of the constants in the drift condition for the family of preconditioned BPS and/or for the family of BPS with different refreshment rate. To this end, we go through the proof of Lemma 7 from [18] to show that this is indeed the case. Then there are A 1 , A 2 > 0 and a class of functions {V M,λr : M ∈ M, λ r ∈ Λ r } such that for each M ∈ M and λ r ∈ Λ r it holds that Proof. We follow the same underlying idea that was used in the proof of Lemma 4.3. It is in fact sufficient that there exist A 1 , A 2 > 0, both independent of M and λ r , such that for all M and λ r there exists a functionṼ M,λr (ξ, θ) such that for all (x, θ) ∈ R d × R d . In order to show the inequality in (24) we rely on [18,Lemma 7]. In particular we begin by showing that the constants c 1 , c 2 , c 3 , c 4 > 0, R > 0 defined in [18, Assumption A8] can be chosen independently of the specific M for the class of standard BP samplers with targetsπ M . This implies that A 1 , A 2 can then be chosen to be the same for every M ∈ M. The dependence on λ r is dealt with as a second step. For ease of notation from now on we denote the corresponding energy function as U M (ξ). Let us restrict our attention to the case ψ = N (0, ½ d ) and thus, using the notation in [18], we take ℓ(ξ) = 1, H( θ ) = η θ 2 for some η ∈ (0, 1) small enough such that R d exp (η θ 2 )ψ(dθ) < ∞, andŪ M (ξ) = U ζ M (ξ) for ζ ∈ (0, 1) chosen as in part (c) of Assumption 3.9. Observe that η, ζ are uniform in M .
We start by considering c 1 , which must be such that ∇ ξŪM (ξ) ≥ c 1 for all points outside of a ball of radius R 1 . In particular we wish to show that there exist c 1 , R 1 that satisfy the property above for all M ∈ M. Observing that for x = M ξ it holds that ∇ ξŪM (ξ) = ζU where C = ζ min M ∈M (1/ M −T ) > 0. Assumption 3.9(c) implies that there existc 1 ,R 1 > 0 such that ∇ x U (x) /U 1−ζ (x) ≥c 1 for any x such that x ≥R 1 , and because M is a compact space, we can choose c 1 = Cc 1 and R 1 =R 1 , which are then independent of M . The constant c 2 must be such that ℓ(ξ) ≤ c 2 , and because for a Gaussian ψ, we may take as written above ℓ(ξ) = 1. Therefore c 2 = 1 for each M ∈ M.
Then, c 3 must be such that ( ∇ ξ U M (ξ) / ∇ ξŪM (ξ) ) ≥ c 3 for all ξ such that ξ > R 3 for some By part (b) of Assumption 3.9 we have that lim x →∞ U (x) = +∞ and since ζ ∈ (0, 1) there exist c 3 and R 3 large enough such that outside of a ball of radius It is left to show that c 4 can be chosen uniform.
for any ξ such that ξ > R 4 (27) for some R 4 > 0. Observing that for θ ∈ A x,M it holds that sup θ∈A In the second to last inequality we used that ≤c 4 for all x such that x ≥R 4 .
As a consequence, because M is a compact space and by (28), there exist c 4 , R 4 independent of M such that (27) holds for all M ∈ M. It is now sufficient to take R = max{R 1 , R 3 , R 4 }, and we have shown that c 1 , c 2 , c 3 , c 4 , R can be picked uniformly for all standard BP samplers with energy function U M . In particular there is no dependence on λ r in all these constants.
It is important to observe that part (c) of Assumption 3.9 is satisfied by all BP samplers with targets in {π M } M ∈M because we have bounds on M and M −1 . Moreover, parts (a) and (b) of Assumption 3.9 are trivially verified because of the transformation scheme in Figure 2 that definesŨ M . It is then possible to apply [18,Lemma 7] to each process to obtain the drift condition (24). Therefore for each M ∈ M and λ r ∈ Λ r we have the Lyapunov functioñ where η ∈ (0, 1) and r depend only the distribution at refreshments ψ, κ λr ∈ (0, 1) depends on the c i 's, and ϕ λr is a positive, non decreasing, continuously differentiable function as defined in [18]. Notice that, as a consequence of the calculations above, κ λr and ϕ λr do not depend on the preconditioner. Applying our usual transformation scheme we obtain the functions This means that for eachṼ M as defined in (29), the same proof of [18,Lemma 7] can be followed and thusṼ M satisfies In particular, A 1 , A 2 both depend on λ r , on c 1 , c 2 , c 3 , c 4 , R, and on other constants that depend only on ψ and thus are independent of M . Therefore A 1 , A 2 can be chosen uniformly in M ∈ M. The last step is now to eliminate the dependence on λ r in A 1 , A 2 . Observe that A 1 (λ r ) is defined in the proof of [18, Lemma 7] as a minimum of several constants. It is enough for our needs to notice that the constants have a continuous dependence in λ r and that we are considering λ r ∈ [λ min , λ max ]. On the other hand, A 2 (λ r ) can be chosen independently of λ r by continuity of V M,λr in λ r . We can therefore choose A 1 , A 2 independently both of M and λ r , thus the wanted simultaneous drift condition follows by (25).

Finalising the proof of Theorem 3.10
Let ∆t > 0 be a discretisation step. Lemma 4.5 gives the drift condition, that is condition (b) in Theorem 3.6. Then Lemma 4.4 implies that a coupling inequality holds for any compact set. Sets of the form V M,λr (x, θ) + V M,λr (x, θ) ≤ C 2 are compact by definition of the class of Lyapunov functions {V M,λr : M ∈ M, λ r ∈ Λ r }. We are in particular free to choose the constant C 2 as large as we wish. Note also that the coupling inequality holds for all t > 0, hence there are no constraints on the choice of ∆t. Moreover µ(V M,λr ) =μ M (Ṽ M,λ R ) < ∞ for all M ∈ M and λ r ∈ Λ r . HereṼ M,λ R is the Lyapunov function of a standard BPS with refreshment rate λ r and targetμ M =π M × Ψ. The containment condition is thus verified as all conditions in part (2) of Theorem 3.6 hold. Proposition 4.6 implies the diminishing adaptation condition, and thus ergodicity follows.

Proving the diminishing adaptation condition
A key part of Theorem 3.1 is condition (b), i.e. the diminishing adaptation condition. For the adaptive scheme described in Section 2.4 the condition can be easily shown to be true as the adaptation happens with diminishing probability. Proposition 4.6. Consider the adaptive schemes in Section 2. In particular, assume that {p n } n≥0 , i.e. the sequence of probabilities of updating the adaptation parameters, is such that p n → 0 as n → ∞. Then the diminishing adaptation holds for any t ≥ 0.
Proof. Consider for example the adaptive BPS. Observe that M n+1 = M n and λ n+1 r = λ n r with probability 1 − p n+1 and thus Thus the diminishing adaption holds. The same reasoning works for the adaptive ZZS.

Numerical experiments
In this section we test the empirical performance of the adaptive schemes we defined in Section 2.4. All experiments are implemented in Julia and the corresponding codes can be found at https://github.com/andreabertazzi/Adaptive_PDMC_samplers. Let us state some settings that hold for all experiments below. The time horizon is set to T = 10 5 for all processes. This is large enough for the adaptive PDMC samplers to learn and take advantage of the covariance structure. When considered fixed, the refreshment rate of the BPS is taken to be λ r = 1. The excess switching rate for ZZS is set to 0 in all experiments. The discretisation step is chosen to be ∆t = 0.5, which in our experiments turns out to be a good choice for a wide range of targets. Moreover, we set adaptation times to be every t adap = 2000 continuous time units. The probability of adapting decays as O(log log n). Finally, no normalisation in the sense discussed at the end of Section 2.3 is employed. The performance measures we consider are the effective sample size per second (ESS/sec) for the mean and for the radius statistic t(x) = d i=1 x 2 i . In Sections 5.1 and 5.2 these are computed in continuous time as discussed in [8] by estimating the asymptotic variance with the batch means method, and the variance of the Monte Carlo estimate on the continuous time trajectories of the processes. On the other hand, in Section 5.3 we compute the mean squared error (MSE) in discrete time for the sample mean and radius statistic and take advantage of the fact that in the large time horizon regime the MSE is approximately given by the asymptotic variance of the observable divided by the number of generated samples. This alternative way to compute the ESS avoids poor convergence of the batch means method in the multimodal case. Finally, in all settings we repeat the same task 20 times and report all the results in box-plots.

Multidimensional Gaussian target
In this section we focus on two different kinds of multivariate Gaussian target distributions. The first one, denoted by MG1, has unitary variances and correlation ρ between all components. Denoting the covariance matrix by Σ, this means that Σ ii = 1 for each i = 1, . . . , d and Σ ij = ρ for all i = j. We study how the adaptive PDMC algorithms compare to their non-adaptive counterparts for different values of ρ and different dimensionalities. In this setting we focus on adaptive algorithms that estimate the full covariance matrix. The second Gaussian target we consider has variances 0.5, 1, 5, 10, 15 repeated depending on the dimension, together with a milder correlation between components. This setting is denoted as MG2 and is useful to compare all kinds of adaptive algorithms we introduced.

MG1 target
In the first experiment we consider a 50-dimensional MG1 target for different values of ρ. In Figure 3 the average ESS/sec and the ESS/sec for the radius statistic are shown. As expected, the performance of the ZZS is degrading as the correlation increases. This behaviour is for the most part caused by the very large number of events that have to be simulated for very narrow targets. The adaptive ZZS successfully improves over this inconvenience and is stable with respect to the increasing correlation. The standard BPS with fixed refreshment rate shows a decaying average ESS/sec, while the ESS/sec for the radius statistic appears to increase as ρ grows up until ρ = 0.4 and then becomes smaller. This behaviour is likely due to the fact that the choice λ r = 1 is more suited for the estimation of the radius in case of a more concentrated target rather than for a standard Gaussian. A similar behaviour is shown by the adaptive BPS's. Overall we notice a marked improvement for the BPS with adaptive preconditioner and fixed λ r . Choosing to adapt only the refreshment turns out to be a detrimental decision when the target is correlated. Indeed the optimality criterion derived [10] assumes a standard Gaussian target. In Figure 4 we study how the adaptive schemes compare to the standard ones for an MG1 target with correlation ρ = 0.8 and increasing dimensionalities of target. The plots show that when the target is strongly correlated the effect of the adaptation shows no sign of diminishing. It also seems clear that the sampler of choice in this case should be the adaptive BPS with fixed, rather than adaptive, refreshment rate. This could be due to the fact that, when the refreshment rate is updated adaptively and the target is anisotropic, a too large λ r is chosen at first due to the high number of reflections, thus slowing down the estimation of the covariance matrix. As suggested by a reviewer, one could avoid this issue by keeping the refreshment rate fixed until the estimate of the covariance stabilises, and only then starting to learn the optimal λ r . It is worth pointing out that the performance of the BPS with adaptive preconditioner and refreshment improves compared to the BPS as the dimension increases. This is according to the theoretical results in [10], which are indeed obtained in the high dimensional   limit. Therefore we expect that for a large d it is reasonable to apply both the transformation scheme and the tuning of λ r .

MG2 target
Let us now consider a 50-dimensional MG2 target with a mild correlation set to ρ = 0.3. Figure 5 shows the results for several adaptive PDMC samplers. The adaptive algorithms that learn the entire covariance matrix show the largest gain in terms of ESS/sec. For the BPS the choice of learning only the variance of each component of the target seems interesting, also in view of larger dimensions. As in the previous section, we observe that the adaptation of the refreshment rate seems to have a bad effect for anisotropic targets.

Logistic regression with correlated data
The next numerical experiment we consider is a Bayesian logistic regression task. In this setting for j = 1, . . . , n obs a binary output value y j ∈ {0, 1} has distribution where {x j } n obs j=1 are known covariates, and β ∈ R d is an unknown parameter. We take a flat prior and thus obtain the posterior We force correlation between some components of the parameter by taking, for j = 1, . . . , n obs and i = 1, . . . , d, (x j ) i = 1 + εN ji , where N ji ∼ N (0, 1) and ε = 0.1. The results of the experiment are reported in Figure 6, in which the samplers are tested with targets as above with d = 2, 4, 8, 16 and n obs = 1000. The adaptation of the refreshment rate follows the alternative scheme discussed in Appendix A.3. This scheme seems more stable as the refreshment rate is update gradually and cannot jump immediately to very large or small values. Although the dimensionality is small and the correlation is limited to a subset of the coordinates, we observe that the adaptive schemes outperform their standard counterparts.

Mixture of Gaussian distributions
Consider a mixture of two 30-dimensional Gaussian distributions N (0 d , Σ) and N (µ, Σ), both with weight 1 2 . Here we take Σ with Σ ii = 1 and Σ ij = 0.25 for j = i. Moreover we take µ = a × (1, . . . , 1), where a is a parameter that determines the distance between the two means. We investigate the performance of the adaptive schemes compared to standard ones as a function of the parameter a. The details on the implementation of this experiment can be found in Appendix A.4. In this experiment we test the robustness of the algorithm in a case in which the target distribution is cigar shaped, but multimodal. The time horizon for ZZS is 10 5 , while for BPS it is 4 × 10 5 . Figure 7 shows the results for the adaptive samplers in this setting. We observe that as the distance between the means increases the performance of the adaptive samplers improves over the standard ones. We remark that for small values of the distance a one cannot expect improvements of the adaptive algorithms as the correlation is small and thus the target is not very anisotropic.

Discussion
In this paper we proposed adaptive schemes to overcome two of the current issues with the BPS and ZZS. We have shown that the refreshment rate and the excess switching rate can be tuned on the fly as long as the updates or the probabilities of updating decrease to 0. With this approach the user does not have to worry about tuning the refreshment rate. A current limitation is that more theory or experiments are needed to determine a criterion that works well with anisotropic targets. In addition to this, we have proposed a way to make the PDMC samplers learn and take advantage of the covariance structure of the target. The theoretical results stated in Theorems 3.8 and 3.10 ensure that the adaptive samplers are ergodic and thus converge to the correct measure π. It is challenging to prove theoretical statements regarding performance improvements of the adaptive schemes over the standard PDMC samplers. However, the numerical experiments we conducted suggest that our adaptive algorithms can lead to a significant performance improvement when there are strong anisotropies. An alternative approach could be to use an optimisation algorithm to obtain an estimate of the covariance matrix by computing the Hessian of the target at its point of maximum. However, this optimisation step can be expensive and moreover for realistic problems it is not a given that this estimator is a good approximation of the posterior covariance. In particular in our theory we do not assume log convexity   of the target, and also we do not assume that we are in the large sample regime where a Bernsteinvon Mises theorem holds. In addition, the adaptive schemes discussed in this paper can be applied to the Boomerang sampler [9]. This would result in elliptical dynamics that are adapted to resemble the (unimodal) target at hand. Naturally, the adaptive algorithms should be run with a time horizon that is large enough to benefit from the adaptation. Two other important settings are the discretisation step and the time between two adaptations. Based on our experience with the experiments, we suggest ∆t = O(10 −1 ) and t adaps = O(10 3 ). For concentrated targets, as for instance posteriors when there is a very large number of data points, it is suggested to choose both values small. We remark that in very high dimensional settings it may be unfeasible to let the samplers learn the full covariance matrix, as the computation of M entails calculating the square root of the empirical covariance matrix. In such cases we suggest either learning only the diagonal elements or blocks of the covariance. We remark that the adaptive PDMC algorithms with subsampling are applicable in the setting of tall data, that is when data-set is made of a large number of observations, but with a moderate dimensionality. In such settings subsampling can be shown empirically to result in an improved efficiency (in terms of ESS per second); a result which is backed by a heuristic argument, based on posterior contraction, i.e., the Bernstein-Von Mises theorem; see [8] for details. More research is necessary to understand in which situations subsampling can lead to improved efficiency, and in particular if improved efficiency is possible in cases for which the Bernstein-von Mises theorem does not apply; see also [6,25] for a critical discussion of subsampling methods. A question that one could naturally ask is how applicable this transformation scheme is in case of a multimodal target. The answer depends on the specific target at hand, but one can design a target as a mixture of Gaussian distributions for which applying the transformation scheme would not speed up the convergence of the sampler. However, when the target is multimodal it is possible for instance to use the adaptive PDMC samplers together with the framework proposed in [32]. In this framework, the adaptive PDMC samplers would be beneficial since the regions around each mode would be explored more efficiently by taking advantage of the covariance structure of the specific mode.
Finally, we remark that the idea of learning the covariance structure of the target on the fly could be applied to obtain adaptive versions of the Hamiltonian Monte Carlo (HMC) algorithm [30] and of the Metropolis Adjusted Langevin Algorithm (MALA) [37]. In particular both the HMC algorithm and the MALA are sensitive to correlation in the target and can thus benefit from a suitable preconditioner, as argued respectively in Section 4.1 of [30] and in [36]. Moreover, the preconditioner could be chosen to take advantage of the geometry of the target, as proposed in [23] for HMC and MALA. The preconditioner could be estimated adaptively with an appropriate adaptation strategy, together with similar ideas presented in this manuscript.

Appendix A: Implementation of adaptive PDMC algorithms
The main issue in the exact simulation of PDMPs lies in the simulation of the switching times. It is generally proposed to use Poisson thinning to overcome this difficulty. The idea is to find upper bounds for the switching rates that are more tractable, then simulate a Poisson process with said rate and finally adjust with an acceptance-rejection step. Applying a preconditioning matrix to PDMPs results in a modification of the switching rates and thus the bound proposed in [8,13] do not directly apply to our proposed algorithms. In the following two sections we give two important examples to motivate that it is possible to find bounds with similar ideas to those used for standard PDMC algorithms. In Appendix A.3 we discuss a different adaptation strategy for the refreshment rate.

A.1. Dominated Hessian of the negative log-likelihood
Let us consider the case in which the Hessian of the negative log-likelihood is dominated by a positive definite matrix. Denoting the Hessian by H U (x) = (∂ i ∂ j U (x)) d i,j=1 , this means that for any x ∈ R d it holds that H U (x)u, u ≤ Qu, u for all u ∈ R d . Then we say that −Q H U (x) Q. Assuming Q is symmetric, we have for all u, v ∈ R d that u, H U (x)v ≤ u 2 Qv 2 . Let us consider the switching rates of a ZZ process with preconditioner M , i.e. λ M i (x, θ) = (θ i M i , ∇U (x) ) + , where M i is the i-th column of M . Remember that the deterministic trajectory of this process is x(t) = x + M θt, where (x, θ) is the initial condition. We have where e i , is the i-th vector of the canonical basis, i.e. with zeros in all components except for a 1 in the i-th component. Taking It is possible to have a tighter bound taking b i = QM i 2 M θ 2 , but this entails having to update b i after every switching time, which would make the overall simulation more expensive. Note that in particular for a diagonal M we have that M θ 2 = ( d i=1 M 2 ii ) 1/2 and we can take It is worth observing that Q 2 can be computed once before running the algorithm and then stored. Terms M 2 and M i 2 can in principle be estimated by taking a maximum over the class of matrices M. However, having loose computational bounds leads to a lower percentage of accepted proposed events and thus to an increased computational burden.
For BPS with preconditioner M we showed that λ M (x, θ) = ( M θ, ∇U (x) ) + . With computations similar to the ones above we obtain

A.2. Subsampling techniques
PDMC algorithms are particularly interesting because they allow for exact subsampling. This means that it is possible to modify PDMC samplers such that the correct invariant measure is maintained but without going through the entire data-set at each iteration. This was illustrated in [8,13] for ZZS and BPS. Subsampling techniques are particularly helpful when the number of data points n is much larger than the dimensionality d of the posterior density function. This is an interesting scenario for the adaptive algorithms that are here introduced, as for a small d the additional computations necessary to learn and take the square root of (part of) the covariance matrix are not overpowering. In this section we explain how adaptive PDMC with subsampling can be implemented.
Suppose the partial derivatives of the posterior distribution can be written in the form where E j i : R d → R d are continuous mappings. This holds for example when the target satisfies This is the case for instance when the target is a posterior density of a model with iid observations. In such cases one can choose U j (x) = − log(π 0 (x)) − n log (l(y j |x)), where π 0 is a prior of the parameters and l is the likelihood associated to observation y j . Then the idea is to define a collection of switching rates along a trajectory.
for all j, then the subsampling procedure is successful. One can simulate IPPs with rates H i (t) and then accept the proposed switch of component i 0 with probability m J i 0 (t)/H i 0 (t) where J ∼ Unif{1, . . . , n}. The next proposition shows that stationarity of π is preserved for the ZZS with subsampling also introducing a preconditioner.
Proposition A.1. The Zig-Zag sampler with subsampling and preconditioner M has invariant distribution µ = π × Unif(Θ). Moreover it coincides with a Zig-Zag process with switching rates Proof. The statements can be derived following the proof of Theorem 4.1 in [8] and is thus omitted.
An interesting situation is that of a negative log-likelihood with Lipschitz partial derivatives, i.e.
In this case we can take a reference point x * and take This allows us to obtain the following uniform bound on the switching rates Therefore we can simulate d IPPs with rates The same estimator for the gradient of the negative log-density can be used for a preconditioned BPS. With computations analogous to above we find Although establishing ergodicity of these adaptive algorithms would be an interesting research question, we leave it for future research due to the lack of theoretical results for their standard counterparts.

A.3. A different adaptation strategy for the refreshment rate
The adaptive scheme for the refreshment rate described in Section 2.4 is a natural implementation of the results in [10]. However, it can be unstable when combined with an adaptive preconditioner. Here we define an alternative adaptive strategy and we prove that the diminishing adaptation condition from Theorem 3.5 is satisfied.
Assume that the refreshment rate of the BPS is constant in the position and velocity spaces. We update it iteratively as follows: in which n refresh (n) and n events (n) are respectively the number of refreshments and the number of events up to time n, q n is a positive, decreasing sequence such that q n → 0, and λ * = 0.7812. For the strategy in Section 2.4 we used that at time n there is a probability of adapting equal to p n , such that p n → 0. In the remainder of this section we show that for the adaptive scheme described in (33) the diminishing adaptation condition holds even if the refreshment rate is updated with probability 1.
Lemma A.2. Let π be a probability distribution. Denote as P t γ the semigroup of a ZZ process with excess switching rate γ : E → R d + . Let γ 1 , γ 2 be two bounded excess switching rates. Then for any Similarly, denote now by P t λ the semigroup of a BPS with refreshment rate λ : E → R + . Let λ 1 , λ 2 be two bounded refreshment rates. Then for any t ≥ 0 Proof. Consider first the case of the BPS. Consider the generators of two BPS's with refreshment rates λ 1 and λ 2 and denote them respectively by L λ 1 and L λ 2 . Now observe that L λ 2 = L λ 1 + B, where B is a bounded operator such that In other words, L λ 2 is a perturbed version of L λ 1 , with bounded perturbation B. Since L λ 1 and L λ 2 generate strongly continuous semigroups (P t λ 1 ) t≥0 , (P t λ 2 ) t≥0 , we can apply the reasoning in the proof of Corollary 1.11 in Chapter 3 of [19] to obtain that, for any t > 0, P t λ 1 − P t λ 2 → 0 as B → 0, which is equivalent to where C 0 (E) is the space of continuous functions on E that vanish at infinity. Inverting the order of the two suprema and applying the Riesz-Markov representation theorem (see e.g. [4]) we obtain the result in (35). Now consider the case of two ZZS with two different bounded excess switching rates γ 1 : The perturbation is now given bỹ Then the result (34) follows by the same arguments as above. Proof. By the triangle inequality for the total variation distance we obtain The inequality above allows us to deal with the two adaptations separately. The second term in the right hand side of (36) is handled by Proposition 4.6. Now focus on the first term in the right hand side of (36), for which we want to apply Lemma A.2. It is sufficient to observe that supremum norm of the perturbation goes to zero in probability when n → ∞. Indeed the sequence of refreshment rates is Cauchy in probability, as for all n we have λ n+1 r − λ n r ∞ ≤ δ n+1 . Convergence in probability follows from the fact that (δ n ) n≥0 is a convergent sequence and therefore the diminishing adaptation condition holds.

A.4. Details on the implementation of a Gaussian mixture target
Consider the target of Section 5.3, that is Here we wish to find upper and lower bounds to the Hessian of the negative log-density, which can be used to simulate the adaptive ZZS and BPS following the approach of Appendix A.1. Start by writing w = 1 2 (µ 2 + µ 1 ) and v = 1 2 (µ 2 − µ 1 ). Then we have and and It remains to find the minimum of m → −4αβ (α exp(−m) + β exp(m)) 2 .
Proof of Proposition 2.2. Following the approach in the proof of Theorem 2.2 in [8], we want to check that for any M ∈ M, and any h ∈ D(L M ), it holds that E L M h(x, v)dµ = 0. By a change of variable and using the same reasoning as in the proof of Proposition 2.1, for any the last equality was obtained by the invariance ofμ M for the standard Zig-Zag process.
Proof of Proposition 2.3. Following the same approach as in the proof of Proposition 2.1, for anỹ f ∈ D(L M ) we define a function h(x, θ) =f (ξ, θ) with x = M ξ. The generator of the transformed process then satisfies Condition (38) provided that In this caseL whereμ M =π M ⊗ ψ and the last equality follows by the invariance ofμ M for the BPS as shown in [13]. Let γ ∈ Y and let C be the uniform (ν γ , δ, n 0 )-small set as defined in Assumption 3.3. As described in [33,38] we define the following coupling between two chains with kernel P γ . The two chains evolve independently each with kernel P γ when they are outside of C. When both chains are inside C, then with probability δ we let the two chains couple after n 0 steps by drawing X n+n 0 = Y n+n 0 from ν γ , or with probability (1 − δ) we independently draw X n+n 0 ∼ (P n 0 γ (X n , ·) − δν γ (·))/(1 − δ) and Y n+n 0 ∼ (P n 0 γ (Y n , ·) − δν γ (·))/(1 − δ). Then by taking h γ (x, y) = (V γ (x) + V γ (y))/2 as suggested in the proof of Theorem 3 in [5] we have for all (x, y) / ∈ C × C.
Define now A(γ) := sup (x,y)∈C×C E γ h γ (X n 0 , Y n 0 ) X 0 = x, Y 0 = y . Then by [38,Theorem 5] we have for each γ ∈ Y the following bound on the total variation distance The dependence on γ in A(γ) can be eliminated by noting that which is independent of γ. In particular by our assumptions we have that B < ∞ and π(V γ ) < ∞ for each γ.

B.2.3. Boundedness in probability of the sequence of Lyapunov functions
In both cases the containment condition (assumption (a) in Theorem 3.1) is thus satisfied if the process {V Γn (Z n )} n≥0 is bounded in probability. This is indeed implied by our assumption that Γ n is updated only if Z n ∈ B, where B is a compact set, as suggested in [15,14,32]. For the sake of completeness we report the main steps. By [34,Lemma 3] it is enough to show that sup n∈N E(V Γn (Z n )) < ∞. By our assumptions and letting D : Taking expectations on both sides one obtains This shows that the sequence is contracting since λ ∈ (0, 1) and this is enough by [34, Lemma 2] to show our boundedness in probability of the sequence.

B.3. Proof of Theorem 3.6
To prove Theorem 3.6 it is sufficient to verify that the discretised process with time step ∆t satisfies the assumptions of Theorem 3.5. Let us first consider the assumptions in alternative (1) of Theorem 3.6. The uniform small set condition of the discretised process (Assumption 3.3(a)) follows by choosing n 0 = t 0 ∆t . The geometric drift condition (Assumption 3.3(b)) is a consequence of Lemma B.3. Now consider the second set of assumptions, that is alternative (2) in the theorem. In this case we wish to show that Assumption 3.4 holds for the discretised process. Part (a) of Assumption 3.4 is immediately verified by the fact that ∆t = t 0 . For part (b), first observe that from Lemma B.3 by our assumptions we have In the notation of Assumption 3.4(b) we have λ = e −A 1 ∆t and C 1 = A 2 A 1 . In particular we have 2C 1 < C 2 since we assumed that The thesis now follows in both cases by Theorem 3.5.

B.4.1. Small set condition for a family of one-dimensional ZZ processes
In this section we show a simultaneous minorisation condition for the ZZ process in the onedimensional case. This result is essential to generalise the result to the multidimensional case (see Lemma 4.1 and its proof).
Assumption B.1 (Assumption 3 in [11]). Let U ∈ C 2 (R). Define the switching rates of a 1-d ZZ process as λ(x, θ) = (θU ′ (x)) + , where x ∈ R and θ ∈ {−1, +1}. There exists x 0 > 0 such that Remark. Assumption B.1 requires that there exists a point x 0 ∈ R such that the process has a strictly positive probability of changing direction if it is outside of and moving outwards of the set [−x 0 , +x 0 ]. The lemma below states that a small set condition follows from this assumption.
Assume either of the two conditions: • γ(x) = 0 for all x ∈ R, but Assumption B.1 holds; • there exists γ min > 0 such that γ(x) ≥ γ min for any x ∈ R, and γ(·) is bounded on compact sets.
Then for any set of the form C = D × V , where D ⊂ R is a compact set and V ⊆ {−1, +1}, there exists t 0 > 0 such that for any t ≥ t 0 there are δ > 0, and a probability measure ν on E such that P t m ((x, θ), ·) ≥ δν(·) for all (x, θ) ∈ C and all m ∈ M.
In particular t 0 depends on V min , and δ depends on V min , V max , γ min , γ max , but neither depends on m or γ. For m ∈ [V min , V max ], the process L m moves with velocity mθ, so either +m > 0 or −m < 0. Consider first the case in which Assumption B.1 is satisfied for some x 0 > 0. Note that as a consequence it is satisfied for any L m . We consider the case x 0 > R because it is the most general. Indeed one is always free to takex 0 > R > x 0 and apply the reasoning below. Alternatively the case x 0 ≤ R follows by the same method of proof.
) be a probability measure on E. Define T = T (ε) := (4x 0 + 2R)/V min + ε where ε > 0. This choice of the time horizon allows the process to reach the region with strictly positive probability of a velocity flip for any initial position (x, θ) with x ∈ [−R, R] and θ ∈ {−1, +1}. The addition of ε > 0 is necessary to introduce a margin where the flip can take place. See Figure 8 for a visual aid.
We show that C is a uniform (ν, δ, T )-small set. Since the proof applies for all ε > 0, although resulting indifferent δ's, one is free to choose t 0 = T (ε) freely. Let us consider the two cases below. First case: let the initial condition be (x, +1) with x ∈ [−R, R] and we want to be at time T in B = A × {−1} with A ∈ B([−R, +R]) any Borel set. We can then use the following inequality P (x,+1) ((X(T ), Θ(T )) ∈ B) ≥ P (x,+1) (X(T ) ∈ A, E 1 ), where E 1 is the event that exactly one velocity switch takes place. We are then in the case of Figure 8(a). Observe that by the choice of T the process has enough time to travel the longest path, i.e. from (−R, +1) to (−R, −1) with the smallest allowed velocity V min . In order to compute the r.h.s. in (41) one can compute P (x,+1) (X(T ) ≤ y, E 1 ) and then differentiate with respect to y. To do this we assume only one velocity switch and impose that X(T ) ≤ y, resulting in the condition where t is the time at which the velocity switch takes place. By rearranging we obtain the condition t ≤ y−x 2m + T 2 =:t(y). Observe that for any m ∈ [V min , V max ] the process has enough time to reach the region where it is assumed there is a strictly positive probability of a velocity flip. Indeed using that y ∈ [−R, R] we find whereC T := [−R − V max T, R + V max T ] is the set of points that can be reached in time T by a process starting in x ∈ [−R, R] with velocity V max (and thus for any value m ∈ [V min , V max ]). Thus λ max is the maximum switching rate achieved within time T .
Second case: again we consider an initial condition with positive velocity, but now also a positive velocity at time T , i.e. we consider a set B = A × {+1}. Taking advantage of the same idea as before we use the bound P (x,+1) ((X(T ), Θ(T )) ∈ B) ≥ P (x,+1) (X(T ) ∈ A, E 2 ), where E 2 is the event that exactly two switches take place. This case corresponds to Figure 8(b). Let t 1 and t 2 be the times of first and second switch respectively. Then X(T ) ≤ y when X(T ) = x + mt 1 − m(t 2 − t 1 ) + m(T − t 2 ) ≤ y, which can be rearranged as t 2 ≥ t 1 + x−y 2m + T 2 =:t 2 (y). We can obtain a bound for t 1 by imposing thatt 2 (y) < T . The resulting condition is t 1 < T /2 − (x − y)/(2m) =:t 1 (y). We observe that by definition of T it follows that for any x, y ∈ [−R, R] which means that there is enough time for the process to reach cyan region in Figure 8 λ(x + ws, 1)λ(x + ws − wt 2 (y), −1) ds dy.
Since the integrand is non-negative we can lower bound this quantity by restricting the domain of integration corresponding to the s variable. A sensible choice is (x 0 − x)/m ≤ s ≤ −(x 0 + x)/m + t 2 (y), as this would imply that both switches take place in the cyan region in Figure 8(b). Indeed for s ∈ [(x 0 − x)/m, −(x 0 + x)/m +t 2 (y)] we have that x + ms ≥ x 0 and x + ms − mt 2 (y) ≤ −x 0 . Observe that where we used that −x 0 ≤ y since y ∈ [−R, R]. Combining this with the fact thatt 1 (y) > x 0 −x m , we have shown that [(x 0 − x)/m, −(x 0 + x)/m +t 2 (y)] ⊂ [0,t 1 (y)]. Therefore on this interval we can bound below the switching rates as follows x 0 +x m +t 2 (y) x 0 −x m λ(x + ms, 1)λ(x + ms − wt 2 (y), −1)dsdy By symmetry the same bounds hold also when the process has initial velocity is −1. Therefore for any ε > 0 we proved that C = [−R, +R] × {−1, +1} is a uniform (ν, δ, T )-small set with Notice that there is an ε dependence also in λ min and λ max , but there is no dependence on the specific value of m. To conclude the argument we observe that it is always possible to incorporate any compact set in a set of the same form as C and therefore all compact sets are uniformly small for the family of Zig-Zag processes. When the switching rate is strictly positive but Assumption B.1 does not hold, the same reasoning can be applied by taking x 0 = 0, λ min = γ min > 0, and T = 2R V min + ε for some ε > 0. The switching rates can be upper-bounded because γ is bounded on compact sets. Thus the only difference is that in this case the process can switch at any time and does not need to escape a compact set to do so. One can again obtain the uniform small set condition with δ defined as in (43), and ν as above.

B.4.2. Simultaneous drift conditions, from discrete time to continuous time
The lemma below is needed to go from simultaneous, continuous time drift conditions to simultaneous, discrete time ones.
for all z ∈ E.
b) the family of processes is such that for K ⊂ E as in the previous point, a constant t > 0, and all z ∈ E, there exists a setK(t) ⊂ E which depends on t such that on the event z / ∈K(t) it holds that Z(t) / ∈ K almost surely for all γ ∈ Y.
Proof. Let t > 0 and γ ∈ Y. As a first step, we apply Dynkin's formula to f (z, t) = exp (c 1 t)V γ (z) as in the proof of Theorem 6.1 in [27] to obtain e c 1 t P t γ V γ (z) = V γ (z) + t 0 P s γ ∂ ∂s + L γ e c 1 s V γ (z) ds.
By the product rule and the drift condition we can write the integrand on the right hand side as In the second inequality we took advantage of condition (b) to conclude that P s γ ½ K (z) = P z (Z(s) ∈ K) ≤ ½K(z). It is then sufficient to take λ = e −c 1 t , κ = c 2 c 1 (1 − e −c 1 t ) to conclude the proof. Observe that λ and κ do not depend on γ (but depend on t).
Remark. Condition (b) in Lemma B.3 is for instance satisfied by a family of preconditioned Zig-Zag processes, if the preconditioner is taken from a compact set of positive definite matrices. Indeed this claim follows from the fact that all processes travel with almost surely bounded velocity. Therefore, it is possible to chooseK(t) by adding a suitable buffer zone around the set K, such that the set K is not reachable in time t starting outsideK(t).
Remark. Instead of drift conditions of the form (44), consider the following case for all z ∈ E.
Without condition (b) in Lemma B.3 we can still conclude by the same line of reasoning that there are λ ∈ (0, 1), κ > 0, both independent of γ, such that for each γ ∈ Y P t γ V γ (z) ≤ λV γ (z) + κ for all z ∈ E.