{'fraction_non_alphanumeric': 0.0423770898857805, 'fraction_numerical': 0.01003559013408376, 'mean_word_length': 4.551866743251005, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Divide-and-Conquer Monte Carlo Fusion

Combining several (sample approximations of) distributions, which we term sub-posteriors, into a single distribution proportional to their product, is a common challenge. For instance, in distributed `big data' problems, or when working under multi-party privacy constraints. Many existing approaches resort to approximating the individual sub-posteriors for practical necessity, then representing the resulting approximate posterior. The quality of the posterior approximation for these approaches is poor when the sub-posteriors fall out-with a narrow range of distributional form. Recently, a Fusion approach has been proposed which finds a direct and exact Monte Carlo approximation of the posterior (as opposed to the sub-posteriors), circumventing the drawbacks of approximate approaches. Unfortunately, existing Fusion approaches have a number of computational limitations, particularly when unifying a large number of sub-posteriors. In this paper, we generalise the theory underpinning existing Fusion approaches, and embed the resulting methodology within a recursive divide-and-conquer sequential Monte Carlo paradigm. This ultimately leads to a competitive Fusion approach, which is robust to increasing numbers of sub-posteriors.


Introduction
In this paper, we are interested in the following d-dimensional (product-pooled) target density (which we term the fusion density), where x ∈ R d , f c (x) for c ∈ {1, . . . , C} represent the individual densities which we wish to unify (termed sub-posteriors), and C represents the total number of distributions to be unified. For convenience, we assume that we have access to independent realisations from each sub-posterior, and that it is possible to evaluate each sub-posterior point-wise. As discussed later, neither of these assumptions form limiting factors for the methodology we develop.
The need to unify several (sample approximations of) distributions, over a common parameter space into a single sample approximation of the distribution in the manner of (1) is surprisingly common. For instance, it is a problem which classically arises in expert elicitation [1,2,13] and meta-analysis [12]. However, it has proven to be challenging methodologically in a number of modern settings due to problem specific constraints. These include when dealing with the privacy constraints of the individual sources [32], in cases where the sheer number of sources is overwhelming, or if the networking constraints of the sources are truly distributed [28]. This in turn has motivated a range of problem specific and pragmatic approximations. These approximations are invariably distributional, and imposed at the level of the individual source. Such approximations limit the applicability of methodological approaches to particular settings, and outside those settings the unified results can be poorly understood, and even misleading. In this paper, we instead focus on developing methodology for an exact Monte Carlo approximation of the unified distribution, one which provides robust inference in a wide range of practical problems, and yet is amenable to use alongside any problem specific constraints.
The majority of the recent methodological developments for representing or sampling from (1) have been focused on tackling distributed 'big data' problems [see for instance 28,24,30,22,29,25]. In this setting, due to the sheer size of the data, the data is split across a number of cores (say C cores), inference is separately conducted on each core (often using MCMC), and then the respective methodologies attempt to unify the sample approximations of the distribution (as per (1), and typically using a convenient approximation). In this paper, we will compare our methodology with a number of the most popular approaches, and so will briefly describe these here. The Consensus Monte Carlo (CMC) approach of Scott et al. [28] represents (1) by means of a weighted average of sub-posterior samples. It can be shown that CMC is exact when each sub-posterior is Gaussian, and can be useful in settings where each sub-posterior is approximately Gaussian, which is often the case in big data settings [20]. However, it has been shown to exhibit large bias in other settings [30]. Neiswanger et al. [24] suggest a strategy (which we term the Kernel Density Estimate Monte Carlo (KDEMC) approach) based on using a kernel density estimate to approximate the sub-posterior densities, and in effect approximating (1) by implicitly sampling from the product of non-parametric density estimates. Finally, the Weierstrass sampler of Wang and Dunson [30] provides an alternative method for approximating (1) by means of using the product of Weierstrass transforms for each sub-posterior.
The Fusion approach to dealing with (1) constructs a sample approximation of f itself, rather than seeking to obtain an ad hoc approximation to f by combining sample approximations of the sub-posteriors, f 1 , . . . , f C . The Monte Carlo Fusion approach of Dai et al. [7] is a rejection sampling approach to sample f by means of sampling from the individual sub-posteriors and a density on an extended space, and so returns iid draws from f . The Bayesian Fusion approach of Dai et al. [8] is an alternative sequential Monte Carlo (SMC) approach to mitigate some of the computational bottlenecks of Dai et al. [7]. Both existing Fusion approaches have key limitations in practice. In particular, and as detailed more fully in Section 4: unifying subposteriors which exhibit strong correlation structure, or sub-posteriors which conflict with one another; and robustness with unifying increasing numbers of sub-posteriors.
In this paper, we begin by reformulating the theory underpinning the existing Fusion approaches of Dai et al. [7] and Dai et al. [8], and introduce a Generalised Monte Carlo Fusion (GMCF) algorithm which is based upon importance sampling (as opposed to rejection sampling) and exploits readily available global information about each sub-posterior. The resulting GMCF methodology has far greater robustness to increasing sub-posterior correlation. We next embed our GMCF methodology within a divide-and-conquer paradigm by combining the sub-posteriors in stages to recover the correct fusion density f , in an approach we term Divide-and-Conquer Monte Carlo Fusion (D&C-MCF). We do this by appealing to the Divide-and-Conquer Sequential Monte Carlo (D&C-SMC) approach of Lindsten et al. [21], and so as a consequence we are able to take advantage of much of what is known about SMC and this particular variant thereof. The D&C-MCF approach addresses both the robustness of Fusion approaches with respect to increasing numbers of sub-posteriors and also by means of tempering methodology, provides a possible resolution to conflicting sub-posteriors. This is shown empirically by means of both illustrative and realistic benchmarks.
The remainder of the paper is organised as follows: in Section 2 and Algorithm 1, we present our Generalised Monte Carlo Fusion approach, allowing us to incorporate within existing Fusion approaches additional information about sub-posteriors, and so greatly improve their efficiency.
In Section 3, we frame our Fusion approach within a divide-and-conquer paradigm, and by introducing a divide-and-conquer sequential Monte Carlo algorithm, address the robustness of Fusion with increasing numbers of sub-posteriors; our tree diagrams for the fusion problem present how the sub-posteriors can be combined in stages, in a flexible problem-specific manner.
We present illustrative examples to study the robustness of our methodology applied to various scenarios in Section 4. In Section 5, we present some applications of our methodology and study the performance of our approach in comparison to competing methodologies. We conclude by outlining a variety of ways the Fusion approaches we introduce could be embedded within other Fusion approaches (for instance, Bayesian Fusion [8]), or separately extended and used in other settings. All technical proofs are collated in the appendices.
Statistical computations for this paper were written in R, C ++ and Rcpp [9]. All code developed in the implementations presented, and for producing numerical results, can be found at https://github.com/rchan26/hierarchicalFusion.

A generalisation of Monte Carlo Fusion
The key observation in the Fusion approach is that if we wish to sample from the density f (C) ∝ c∈C f c , where C is an index set representing the sub-posteriors we want to unify, then this is simply a marginal of a density on an extended state space (g C given in Proposition 2.1). For the purposes of simplifying the notation in Proposition 2.1, we denote by x (C) ∈ R |C|×d a vector composed of x (c) c∈C ∈ R d (in particular, we have x (C) := x (c 1 ) , . . . , x (c |C| ) , with c i denoting the i th element of the index set C). Recall, from the introduction we assume that we can evaluate each sub-posterior point-wise. Although not required for Proposition 2.1, in the subsequent methodology we develop (in common with Dai et al. [7]) we further assume that for c ∈ {1, . . . , C}, f c is nowhere zero and everywhere differentiable, and that we can compute A c (x) := log f c (x), ∇A c (x), and div∇A c (x) point-wise (where ∇ is the gradient operator and div is the divergence operator).
Proposition 2.1. Suppose that p c is the transition density of a Markov chain on Ê d with a stationary probability density proportional to f 2 c . Then the (|C| + 1)d-dimensional probability density proportional to the integrable function admits marginal density Dai et al. [7] exploited Proposition 2.1 by noting that if the index set C := {1, . . . , C}, then we recover the target fusion density f (as given in (1)). In the theory we develop in this section, we consider the abstraction to more general index sets, as this facilitates the divide-and-conquer approach we introduce in Section 3.
Since g C will not typically be accessible directly, Dai et al. [7] proposed sampling from g C by constructing a suitable (|C| + 1)d-dimensional proposal density (say, h C ) for use within a rejection sampling algorithm [7, Algorithm 1], and then simply retaining the y (C) marginal of any accepted draw as a realisation of f (C) . Dai et al. [7] showed that if p c in Proposition 2.1 was chosen to be the transition density of a constant volatility Langevin diffusion at time T with invariant measure f 2 c for each c ∈ C respectively, then for a (easily accessible) proposal h C constructed by sampling a single draw from each sub-posterior (x (c) ∼ f c for c ∈ C), and then a single Gaussian random variable parameterised by the sub-posterior realisations (corresponding to the y (C) -marginal), the acceptance probability was readily computable.
As noted in the introduction, one of the principle shortcomings of the rejection sampling Monte Carlo Fusion (MCF) approach of Dai et al. [7] is its lack of robustness with increasing subposterior correlation. In this section, we develop theory and methodology to improve upon the approach of Dai et al. [7] (even in the setting where C := {1, . . . , C}), by providing a principled way to weight the contribution from each sub-posterior in the construction of the proposal h C , and then embed the approach in more robust Monte Carlo methodologies (such as importance sampling and sequential Monte Carlo). Having the flexibility to weight the contribution from each sub-posterior allows us to incorporate within the algorithm global information for each sub-posterior which will often be readily obtainable, such as information about sub-posterior means or covariance structures (x c andΣ c for c ∈ C respectively).

Theory
Although Monte Carlo Fusion [7] is a rejection sampler and outputs independent realisations from (1), the computational efficiency of the approach is limited by the quality of the proposal used. The construction of the proposal suggested by Dai et al. [7] closely matches the extended target when each sub-posterior has a covariance structure which is approximately a scaled identity (i.e.Σ c ≈ (1/T )I d for c ∈ C and user chosen constant T > 0). However, there is a lack of robustness in the approach (in the sense that the acceptance probabilities rapidly degrade) when the sub-posteriors have non-identity covariance structure. In this section, our aim is to improve the quality of the proposal to address this lack of robustness. As we assume we are able to simulate directly from each sub-posterior x (c) ∼ f c for c ∈ C, it is possible to alter the proposal of Dai et al. [7] via the manner in which the |C| sub-posterior realisations paramaterise the Gaussian random variable corresponding to the proposed y (C)marginal. A natural choice for the proposal parameterisation would be one which incorporated global information for each sub-posterior, for example using estimated sub-posterior covariance structures, and estimated sub-posterior means. In particular, we propose the proposal density for the extended fusion target (2) to be proportional to the following function (where we let Λ c be the user-specified matrices associated with sub-posterior f c for c ∈ C): wherex Simulating from the proposal h C in (3) is possible by first simulating a single draw from each sub-posterior (which are assumed to be accessible) to obtain {x (c) } c∈C . These draws are then used to computex, and then simulate y (C) ∼ N x (C) , T Λ C . This proposal for y (C) is a Gaussian perturbation of the approach suggested in Consensus Monte Carlo (CMC) method [28] for combining sub-posterior draws (namelyx (C) ), on a scale commensurate with that of the posterior (modulated by a user-specified parameter T ).
Now considering the choice of the |C| stochastic processes on Ê d with stationary densities proportional to f 2 c for c ∈ C (as required in Proposition 2.1), we similarly want to exploit the availability of global sub-posterior information. As such, we choose p c in Proposition 2.1 to be the transition density over [0, T ] of the following stochastic differential equation: where in this case Λ c can be interpreted as a preconditioning matrix with Λ 1/2 c being the (positive semi-definite) square root of Λ c where Λ Brownian motion. Note that for the purposes of our numerical simulations we used the Schur decomposition.
Having constructed a more informative proposal h C , we now require an expression for the ratio of the target to proposal densities. Proposition 2.2. Letting p c y (C) |x (c) be the transition density of the diffusion given in (5) we have where and where Tr(·) denotes the trace of a matrix, and W Λc denotes the law of a Brownian bridge {X (c) T := y (C) and covariance matrix Λ c .
Proof. See Appendix A.

Methodology
It would be possible to develop Monte Carlo methodology based directly upon Propositions 2.1-2.2 provided we could simulate from h C and compute the quantities ρ 0 and ρ 1 . Alternatively, viewed from an importance sampling perspective, Proposition 2.2, allows us to use ρ 0 and ρ 1 as importance weights. The quantity ρ 0 in (7) in essence penalises proposals based on the distance the weighted average of the sub-posterior draws are from one another, whereas ρ 1 penalises how far the proposal y (C) is from each sub-posterior draw under the corresponding sub-posterior.
Considering the unnormalised importance weight, the computation of ρ 0 is direct, whereas computing ρ 1 in (8) is not direct as it involves the evaluation of path integrals of functionals of Brownian bridges. This is the key complication in this methodology.
For the purposes of the methodology we will subsequently develop, it is sufficient to find a non-negative unbiased estimator of ρ 1 with finite variance which can be obtained with finite cost. To do so, we require for a given sample path X (c) [0,T ] ∼ W Λc that we have absolute bounds on φ c for each c ∈ C. In general, it is not possible to find global bounds for φ c , and so following the approach of Beskos et al. [3] and Pollock et al. [26], we simulate Brownian bridges in such a way that we can determine a compact region which almost surely constrains a given path (the details of which appear in Pollock et al. [26,Sections 8.1 and 8.5]), which enables us to instead find local bounds on φ c . In particular, we let . We note that it is possible to partition the sample space into disjoint sets and simulate from associated distribution function (without having to sample the underlying path), R c ∼ R c conditional on the user specified partitioning of the space. Although it is possible to find tight local bounds for φ c in a problem specific manner by exploiting specific structure, it is helpful for practitioners to note that it is possible to find generic (less tight) bounds (which we denote by L (c) where d denotes the dimension of x, · is the Euclidean norm,x (c) is a user-specified point central to R c , and where P Λc n is a quantity such that with γ denoting the spectral radius.
Once local bounds for φ c are obtained, it is possible to be unbiasedly estimate ρ 1 using auxiliary diffusion bridge path-space samplers developed in (for instance) Beskos et al. [4,3], Fearnhead et al. [10], Pollock et al. [26]. In particular, where the subscript for each expectation denotes the law with which we are taking the expectation,R denotes the law of {R c ∼ R c : c ∈ C},W denotes the law of the |C| Brownian bridges (conditional on the |C| distinct starting points, but with common endpoint y) Corollary 2.1.ρ 1 is a positive unbiased estimator of ρ 1 .
Proof. X which bound φ c as per Proposition 2.3, there are two natural choices of unbiased estimator that could be used by making particular choices for the distribution of the discrete random variable used to simulate κ 1 , . . . , κ C . We denote theseρ X )T for c ∈ C leads to the following unbiased estimator of ρ 1 : leads to the following unbiased estimator of ρ 1 : The estimatorsρ which are positive with finite variance. In addition,ρ (a) for all x and c ∈ C.
In developing a rejection sampling approach to sampling from f (C) by means of the extended density g C with our proposal h C (as suggested by Proposition 2.1), it is natural to select the estimatorρ we recover the MCF approach of Dai et al. [7].
Our approach embeds the Fusion methodology within the Divide-and-Conquer Sequential Monte Carlo (D&C-SMC) approach of Lindsten et al. [21] in order to address the robustness of Fusion with increasing numbers of sub-posteriors. In this setting, we are not restricted to obtaining independent realisations of f (C) in (1), and so we can instead employ the estimatorρ 1 (as we no longer require our estimator to be bounded above by a constant), which has better asymptotic properties [10]. The theoretical properties of D&C-SMC are increasingly-well characterized and include a strong law of large numbers, finite sample L p errors bounds as well as a √ N -central limit theorem under mild conditions (see Kuntz et al. [19]).
A (self-normalised) importance sampling approach would proceed as follows: approximate the extended (|C| + 1)d-dimensional density g C in Proposition 2.1 by sampling from the (|C| + 1)ddimensional proposal density h C (i.e. ( x (C) , y (C) ) ∼ h C ) as given in (3) and (4). Thus far, we have assumed we have access to iid draws from each sub-posterior, but this is not necessary for our methodology. In particular, in many settings (including that in Section 3) we will have instead have access to importance weighted draws from each sub-posterior (w (c) for c ∈ C for instance, although in the simplest setting w (c) ∝ 1 for c ∈ C). As such, the proposals are then assigned an unnormalised importance weight w ′ ( x (C) , y (C) ) := ( c∈C w (c) ) · (ρ 0 ·ρ 1 )( x (C) , y (C) ).
In practice, and to better relate to our methodology developed later in Section 3, we adopt a variant of the importance sampling approach described above. We term the following approach Generalised Monte Carlo Fusion (GMCF), and summarise it in Algorithm 1. We assume for simplicity in presenting our algorithm that we have access to M importance weighted samples from each sub-posterior (each of which form a Monte Carlo representation for their respective sub-posterior). In addition, we note that the marginal importance sampling weight ρ 0 (in (7)) only depends upon the sub-posterior realisations (and not y (C) ). To exploit this fact, we compose M partial proposals by pairing the draws from each sub-posterior { x We then sample with replacement N times from this collection of M draws with associated partial weights w j · ρ 0 ( x (C) 0,j ). For each of the N partial proposals, we then complete the proposal by simulating a corresponding endpoint y (C) i . Taking account of our earlier sampling, each of the N particles will be given (an appropriately normalised) weight proportional to ρ 1 ( x i ). By retaining for each of the N weighted particles the marginal for y (C) (i.e. {y , we simply have an approximation to the desired density, for c ∈ C, the user-specified matrices, {Λ c : c ∈ C}, the number of particles required, N, and time horizon T > 0. . . , M} as per (7).
Of course, in general in the Input step of Algorithm 1, we may have access to different numbers of samples from each sub-posterior: say M c importance weighted samples for sub-posterior f c (for c ∈ C). In order to compose our M partial proposals in Step 1, there are a number of approaches we could take. As presented above, if M c = M for c ∈ C, we simply pair the sub-posterior draws index-wise. This is a basic merging strategy of the sub-posterior realisations and has the advantage that it can be implemented in O(M ) cost (and if M c = M for every c ∈ C one could simply sub-sample to obtain a common number of samples from each sub-posterior). However, as noted in Lindsten et al. [21], while this approach has a low computational cost, it can lead to high variance when the product c∈C f c (x (c) ) differs substantially from the corresponding marginal of f (C) -which one might expect to be the case in our setting when the sub-posteriors disagree.
We found this simple approach more than adequate in our simulations, but there are more sophisticated options available should they be required in still more challenging settings. In particular, as described in Lindsten et al. [21,Section 4.1], at the expense of a computational cost O( c∈C M c ), one could instead compose all possible permutations of the samples from each sub-posterior before weighting and then resampling to reduce the number of points in the approximation back to a pre-specified number, arriving at a better approximation at a greater cost. They termed this approach "mixture resampling" and also detailed a "lightweight mixture resampling" approach in which more than one permutation, but not all possible permutations, are used and found it to work well; as noted by Kuntz et al. [18] such a strategy can be connected directly with the theory of incomplete U -statistics and consequently one might hope to realise much of the benefit of mixture resampling at a much reduced cost [see, e.g. 17]. Although Step 2a corresponds to a multinomial resampling of the partially composed proposals, an approach we followed in the interest of simplicity, we can of course use other resampling methods and might expect better performance by choosing a lower-variance approach [see 14, for a recent investigation of the properties of many such schemes].

A divide-and-conquer approach to Fusion
A key drawback of the Monte Carlo Fusion approach of Dai et al. [7] is that it lacks robustness with increasing number of sub-posteriors, |C|. This is unsurprising as the extended target and proposal densities (g C and h C ) of Proposition 2.1 are (|C| + 1)d-dimensional, and these become increasingly mismatched with increasing dimension. In particular, as a consequence of the definition of ρ 1 in (8) of Proposition 2.2, the acceptance probability of any rejection-based scheme will decrease geometrically with increasing |C|.
As presented both in Dai et al. [7] and Section 2, Fusion is an example of a fork-and-join approach which unifies all of the sub-posteriors in a single step. In particular, using the notation of Section 2, all elements of the index set C := {1, . . . , C} are combined simultaneously. This is illustrated in the tree diagram of Figure 1, where the leaves of the tree represent the available sub-posterior densities, the directed edges are used to illustrate the computational flow of Monte Carlo Fusion, and the root vertex of the tree is the desired fusion density, f (as given in (1)). As the goal of the methodology is to approximate f in (1), one could envision a recursive divideand-conquer approach in which the sub-posteriors are combined in stages to recover f . There are a number of possible orderings in which we could combine sub-posteriors, and so we represent these orderings in tree diagrams, and term these hierarchies (see Figure 2). For instance as illustrated in Figure 2a, one approach would be to combine two sub-posteriors at a time (we term this a balanced-binary tree approach). In Figure 2a, the intermediate vertices represent intermediate (auxiliary) densities up to proportionality. The approximation of the distribution associated with any non-leaf vertex is obtained by an application of Fusion methodology to the densities of the children of that vertex. A balanced-binary tree approach is perhaps the most natural way to combine sub-posteriors in a truly distributed setting (where the simulation of each sub-posterior has been conducted separately, and so the inferences we wish to combine are distributed). Another approach one might employ is given in Figure 2b, whereby sub-posteriors are fused one at a time (and so we term this the progressive tree approach). This is perhaps the most natural approach in an online setting. We focus on fusion methodology applied to these two natural hierarchies for the remainder of this paper, although other hierarchies are certainly possible within our framework, and there is no limitation in unifying more than two vertices at any level of a tree (as suggested by both Section 2 and Figure 1).
From this recursive perspective, sample approximations of auxiliary densities obtained at one level of any tree are themselves treated as sub-posteriors at the next level up. As such, one can iteratively apply the Fusion methodology of Section 2, working through the levels of the tree from the leaves to the root, using at each stage the output of one step as the input for the subsequent step. An advantage of our divide-and-conquer approach is that as fewer sub-posteriors are combined at each stage, we avoid (at each stage) the rapidly diminishing and variable importance weights. By utilising the importance sampling approach to Fusion we developed in Section 2.2, we can embed Fusion within Sequential Monte Carlo to address the robustness of Fusion with increasing |C|, albeit this being a trade-off with the cost of the repeated application of the methodology.
A divide-and-conquer variant of Sequential Monte Carlo (D&C-SMC) was recently introduced in Lindsten et al. [21]. D&C-SMC generalises the classical SMC framework from sequences (or chains) to trees, such as those in Figures 1 and 2 For the purposes of utilising the methodology developed in Section 2.2, we define the following notation for non-leaf vertices v / ∈ Leaf(T): let C v := ∪ u∈Ch(v) C u denote the index set representing the sub-posteriors that we want to unify for vertex v / ∈ Leaf(T). In addition, to simplify the notation and avoid an unnecessary level of subscripts, we index densities and other quantities by v rather than C v when it is clear what is intended. In particular, let Λ v := Λ Cv ,x (v) :=x (Cv ) as per (4) The extended target and proposal densities for vertex v / ∈ Leaf(T) are denoted g v := g Cv and h v := h Cv , respectively. Lastly, the importance sampling weights for v / ∈ Leaf(T) are given by ρ To describe the Divide-and-Conquer Monte Carlo Fusion approach, we specify an algorithm that is carried out at each vertex v ∈ V which leads to a recursive procedure; an initial call to d&c.fusion(Root(V), . . .) carries out the overall approach. For v ∈ V, we define a procedure (as given in Algorithm 2), which returns a weighted particle set { x denotes the normalised importance weight of particle i at vertex v ∈ V. From this particle set, we can take the marginal weighted samples for y (v) to approximate the fusion density f v ∝ posteriors. If independent sampling is not feasible, one could use MCMC to obtain unweighted sample approximations at the leaves formal arguments under appropriate regularity conditions could in principle follow an approach analogous to that in [11]. If v is a non-leaf vertex, we simply call Algorithm 1 by inputting the importance weighted samples {y for u ∈ Ch(v). As in standard SMC although the auxiliary distributions are defined on larger spaces, we do not need to retain sampled values which are not subsequently used; to obtain a more computationally manageable algorithm, we can choose to retain only the final parameter space marginal at each vertex (i.e. only returning the weighted particle set {y ) since the computation of the importance weights ρ Note when calling Algorithm 1 (in Step 3 of Algorithm 2), we input N importance samples from each sub-posterior (i.e. M u = M for u ∈ Ch(v)). As noted in Section 2.2, for simplicity, we use these to obtain partial proposals by pairing the sub-posterior realisations index-wise. If we use this simple merging scheme then M = N , and so resampling in Algorithm 1 Step 2a might not be necessary. As standard within the SMC literature, in practice we choose to re-sample only when we observe weight degeneracy [16]. We monitor weight degeneracy by computing the effective sample size (ESS := ( N j=1 w 2 j ) −1 ) of the particle set, and if it falls below some userspecified threshold then resampling. There are two steps where we can compute the ESS (and if necessary resample): (1) after computing the weights for the partially composed proposals, . . , N } (in Step 1); and (2) after computing the weights ρ

Illustrative comparisons with Monte Carlo Fusion
As discussed in the introduction, although direct Monte Carlo approximations to (1) (without approximating the individual sub-posteriors) offered some compelling advantages, such existing exact methods lack robustness in various key practical settings. In this section, we re-visit these settings, and with the aid of illustrative examples, show that the Generalised Monte Carlo Fusion and Divide-and-Conquer Monte Carlo Fusion approaches we introduced in Sections 2 and 3 respectively, address these key bottlenecks when contrasted with the earlier Monte Carlo Fusion of Dai et al. [7]. In particular, in Section 4.1 we consider the effect of increasing subposterior correlation, in Section 4.2 we consider the robustness with increasing numbers of sub-posteriors, and in Section 4.3 we consider how to address conflicting sub-posteriors. When applying Generalised Monte Carlo Fusion (or calling it as an embedded algorithm in Divide-and-Conquer Monte Carlo Fusion), in all cases we use the GPE-2 variant of Algorithm 1 discussed in Definition 2.2 of Section 2, and follow the direction of Fearnhead et al. [10] by estimating the mean of the Negative Binomial distribution in (16) by using the Trapezoidal rule and setting β c = 10. We note that different choices for the parameters of the Negative Binomial distribution in Definition 2.2 may affect the efficiency and variance of our algorithm, but does not introduce any bias to the estimator.
In order to compare Monte Carlo Fusion, Generalised Monte Carlo Fusion and Divide-and-Conquer Monte Carlo Fusion, we compute both the computational run-times of each methodology and a metric which we term the Integrated Absolute Distance (IAD). To compute the IAD we average across each dimension the difference between the true target (fusion) density (f ), and a kernel density estimate of the density using draws realised using our chosen methodology (f ). In particular, In the case where the true density is not available analytically, we take as a proxy for the true density a kernel density estimate using realisations from an MCMC run.

Effect of correlation
In this example, we consider the illustrative case in which we wish to recover a bi-variate Gaussian target distribution, f ∝ f 1 f 2 , where f c ∼ N 2 (0, Σ) with, Σ = 1.0 ρ corr ρ corr 1.0 .
As we are only considering combining two sub-posteriors in this section, we in effect consider only the Generalised Monte Carlo Fusion (GMCF) approach of Section 2. To study the impact of sub-posterior correlation on the robustness of Monte Carlo Fusion (MCF) and GMCF we can simply consider varying the single parameter ρ corr , and (in this case) compute the Effective Sample Size (ESS) per second averaged across 50 runs in order to compare the efficiency of each methodology. For simplicity, we assume we are able to sample directly from each subposterior, and for both methodologies we set T = 1. For the purposes of the GMCF approach of Algorithm 1, we simply set Λ c =Σ c , whereΣ c is the estimated covariance matrix from the sub-posterior samples for c = 1, 2 (and so in effect we have incorporated global information into our proposals), and use a particle set size of N = 10000. The results are presented in Figure  3, which clearly show that GMCF is robust to increasing sub-posterior correlation, and offers a significant computational advantage over MCF (which in this case exhibits a strong degradation in efficiency and performance).

Effect of hierarchy
In this example, we consider the illustrative case of attempting to recover a univariate standard Gaussian target distribution. In particular, we have f ∝ C c=1 f c , where f c ∼ N (0, C) for c = 1, . . . , C. By simply varying C, we can study the robustness with increasing numbers of sub-posteriors of Monte Carlo Fusion (MCF) (in effect the fork-and-join approach illustrated in Figure 1), and both our suggested versions of Divide-and-Conquer Monte Carlo Fusion (the balanced-binary tree approach illustrated in Figure 2a, and the progressive tree approach illustrated in Figure 2b). Note that in our chosen idealised setting, there is no advantage conferred with our embedded Generalised Monte Carlo Fusion methodology of Section 2, and so we are simply contrasting hierarchies. In all cases we use a particle set of size N = 10000 with resampling if ESS < N/2, set T = 1, use an appropriately scaled identity as the preconditioning (scalar) matrix, and average across 50 runs. The results are presented in Figure 4, which clearly show that, in contrast to the fork-and-join tree approach, both the balanced-binary tree and progressive tree approaches are robust in recovering the correct posterior distribution in the case of increasing C at the cost of modestly increased computational cost.

Dealing with conflicting sub-posteriors
Directly unifying C conflicting sub-posteriors (sub-posteriors which have little common support and have high total-variation distance) using a fork-and-join approach as in Monte Carlo Fusion [7] and Figure 1 is impractical. This can be understood with reference to (7) and (8), which indicates that importance weights will degrade rapidly in this setting.
An approach to deal with conflicting sub-posteriors is to temper the sub-posteriors (to an inverse temperature β ∈ (0, 1] such that there is sufficient sub-posterior overlap), and then propose a suitable tree for which the recursive Divide-and-Conquer Monte Carlo Fusion approach we introduced in Section 3 could then be applied to recover (1). In particular,   One such generic tree is provided in Figure 5, in which the tempered sub-posteriors are first unified into 1/β ∈ N tempered posteriors, which are then again unified into f . : : : · · · f β C ··· d d Figure 5: Illustrative tree approach for the Fusion problem in the case of conflicting sub-posteriors as in Section 4.3. 1/β copies of the C tempered (and over-lapping) subposteriors represent the leaves of the tree, which are unified into 1/β tempered versions of f (using a suitable tree and D&C-MCF as in Section 3), and then unified again (using another tree, and D&C-MCF) to recover f .
To illustrate the advantage of our D&C-MCF and tempering approach in the case of conflicting sub-posteriors, we consider the scenario of unifying two Gaussian sub-posteriors with the same variance (1), but with different mean (±µ). In particular, we have f ∝ f 1 f 2 where f 1 ∼ N (−µ, 1) and f 2 ∼ N (µ, 1). By simply increasing µ we can emulate increasingly conflicting sub-posteriors and study how the Monte Carlo Fusion approach of Dai et al. [7] (which is equivalent to the fork-and-join approach of Figure 1), behaves in terms of the IAD metric and computational time. We contrast this with our tempering approach, considering a range of temperatures 1/β ∈ {2, 4, 8, 16}, and then following the guidance of Figure 5. In particular, we use our D&C-MCF approach to unify the tempered sub-posteriors with the balanced-binary approach of Figure 2a for both the first and second stage in Figure 5. In all cases, we use a particle set size of N = 10000 with resampling if ESS < N/2, set T = 1, and average across 50 runs. The results are presented in Figure 6, and show clearly that our D&C-MCF approach is significantly more robust to conflicting sub-posteriors than the MCF approach where no tempering is applied. A natural trade-off arises when applying the tempering approach suggested, in that decreasing β results in tempered sub-posteriors which are less conflicting and are easier to combine, but there is an increased computational cost in recovering f as an increased number of levels are added to the resulting tree.

Examples
In this section, we consider a logistic regression model applied to a simulated data set (Section 5.1), and a real credit card data set (Section 5.2). For each data set, we compare the performance (in terms of computational run-time and the Integrated Absolute Distance) of the Divide-and-Conquer Monte Carlo Fusion (D&C-MCF) method we introduced in Section 3, against the other established approximate methodologies we discussed in the introduction. As in Section 4, we use the GPE-2 variant of Algorithm 1 discussed in Definition 2.2 of Section 2 with the same guidance on parameterisation.
For both examples, we consider the setting where the data is first split into C disjoint subsets and the simulation of each sub-posterior is conducted separately. For this reason, we focus on the balanced-binary tree approach (Figure 2a), as for a fixed data set this is the most natural hierarchy (whereas, the progressive tree approach in Figure 2b is more naturally suited to an online setting). The established methodologies we contrast our implementation (which can be found at https://github.com/rchan26/hierarchicalFusion) against are Consensus Monte Carlo (CMC) [28] (implemented using the parallelMCMCcombine package in R available from https://github.com/cran/parallelMCMCcombine [23]), the kernel density averaging approach of Neiswanger et al. [24] (which we term KDEMC, and also implemented using the parallelMCMCcombine R package), and the Weierstrass Rejection Sampler (WRS) [30] (implemented using the R code available at https://github.com/wwrechard/weierstrass). For each example, as a benchmark for f (in terms of a reference in computing IAD), we use Stan [5] on the entire data set, together with an appropriate choice of prior, to find a reference sample approximation of the desired f .

Simulated data example
We begin by considering a variation of the example in Scott et al. [28,Section 4.3], we construct a small data set of size m = 1000, to which we applied a logistic regression model: Each record of the simulated design matrix contained four covariates in addition to an intercept. The ith entry of the design matrix is given by are random variables generated from a mixture density with a point-mass at zero (and so are either activated or not). In particular, we have for j ∈ {1, . . . , 4} that ζ i,j ∼ p j N (1, 1)+(1−p j )δ 0 . To study our methodology, we begin by splitting the data set of size m = 1000 equally between a number of cores. Here we consider C ∈ {4, 8, 16} subsets. On each data set on each core, we fit the logistic regression model using MCMC with Stan using a Gaussian prior distribution with mean 0 and variance 1/C on each of the parameters, resulting in the required C (sample approximation) sub-posteriors. The resulting sub-posteriors in the case of C = 16 naturally resulted in conflicting sub-posteriors (see Section 4.3) due to the small amount of data on each core. We then applied our D&C-MCF approach using a balanced-binary tree with N = 10000 and T = 0.5, and compared our approach to CMC, KDEMC and WRS. The results are shown in Figure 7, with reference to the Stan benchmark.
It is clear from Figure 7a that D&C-MCF achieves the best sample approximation in terms of IAD of any of the approaches considered. Furthermore, the quality of the sample approximation is robust to the increasing sub-posteriors that we consider. In terms of computational cost, CMC outperforms all other methodologies, but has by far the poorest performance in terms of IAD. Although more expensive than the WRS, D&C-MCF has a comparable computational cost to KDEMC. Further D&C-MCF appears to achieve a quality sample approximation at a cost which grows at roughly the same rate as the other approaches.

Credit-card data example
In this example, we consider the 'Default of credit card clients' data set available from the UCI Machine Learning Repository at https://archive.ics.uci.edu/ml/datasets [31], and again fit the logistic regression model of (21). The data set comprised m = 30000 records, each of which contained response variable of whether a default had occurred (which we treated as y i = 1), or not (in which case y i = 0). From the data set we used the X2: Gender attribute, treating it as a binary covariate with 1 being male and 0 female. In addition we used the X3: Education attribute to create three further binary covariates: a binary corresponding to whether the individual had completed high school; a binary corresponding to completion of university; and a further binary corresponding to completion of graduate school. These three education covariates were chosen to induce strong correlation in the resulting inference.
We again split the data set of size m = 30000 equally between C cores, and again used Stan  together with a Gaussian prior distributions with mean 0 and variance 1/C on each of the parameters, to arrive at our C sub-posteriors. In this example, we consider C ∈ {32, 64, 128} subsets, which we then attempt to unify. This example is particularly challenging as the data and sub-posteriors exhibit large correlation (due to the covariates related to X3: Education), especially for large C. D&C-MCF using a balanced-binary tree with N = 30000 and T = 0.5, together with CMC, KDEMC and WRS were then used to unify the C sub-posteriors. The results are shown in Figure 8 with reference to the Stan benchmark.
The results in Figure 8 are comparable to those of Section 5.1: D&C-MCF achieves the best IAD of all methodologies, and is robust to increasing C. This comes at a moderate fixed computational cost, which scales no worse than any other methodology (including CMC), and indeed in this example has a computational cost which is lower than KDEMC.

Conclusion
The Fusion approach to unifying sub-posteriors into a coherent sample approximation of the posterior (as in (1)), offers fundamental advantages over approximation based approaches. In particular, Fusion avoids having to impose any distributional approximation on the sub-posteriors, and so is more robust to a wider range of models, and circumvents needing to understand the impact of imposed approximations on the unified posterior. To date, Fusion approaches (such as the Monte Carlo Fusion approach of Dai et al. [7]) have had impractical computational cost in realistic settings, lacking robustness when considering: the number of sub-posteriors being unified; when unifying highly correlated sub-posteriors; and when considering conflicting subposteriors. In this paper, we have substantially addressed the practical issues of existing Fusion approaches by means of a number of theoretical and methodological extensions of the original framework. In particular, in Section 2 we introduced Generalised Monte Carlo Fusion (GMCF), which is an importance sampling approach to Fusion that incorporates available global information for each sub-posterior in order to construct informative proposals. As shown in Section 4.1, GMCF addresses the lack of robustness when the sub-posteriors have strong correlation structure. By embedding GMCF within a Divide-and-Conquer Sequential Monte Carlo (D&C-SMC) framework [21,19], in Section 3 we introduced Divide-and-Conquer Monte Carlo Fusion (D&C-MCF), together with a number of tree hierarchies, which allow the sub-posteriors to be combined in stages to recover f . As demonstrated in Section 4.2, D&C-MCF is a robust approach to unifying large numbers of sub-posteriors. Furthermore, as shown in Section 4.3, even in the setting of conflicting sub-posteriors, D&C-MCF together with tempering and an appropriate hierarchy can result in a practical Fusion algorithm. Finally, in Section 5 we applied our D&C-MCF approach to realistic data sets and compared its performance with competing approximate methodologies. In all of these settings, our implementation of D&C-MCF offered the best performance in terms of Integrated Absolute Distance to an appropriate benchmark, at a modest computational cost.
There are a number of interesting avenues for extending the work of this paper. Perhaps most interesting is to adapt the D&C-MCF approach to constraints in practical settings. As discussed in the introduction, one particularly promising direction is when considering (1) under privacy constraints of the individual sources [32]. In this setting, we may have a number of parties that wish to combine their distributional analysis on a common parameter space and model but cannot reveal their distribution due to confidentiality. This is an active area of research of the authors, and motivates variant tree hierarchies in D&C-MCF. Another application instance is when considering a truly distributed 'big data' setting where communication between different cores is expensive [28]. Unlike Monte Carlo Fusion [7], where communication between cores can be reduced to a single interaction per independent realisation from (1), in the case of D&C-MCF, methodology would need to be developed to mitigate the additional communication from a naive implementation. One approach would be to investigate embedding a sub-sampling approach [27].
From a theoretical perspective, current Fusion methodologies only consider sub-posteriors on a common parameter space. One direction of interest is extending Fusion methodology to combine sub-posteriors with varying dimension. The Markov Melding framework of Goudie et al. [15] where separate sub-models (potentially of differing dimension) are fitted to different data sources and then joined, is promising. In this setting, the tree hierarchies could be defined by the model itself. To mitigate computational robustness of Fusion with increasing dimension in this setting, it may be possible to further utilise the methodology in Lindsten et al. [21].
From a methodological perspective, there are two clear directions. One is to embed within D&C-MCF the Bayesian Fusion (BF) approach of Dai et al. [8]. BF is a sequential Monte Carlo (SMC) approach which reduces the cost of Fusion approaches to linear in T , and so may further improve the computational scaling of D&C-MCF. The second is to pursue the tempering within D&C-MCF approach of Section 4.3, and explore its practical usage as an alternative to parallel tempering approaches within truly multi-modal problems.