{'fraction_non_alphanumeric': 0.041047543541085904, 'fraction_numerical': 0.008516406230019693, 'mean_word_length': 4.361579596873715, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Federated Bayesian Computation via Piecewise Deterministic Markov Processes

When performing Bayesian computations in practice, one is often faced with the challenge that the constituent model components and/or the data are only available in a distributed fashion, e.g. due to privacy concerns or sheer volume. While various methods have been proposed for performing posterior inference in such federated settings, these either make very strong assumptions on the data and/or model or otherwise introduce significant bias when the local posteriors are combined to form an approximation of the target posterior. By leveraging recently developed methods for Markov Chain Monte Carlo (MCMC) based on Piecewise Deterministic Markov Processes (PDMPs), we develop a computation -- and communication -- efficient family of posterior inference algorithms (Fed-PDMC) which provides asymptotically exact approximations of the full posterior over a large class of Bayesian models, allowing heterogenous model and data contributions from each client. We show that communication between clients and the server preserves the privacy of the individual data sources by establishing differential privacy guarantees. We quantify the performance of Fed-PDMC over a class of illustrative analytical case-studies and demonstrate its efficacy on a number of synthetic examples along with realistic Bayesian computation benchmarks.


Introduction
In the problem of Federated Bayesian learning we are faced with the unique challenge that, either due to privacy or scalability, the model and its data are distributed across a federation of workers. In this setting, the model and/or data owned by the individual worker must not be disclosed to the other workers, and neither to any coordinator. While this problem has been studied in previous works, much of the proposed methodology has involved sacrificing the asymptotic exactness, which is characteristic of MCMC-based sampling algorithms, to facilitate the distribution and federation of the data, or alternatively being focused on a very narrow class of models.
The problem of federated learning has largely been studied in the optimisation setting, where a global loss function which can be decomposed into local worker contributions must be optimised. Classical strategies including Fed-SGD and Fed-Avg [MMR + 17], FedAc [YM20], and subsequent extensions, build on the general idea that each worker seeks to locally move towards the optimiser of the model based on their own local data contribution, computing a local gradient, and a central server aggregates these local gradient candidates in an appropriate fashion to obtain an estimate of the global minimiser.
The validity of these approaches typically hinge on strong assumptions. Firstly, it is assumed that the data across workers is independent and identically distributed, though a number of more recent federated learning methods seek to weaken these requirements, e.g. through knowledge distillation [ZHZ21] or Bayesian non-parametric modelling, e.g. [YAG + 19]. Secondly, it is assumed that workers use the same local model, though recent work on model personalization has suggested some strategies to address this [MMRS20].
In this work, we address the problem of federated learning in a Bayesian setting, i.e. we seek to generate samples from a global posterior probability distribution obtained as a multiplicative composition of local posteriors distributed across the workers, and without sharing of model and/or data. This is an inherently more challenging problem due to the fact that far more information pertaining to the local posterior distribution must be somehow communicated with other workers while ensuring privacy of data/model. Previous works have sought to lift methodology from the federated learning to the Bayesian setting, employing Stochastic Gradient Langevin Dynamics (SGLD)-based generalisations of Federated Learning counterparts, e.g. [ZLZ + 19, EMMBK21, VPD + 22]. Similarly, the Langevin-type algorithm proposed in [SSR22] combines distributed MCMC with compression techniques to reduce the burden of communicating large gradients.
Related approaches seek to employ approximations of the local posterior contributions which are used to communicate information to the central server. Such approaches include distributed variational inference [ZBKM18, HWL + 17], using Gaussian approximations [ASGXR20] and ensemble approaches [LAD + 21]. In [BGR22], local predictive posterior contributions are distilled and stored into a neural network which is communicated with the central server.
Some works seek to reformulate the federated learning problem through the lens of Bayesian model averaging, where the local model contributions are combined into an accurate global approximation as a model ensemble [CC20,TG20], building on other Bayesian uncertainty quantification methods used in deep learning such as [MIG + 19]. Related to this are approaches which adopt a Bayesian hierarchical modelling view of Federated Learning, introducing hierarchical priors and fixed and random effects to share global information across the different federated workers, [KVMD22].
All of the approaches discussed above either employ local posterior approximations to enable effective communication, and/or are contingent on very strict approximations on the structure of the model. To our knowledge, there is no approach which can perform a full, (asymptotically) exact Bayesian inference in this context for a general Bayesian model. In this paper we provide a federated (or distributed) approach to Markov Chain Monte Carlo with the following properties: (i) the correct posterior distribution is retained; (ii) the observational data may be distributed among workers with no requirement to exchange information other than the algorithmic output; (iii) the observational data amongst different workers does not have to be identically distributed, nor do the local prior distributions have to be the same; (iv) the efficiency of the federated approach compares favourably to the non-federated approach in the sense that the algorithmic slowing down is compensated by the fact that computation is distributed among workers; (v) the amount of information that is communicated between the workers and the server respects the privacy requirements of the worker, which can be quantified from a differential privacy viewpoint.
We will base our approach on the framework of Piecewise Deterministic Monte Carlo [BVD17,BFR19], which we will introduce in Section 2. As discussed in Section 3 this framework can be easily extended to allow for a federated (or distributed) approach while retaining the correct stationary distribution. The computational efficiency of our method is discussed in Section 4. We will also consider our approach from the viewpoint of differential privacy in Section 5. We provide numerical experiments for several examples to establish proof of concept and investigate efficiency properties in Section 6.

Piecewise Deterministic Monte Carlo
Here we concisely describe the framework of Piecewise Deterministic Monte Carlo (PDMC) in some generality. Essentially, a PDMC sampler is based on a piecewise deterministic Markov process. This is a continuous time Markov process which moves along continuous deterministic trajectories, until at random times, a jump within the state space is made. In PDMC the state space consists of a position process X(t) ∈ R d and a velocity process V (t) taking values in a set V ⊂ R d . The jumps (or events) will only affect the velocity. The process will be designed to have a particular stationary probability distribution µ(dx, dv) with marginal position distribution π(dx) = v∈V µ(dx, dv). Here π may be considered to be a Bayesian posterior distribution of interest.

Deterministic dynamics
In the general setting, the deterministic dynamics are described as the solution of an ordinary differential equation where ψ : R d → R is a sufficiently regular function so that solutions to (1) are defined uniquely, e.g., ψ may be assumed to be globally Lipschitz. The deterministic dynamics are assumed to preserve a 'reference' stationary measure µ 0 (dx, dv) = π 0 (dx) ⊗ ν(dv), where π 0 (dx) = exp(−U 0 (x)) dx for a suitable function U 0 : R d → R and ν is a probability measure on V. This means that for a solution φ(t; An interesting special case is when π 0 (dx) is chosen to be the prior distribution in a Bayesian inference problem, but this is not necessary.
Example 2.2 (Bouncy Particle Sampler and Boomerang Sampler). Let V = R d equipped with a Gaussian measure ν = N (0, Σ). where Σ is a positive definite matrix. The dynamics (1) preserve µ 0 by taking ψ(x) = −Σ∇U 0 (x). In particular, for the Bouncy Particle Sampler (BPS, [BVD17] we take U 0 (x) = 0 and thus ψ(x) = 0, and usually Σ = I d so that . In contrast to the ZZS, the BPS has a continuous space of possible velocities, but as for ZZS, the velocities do not change under the deterministic dynamics. For the Boomerang Sampler the deterministic dynamics correspond to a (skewed) harmonic oscillator.

Jumps
Next we specify the jumping mechanism which changes the velocity at random times. This is governed by a jump intensity λ(x, v) and a Markov jump kernel Q(x, v, dv ) : More generally we may have multiple types of jumps with multiple types of rates (λ i ) k i=1 and jump distributions (Q i ) k i=1 , competing for which event occurs first. Suppose we start from time 0 at position (x 0 , v 0 ) and recall that we have deterministic dynamics t → φ(t; x 0 , v 0 ). The distribution of the inter-event times τ i are given by The event that actually takes place is specified by setting i 0 = arg min i τ i . At time τ i0 we make a transition according to the selected jump kernel Q i0 , so that the distribution of the velocity after the jump is given by Remark 2.3. It is always possible to write a combination of jump mechanisms (λ i , Q i ) k i=1 as a single jump mechanism (λ, Q) by defining This provides a convenient notational simplification which we will use whenever this does not cause confusion.
Using the notation of the previous remark, Algorithm 1 describes a general PDMC sampler. In practice it may be challenging to simulate τ satisfying (2). We discuss the usual approach of Poisson thinning in the Appendix.

Stationary distribution
It is possible to formulate conditions on λ i and Q i in order to have a prespecified stationary distribution. For a function f : Suppose we wish to have the distribution exp(−U (x))π 0 (dx) as (marginal) stationary distribution. In order to achieve this we impose the following conditions (understood to hold for all bounded measurable f : R d × V → R): (i) Invariance of ν under the jump kernels: for all 2 B(V) denotes the σ-field of Borel subsets of V.

4:
Set Simulate V (T k+1 ) ∼ Q(X k+1 ,Ṽ k+1 , ·) 6: end for (ii) Effective sign reversal under jumps: for all i = 1, . . . , k and and (iii) Event intensity condition: for all ( Under the stated conditions it follows that the process with deterministic dynamics (1), and jumps according to (λ i , Q i ) m i=1 , has stationary distribution µ(dx, dv) ∝ exp(−U (x)) π 0 (dx) ⊗ ν(dv). The proof of this result depends on the notion of the Markov process generator and is beyond the scope of this work; see e.g. [BFR19,BVD17,BGKR20] Example 2.4 (Zig-Zag Sampler). For Zig-Zag, we consider for i = 1, . . . , d, where γ i is a non-negative function (called the excess switching rate or refreshment rate This corresponds to flipping the ith direction of the velocity at a rate which depends on the ith partial derivative of U as indicated. Example 2.5 (Bouncy Particle Sampler and Boomerang Sampler). For BPS and Boomerang, recall that ν(dv) ∝ N (0, Σ) for a positive definite matrix Σ. We take and This corresponds to a reflection of the velocity in the direction of the gradient of U . In addition, we require a refreshment jump at rate λ 0 (x), which independently draws a new velocity from the distribution ν: without this refreshment the process will in general not be ergodic, i.e., it will not explore the full state space.

The output of a PDMC algorithm
In order to determine the full continuous time trajectory, it is sufficient to determine the positions and velocities (X k , V k ) immediately after jumps. These points are called the skeleton points. The continuous time trajectory is obtained by the deterministic dynamics originating from the skeleton points, as Provided the piecewise deterministic process is ergodic (as discussed in e.g. [BVD17, BRZ19] for BPS and Zig-Zag respectively) we have the following approximation for our computation of interest: with probability one, were (X(t), V (t)) is any random realization of the piecewise deterministic process with characteristics (φ, λ, Q) and arbitrary initial condition. Due to the piecewise linear nature of the trajectories of (X(t)) t≥0 it is often very straightforward to evaluate the one-dimensional integrals in this expression. Alternatively, one can obtain a discrete set of samples ( X k ) k∈N by setting ( X k , V k ) = (X(kδ), V (kδ)), for some arbitrary δ > 0. In this case the usual MCMC approximation formula is satisfied with probability one, because ( X k , V k ) can be seen as a discrete time ergodic Markov chain in R d × V with marginal invariant density on R d given by π.

Federated Piecewise Deterministic Monte Carlo
Now consider the setting in which π(dx) admits the factorization We will distribute the simulation of π over M workers, where we assume that the function U 0 (with its gradient) is available to every worker, whereas for each m = 1, . . . M the function U m (with its gradient) is only available to the mth worker.

Federated computation of the first event
The essential idea of Federated PDMC is that every worker m ∈ {1, . . . , M } proposes a switching time associated to their own rate function λ m . This proposed switching time is communicated (along with the proposed change in velocity) to the coordinating server which selects the minimum of the proposed switching times and the proposed switch. From this time and the new combination of position and velocity, the process is repeated. Under simple conditions this approach can be seen to have the correct stationary distribution.
We suppose that every worker has its own jump mechanism consisting of jump intensity (The deterministic dynamics φ are identical for all workers.) Furthermore every worker computes a new choice of velocity v m according to their individual jump distribution Q m , i.e. every worker Then every worker sends its proposal (τ m ,Ṽ (m) ) to the server. The server determines the minimum switching time and associated new velocity. Therefore the effective switching time for Federated PDMC is the first arrival time of an inhomogeneous Poisson process with rate The above procedure, described in detail in Algorithm 2, provides a genuinely federated algorithm since every machine only requires access to U m , without affecting the invariant probability distribution. Indeed it is straightforward to verify that the conditions of Section 2.3 are satisfied for (λ fed , Q fed ). Moreover any ergodicity properties of the non-federated PDMC algorithm carry over to the Federated PDMC algorithm, since the effective event rate λ fed of Federated PDMC is increased relative to non-federated PDMC, as discussed in Section 4.1. Therefore the obtained skeleton points may be used as discussed in Section 2.4.
Remark 3.1. Each worker only requires access to the component of the global variable which the function U m depends on. It is straightforward to reformulate Algorithm 2 so that each worker only acts on a set of local variables, which are subsequently mapped onto the global state X n by the central server, and vice versa.
Remark 3.2. There is flexibility in how to accommodate for the prior distribution: it can be absorbed in the function U 0 or distributed amongst the workers through the functions U m . Alternatively it could be handled by an artificial extra worker which simulates event times associated with the prior distribution. Also combinations of these approaches are possible.

Computational efficiency of Federated PDMC
Consider a Bayesian context in which we have N observations, distributed among M workers, each worker having access to a batch of size n m points, with N = M m=1 n m . In many relevant Locally: Every worker simulates τ m andṼ (m) such that Centrally: Set

5: end for
cases the computation of the proposed switching times τ m will be a bottleneck factor for the total computational effort, and it is reasonable to expect that this effort is linear in the size of the data. Since the workers operate in parallel, we find that the computational effort required is of order max m=1,...,M n m . As a special case, if all batch sizes are equal, i.e., n m = N/M for all m, the computational effort required per iteration is N/M . We see that the computational effort per proposed switch is reduced by a factor M relative to the case in which all N observations would be processed by a single worker.
This simple computation does not yet paint the full picture: Although the invariant distribution is not affected by employing Federated PDMC, the event rate is modified since the operation of taking the positive part a → (a) + occurs at each individual worker. We will consider the effect of this in some detail below.

Expected event rate for the exponential family
Consider the 'canonical rate' λ can , corresponding to a single machine generating the switches associated with the potential function U (x) = M m=1 U m (x). In the context of the BPS and Boomerang Sampler (as an example), We further investigate the expected switching rate for the federated one-dimensional Zig-Zag sampler for exponential family models indexed by a parameter x. We assume the data is generated from this model for a fixed (unknown) parameter value x 0 .
For simplicity we assume a flat (Lebesgue) prior measure. We show in the Appendix that in this situation the posterior expected event rate under the federated intensity and under the distribution of the data is magnified by a factor √ M , compared to the canonical rate: Remark 4.1. The larger (expected) effective switching rate results in more simulated switches per unit time interval, so that a larger computational effort is required to simulated such an interval. There is another aspect which affects computational efficiency: As the switching rate increases beyond the canonical rate, the process trajectories become more diffusive (see [BD17]), resulting in an increased Monte Carlo error per simulated unit time interval. The quantification of this error is beyond the scope of this work.
Remark 4.2. For an exponential family the analysis breaks down into easily manageable parts but it may well be possible to generalize these results beyond this setting. Also the scaling dependence on dimension in multivariate settings is left for further research.
Remark 4.3. The increase of the event rate λ fed compared to λ can can be reduced by taking a control variates approach, in a similar approach as discussed in [BFR19].

Differential privacy
A relevant aspect of federated inference is the amount of privacy achieved by taking a federated approach. In Federated PDMC, every worker only computes its proposed switching time along with the corresponding proposed change of velocity. This is a very limited amount of information and as such at an intuitive level we may be optimistic about the privacy achieved by the Federated PDMC approach. For a theoretical approach, we may employ the concept of differential privacy [DKM + 06]. In differential privacy the privacy of a stochastic algorithm is quantified as follows. For ε, δ > 0, an algorithm gives δ-approximate ε-indistinguishability or simpliy (ε, δ)-differential privacy if for outputs τ andτ corresponding to data sets differing at most one row, for any measurable set S. We will now investigate the differential privacy of Federated PDMC. We make the assumption that where a single change of the ith observation in batch m results in a change U

Differential privacy for communication of switching times
We suppress the dependence on m, so let λ and λ denote the switching intensities of a single machine with corresponding proposed switching times τ andτ . We have the following result.
The proof of Theorem 5.1 is provided in the Appendix. Theorem 5.1 establishes that in order to achieve a certain level of differential privacy, we may tune the refreshment rate of the algorithm to be at least ρ as specified by (7). This can be achieved by setting the switching rate for each machine m as where λ m (x, v) is any valid switching intensity and with ρ as desired. Indeed, the conditions of Section 2.3 remain to be satisfied if a constant 3 is added to a valid event intensity.

Experiments
In this section, we demonstrate how the the Federated Zig-Zag sampler could be used for a number of classical Bayesian computation benchmarks. In each case, we study the influence of the number of distributed workers on both the computational efficiency as well as the speed of convergence of the PDMP to the target equilibrium distribution. Further implementation details, including the derivation of dominating Poisson processes for computing the next switch time can be found in the Appendix. An additional example demonstrating the Federated Zig-Zag sampler on a spatial log-Gaussian Cox model is provided in the Appendix.

Multivariate Gaussian Distribution
In this scenario, we assume that N = 50 independent, identically distributed observations y 1 , . . . , y N , are made of a d = 10 dimensional multivariate Normal distribution with unobserved µ and covariance matrix of the form Σ = α 2 I d×d , where α = 1, for simplicity. We choose this simple example to be able to empirically validate the results of Section 4.1. Assuming a uniform improper prior on µ for simplicity, the goal is to characterise the posterior distribution of µ|(y 1 , . . . , y N ), where it is assumed that the N observations are distributed equally amongst M workers. We assume that y is generated with true mean µ 0 = (0.5, . . . , 0.5). The posterior distribution is then given by N (ŷ, Σ/N ), whereŷ = 1 N N i=1 y i . In Figure 1 we plot the estimated mean switching rate for 20 independent runs of the Federated Zig-Zag Sampler, as a function of the number of processors, and compare them against the theoretical prediction obtained in Section 4.1. The solid curve demonstrates the theoreticallyderived expected switching rate of the central server with leading order σ 0 √ N M , where σ 0 = αd. We can clearly see a good agreement between experiment and theory. To better understand the effect of the model federation on the convergence of the PDMP we estimate the effective sample size (ESS) of the process across multiple, independent runs. In Figure 2 we plot the average ESS per gradient evaluation of the full potential. In the fully sequential setting, we observe that the ESS per evaluation is decreasing as we add workers, due to the diffusion arising from increased switching induced by the Federated Zig-Zag algorithm, as discussed in Remark 4.1. However, if we take into account the fact that the gradient evaluations are in fact parallelised (so that effort is reduced proportionally to the number of workers), one observes a benefit in increasing the number of workers.

Logistic Regression
We now consider a Bayesian logistic regression problem. Given N observations {(ξ 1 , η 1 ), . . . , (ξ N , η N )}, where ξ i ∈ R d and η i ∈ {0, 1}, we postulate that η i ∼ Ber(p i ), such that logit(p i ) = ξ i x, for an unknown x ∈ R d . The first component of each ξ i is taken to be equal to one to allow for an intercept in the model. The posterior distribution for x given the observations is given by where π 0 is the prior on x assumed to be standard normal, independent Gaussian distribution on all the components. Implementation details of the implementation of the Zig-Zag sampler for this model can be found in the Appendix. To demonstrate the Federated Zig-Zag method we generate N = 1000 synthetic observations with d = 6, and distribute them over M workers.
To show that the target posterior distribution is well approximated, we compared the output of the Federated Zig-Zag with a large MCMC sample for the same posterior, generated using Hamiltonian Monte Carlo (HMC), implemented in Blackjax [LL20]. In Figure 3 we plot the marginal 1-Wasserstein distances between 20 independent runs of the Federated Zig-Zag scheme, run up to T = 100, and a reference sample from the global posterior distribution. It is clear that the overall error is small, and is not affected by the distribution of data across the workers.

Time Series Model
Consider the AR(1) model in which a time series (Y i ) satisfies the stochastic update rule Y i = c + xY i−1 + ε i , for i = 1, 2, . . . . Here c and x are unknown constants, and (ε k ) are i.i.d. random variables in R drawn from a distribution with density function g.
We assume that we make N independent observations of trajectories of the time-series, each observed at K points, i.e. we observe y (1) , . . . , y (N ) , where y (i) = (y We consider a robust inference setting, where we choose g to be the density of a heavy-tailed distribution. In this example, we choose g to be a Student-T distribution with ν degrees of freedom. Recall that the heaviness of the tails increases as ν → 0. Mathematical details on the implementation of this model using the Zig-Zag sampler are provided in the Appendix. To demonstrate the Federated Zig-Zag algorithm we suppose that N observations are evenly distributed amongst the M workers. In Figure 4 we plot the marginal Wasserstein distance between 20 independent runs of the Federated Zig Zag, output obtained after T = 100 process time units of simulation, and a reference MCMC sample obtained using HMC. We clearly see that the sampler is able to correctly approximate the correct posterior, and that this remains stable as the number of workers increases).

Discussion
In this work we have introduced a generic method for Bayesian computation aimed at federated or distributed multiplicative compositions of local posterior distributions.
Our approach hinges strongly on the notion of Piecewise Deterministic Monte Carlo, a field in Bayesian computation that is currently very much under development. An intrinsic challenge of PDMC is the simulation of event times which in practice relies upon the use of a-priori bounds on the (local) gradients of the log posterior. Other numerical approaches to the computation of event times exist; see e.g. [CSR22, PCP + 20]. The analysis and numerics presented in this paper have focussed on a federated version of the Zig-Zag Sampler. Similar generalizations are possible for other PDMP based samplers with distinct advantages and disadvantages which we hope to study in future work.
Further research on the numerical efficiency of the distributed simulation of switching times in general multivariate settings for different PDMC algorithms is necessary. Such an understanding would be instrumental in having rules of thumb for designing optimal architectures for federated or parallel Bayesian inference, such as the optimizing the number of parallel workers. Efficiency gains may be achieved using control variates in similar spirit to [BFR19].

A Event time simulation
An important practical aspect of PDMP simulation is drawing the random times τ satisfying Equation (2) of the manuscript. We wish to simulate τ such that where (x, v) is any initial position of the trajectory.
To achieve this using the method of Poisson thinning, we assume that for every (x, v) there is a bounding function t → λ(t; x, v) such that λ(φ(t; x, v))) ≤ λ(t; x, v) for all t ≥ 0. We furthermore assume that the functions λ(·, x, v) are chosen in such a way that there is an explicit formula for the inverse function In words, σ is distributed according to the first jump time of a inhomogeneous Poisson process with rate function (λ(t; x, v)) t≥0 .
In order to obtain a switching time with the desired distribution, we follow an iterative procedure. We sample a proposed switching time τ satisfying P(τ ≥ t) = exp(− t 0 λ(s; x, v) ds), which we accept as true switching time with probability λ(φ(τ ; x, v))/λ(τ ; x, v). If we do not accept the proposed switching time, we increase the time variable by τ and repeat with new starting point φ(τ ; x, v). The full procedure is given in Algorithm 3. Set S = true with prob. λ(φ(σ; x, v))/λ(σ; x, v)

6:
Set τ = τ + σ and (x, v) = φ(σ; x, v) 7: end while In many settings we have that λ(x, v) = ( v, ∇U (x) ) + , (or its one-dimensional variant, , where U has a bounded Hessian. If we also assume linear trajectories as used in the Zig-Zag Sampler and the Bouncy Particle Sampler, using Lipschitz continuity of a → (a) + , we find using the mean value theorem that where · denotes the Euclidean norm, and we can use any vector norm · p , and associated induced matrix norm · p , for p ∈ [1, ∞], in the estimate for λ i .
In the context of Federated PDMC as discussed in Section 3 the procedure outlined in Algorithm 3 can be used for the simulation of the switching times τ m for each machine: in this case we just replace the global switching rate λ by the switching rates λ m of the individual machines and make sure we find a suitable upper bound λ m with accompanying inverse H m .

B Expected switching rate for exponential families
In this section we consider the expected switching rate for the federated one-dimensional Zig-Zag sampler for data (y i ) generated from an exponential family model with parameter x, f (y; x) = h(y) exp (c(x) + xt(y)) . (9) We consider the situation where the data (y 1 , . . . , y N ) is generated from (9) for a fixed 'true' parameter x 0 . We assume that the data (y 1 , . . . , y N ) is partitioned into M batches m = 1, . . . , M , where the mth batch consists of n m elements, denoted by y m,i , i = 1, . . . , n m . For simplicity we assume a flat (Lebesgue) prior measure.
t(y m,i ), and using that We consider the Zig-Zag Sampler in R. Write µ(dx, dv) ∝ exp(−U (x)) dx⊗Uniform({−1, +1}). For the federated learning switching intensity we estimate, using the 1-Lipschitz property of a → (a) + , We analyze the two terms separately. First we consider the posterior expectation of the canonical switching rate, E µ λ can (x, v). We have E π U (x) = 0. Assuming posterior contraction (see the Bernstein-von Mises theorem, [vdV98, Section 10.2]), we have asymptotically that wherex =x(y 1 , . . . , y N ) denotes the maximum likelihood estimator for x and I(x) denotes the Fisher information associated with the parametric model f (·; x). (The expression (11) should be interpreted in an appropriate asymptotic sense as in [vdV98].) By the Delta method, asymptotically, Next we consider the effective excess switching rate, bounded by M m=1 |Z m |. We have Recall that y i are i.i.d. according to f (·; x 0 ) for a fixed parameter x 0 , and denote σ 2 0 for the variance of t(y i ). We see that the mean of Z m is zero, and its variance is Therefore we have Using Jensen's inequality, Combining the terms in (10) we find that According to this analysis the switching rate of the Federated PDMC Sampler (Algorithm 2 in the manuscript) is increased by a term of O( √ M N ) relative to the canonical rate of the standard PDMC Sampler (Algorithm 1 in the manuscript), which is O( √ N ).

C Differential privacy
For simplicity we write the proof of Theorem 5.1 in the manuscript in terms of time dependent switching rates λ(t). These may be interpreted as λ(t) = λ(φ(t; x, v)) in the context of Federated PDMC.
Proof. We have from which it is straightforward to obtain the stated result.
We see that under reasonable conditions we have ε-differential privacy, provided that we restrict the maximal switching time to a finite time interval. As an extension, we may obtain an upper bound on the probability to have a switching time larger than this time, in order to obtain (ε, δ)-differential privacy.
Proof. This follows since We are now ready to provide the proof of Theorem 5.1.
Proof. Using Lemma C.2 it suffices to show that, for some t 0 > 0, where f τ denotes the density function of τ . Let t 0 = ε−log γ K with γ = K ρ + 1. Indeed we have Furthermore we have that and using Lemma C.1 we find that the ratio of densities f τ /fτ is bounded by exp(ε) for t ≤ t 0 .
Finally, for fixed ε > 0 and δ > 0, if we take ρ satisfying the indicated inequality, then we must verify that ε > log 1 + K ρ and that Indeed,

D Increased Privacy via Dynamic Prior Switching
Suppose the communication between the server and one the nodes is compromised by an attacker. Then over a long period of time, it would be theoretically possible for the attacker to reconstruct the global potential function by observing the switches over long periods of time. To mitigate this we propose a heuristic strategy which obfuscates the individual worker contributions to the likelihood by weighting them with a random piece of the prior, which evolves dynamically. To be more specific, suppose we are targeting the following global posterior: where U m is the potential contribution for the m th worker, and π 0 ∝ e −U0 is a prior. The proposed strategy is to distribute the prior across the workers so that the local potential for the m th worker becomes where the weights α 1 , . . . , α M are randomly chosen scalars such that M m=1 α m = 1. Clearly, exp(− m U m (x)) ∝ π(x). To dynamically change the weights, we introduce a constant redistribution rate λ redist and assume that prior re-distributions (i.e. resampling of the α's) occurs at discrete times determined by an independent Poisson process with constant rate λ redist . The new algorithm is expressed in Algorithm 4. Note that the worker routine remains unchanged from the standard Federated Zig-Zag algorithm.
The dynamic re-distribution of the prior introduces time-inhomogeniety into the process, and it does not formally follow from previous results that the resulting process has the correct unique stationary distribution, though it is intuitively clear that would be the case. We leave the analysis of this new PDMP variant as a subject for future work. To demonstrate that the redistribution process does not affect convergence to equilibrium, we repeat the experiment from the Logistic Regression example from Section 6.2 of the main paper, under the same conditions, with λ redist = 0.1. We plot the marginal Wasserstein distances against a HMC reference sample in Figure 5. Compared to the results for the Federated Zig-Zag in Figure 3 of the main text we observe that the introduction of the prior re-distribution process does not adversely affect the convergence behaviour of the continuous-time process.

5:
if τ redist < τ m0 then 6: Resample α 1 , . . . , α M , such that M m=1 α m = 1. To demonstrate a more complex example of federation across multiple workers, we consider a Log-Gaussian Cox model, similar to the one considered in [Gal16] and [WR20] for the Zig-Zag and Coordinate samplers, respectively. We assume that the observations y = {y ij } d i,j=1 are Poisson distributed and conditionally independent given the intensity λ = {λ ij } d i,j=1 where λ ij = exp(x ij ). The latent process x = {x ij } d i,j=1 is a Gaussian process defined on the d × d grid with mean zero and precision matrix P uv = β(δ u,v − αA u,v ), where A = (A u,v ) is the adjacency matrix of the grid and where u = (i, j) and v = (i , j ) are grid coordinates of the nodes. The posterior distribution for x given the observations y is given by We assume that the nodes are distributed spatially across M = 4 workers, see Figure 6. Let V 1 , . . . , V M be the nodes assigned to each worker. The m th worker will target local potential: The prior term is handled by the central server which has a potential of the form Note that the m-th worker only needs to observe the latent variables x ij for (i, j) ∈ V m which are relevant to its observations. The interactions between the latent variables through the Gaussian process prior are handled entirely the central server. We choose α = 0.1 and β = 1. For clarity of presentation we demonstrate the scheme for d = 4, i.e. a grid with 16 nodes, noting that our experiments run fine on much larger grids. We run the Federated ZigZag sampler for the 4 workers until time T = 100.0. We extract samples from the continuous time process by extracting samples from the resulting trace at 10 −3 time-steps obtaining approximately 10 5 samples. In Figure 7 we plot the marginal posteriors for the output of the Federated ZigZag process, compared against a reference sample generated by running HMC on the same target distribution. We see excellent agreement between the two sets of distributions.

F.1 Multivariate Gaussian Distribution
Consider a multivariate Gaussian distribution N (µ, Σ) where µ ∈ R d and Σ is a symmetric, positive definite d × d covariance matrix. Then the associated potential is given by where P = Σ −1 , and the intensity functions for the Zig-Zag process take the form so that Given that this intensity is of the form max(0, b + at) the distribution of the next event can be sampled directly by computing the associated inverse function H(u; x, v) in (8) exactly, so that for V ∈ U [0, 1], the random variable H(− log V ; x, v) will be distributed according to the first jump time of the Poisson process with rate function max(0, b + at).

F.2 Logistic Regression
We assume that worker m, with m ∈ {1, . . . , M }, has access to independent observations (ξ All observations are mutually independent, also between different workers. Here x ∈ R d is the parameter which we wish to infer. We assume a prior distribution π 0 (x) ∝ exp(−U 0 (x)) over the unknown parameter x. The posterior distribution is then specified as It can be written as π(x) ∝ exp(−U (x)), where In particular, the Hessian of U m is bounded for each m provided that the Hessian of U 0 is bounded.
By the chain rule, We may therefore obtain a bound on the Hessian norm as sup x,c (1 + y (i) k−1 ) 2 .

F.4 Log Gaussian Cox Model
The potential function for a single worker is given by U m (x) = (i,j)∈Vm (−y ij x ij + exp(x ij )), so that ∇ xij U m (x) = −y ij + exp(x ij ), for (i, j) ∈ V m .
For simplicity we flatten the index, so that k = d(i − 1) + j, and re-index x k = x (i,j) , so that ∇ x k U m (x) = −y k + exp(x k ).
The associated switching intensity is given by We can bound this above as follows: The last upper bound provides a means of exactly simulating a dominating inhomogeneous Poisson process which can be subsequently thinned to simulate the local worker Zig-Zag process for the Log-Gaussian Cox model. Indeed, we can perform the following steps to simulate the next switching local Zig-Zag sampler event from (x, v). 5. Set x = x + τ m0 v.