# The divide-and-conquer sequential Monte Carlo algorithm: theoretical
properties and limit theorems††thanks: JK and AMJ acknowledge support from the
EPSRC (grant # EP/T004134/1) and the Lloyd’s Register Foundation Programme on
Data-Centric Engineering at the Alan Turing Institute. FRC acknowledges
support from the EPSRC and the MRC OXWASP Centre for Doctoral Training (grant
# EP/L016710/1). FRC and AMJ acknowledge further support from the EPSRC (grant
# EP/R034710/1).

Juan Kuntz Department of Statistics, University of Warwick, Coventry, CV4 7AL,
UK. Email: $\\{$juan.kuntz-nussio, francesca.crucinio,
<EMAIL_ADDRESS>Francesca R. Crucinio22footnotemark: 2 Adam M.
Johansen22footnotemark: 2

###### Abstract

We revisit the divide-and-conquer sequential Monte Carlo (DaC-SMC) algorithm
and firmly establish it as a well-founded method by showing that it possesses
the same basic properties as conventional sequential Monte Carlo (SMC)
algorithms do. In particular, we derive pertinent laws of large numbers,
$L^{p}$ inequalities, and central limit theorems; and we characterize the bias
in the normalized estimates produced by the algorithm and argue the absence
thereof in the unnormalized ones. We further consider its practical
implementation and several interesting variants; obtain expressions for its
globally and locally optimal intermediate targets, auxiliary measures, and
proposal kernels; and show that, in comparable conditions, DaC-SMC proves more
statistically efficient than its direct SMC analogue. We close the paper with
a discussion of our results, open questions, and future research directions.

##  1 Introduction

Distributed inference algorithms that scale efficiently to large models and
data sets are of widespread interest at the present time. Divide-and-conquer
methods, which recursively decompose an inference problem into smaller
problems whose solutions can be combined in a principled manner (e.g. see
[43]), provide a promising approach in this setting. In computational
statistics many methods have been developed which allow for distributed
inference either by performing inference for a subset of the data at each
processor and ultimately combining their results [16, 45, 47, 54, 58] or by
carrying out part of the computation associated with each step of an algorithm
at a separate processor and minimizing the necessary communication [29, 50,
52, 56]. Another strategy which is beginning to find applications ranging from
general distributed Bayesian inference [8] through inference in broad model
classes such as general state space hidden Markov models [13, 20, 21] to
inference in specific domains such as phylogenetics [30] is the _divide-and-
conquer sequential Monte Carlo_ (_DaC-SMC_) algorithm introduced in [41]. The
DaC-SMC algorithm recursively breaks down large inferential problems into
smaller, more manageable ones and merges their solutions in a systematic
manner. It is particularly amenable to distributed implementation.

The goal of this paper is to provide a comprehensive theoretical foundation
for this algorithm. To achieve this, we view the DaC-SMC algorithm here as a
merger of product-form estimators [37] and sequential Monte Carlo (SMC), two
classes of Monte Carlo estimators that are (separately) well-understood.

Product-form estimators encompass a class of generalized U-statistics [39, 36]
that exploit conditional independences present in the measure targeted, or the
proposal distribution employed, to achieve lower variances than standard Monte
Carlo estimators [37]. The idea underpinning them is very simple: if
$(X^{1}_{1},\dots,X^{N}_{1}),\dots,(X^{1}_{K},\dots,X^{N}_{K})$ are $K$
independent sequences of samples respectively drawn from probability
distributions $\mu_{1},\dots,\mu_{K}$, then every ‘permutation’ of these
samples has the same product law:

$(X^{n_{1}}_{1},\dots,X^{n_{K}}_{K})\sim\mu:=\mu_{1}\times\cdots\times\mu_{K}\quad\forall(n_{1},\dots,n_{K})\in\\{1,\ldots,N\\}^{K}.$

In other words, using only $N$ samples from each marginal
$\mu_{1},\dots,\mu_{K}$, we can construct $N^{K}$ samples from their product
$\mu$. Product-form estimators approximate $\mu$, and averages thereof, using
the empirical distribution of these $N^{K}$ samples (as opposed to the default
choice involving the empirical distribution of the $N$ ‘unpermuted’ samples
$(X^{n}_{1},\dots,X^{n}_{K})_{n=1}^{N}$). This increase in sample number leads
to a decrease in estimator variance, and we say that product-form estimators
are more _statistically efficient_ than their conventional counterparts. The
empirical distribution $\mu_{\times}^{N}$ of the $N^{K}$ samples can be
expressed as the product of the ‘marginal’ empirical distributions
$\mu_{1}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{X^{n}_{1}},\dots,\mu_{K}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{X^{n}_{K}}$
(hence, the term ‘product-form estimators’). For this reason, product-form
estimators embody the fork-and-join model [12]: they ‘fork’ the problem of
approximating $\mu$ into that of approximating its marginals, solve these
smaller problems by computing $\mu_{1}^{N},\dots,\mu_{K}^{N}$ separately, and
then ‘join’ their solutions by taking their product, so obtaining a solution
($\mu_{\times}^{N}$) for the original problem.

Product-form estimators find their greatest use not in isolation, but when
embedded within more sophisticated Monte Carlo methodology to tackle the
aspects of the problem exhibiting product structure or conditional
independences (e.g. [55, 1, 53, 37]). Here, we study such embeddings within
SMC algorithms, a class of Monte Carlo methods that approximate sequences of
distributions defined on state spaces of increasing dimension by iteratively
mutating and selecting weighted populations of samples or _particles_ (see,
for example, [11] for a recent book-length treatment). These algorithms are
very widely employed in the context of hidden Markov models [23, 33, 25] and,
following [9, 19] and related works, have been shown to be a useful
alternative to Markov chain Monte Carlo for approximating general
distributions.

Incorporating product-form estimators into sequential Monte Carlo methodology
one obtains the DaC-SMC algorithm introduced in [41]. Its purpose is to
estimate targets on medium-to-high dimensional spaces that are products of
collections of lower dimensional spaces. In order to accomplish this, an
auxiliary tree is defined. Each node of this tree has some subset of the
variables within the original model assigned to it and each subtree of the
tree has associated with it an intermediate model over the variables assigned
to the nodes in the subtree. The algorithm progressively evolves clouds of
weighted particles up the tree (from leaves to root) so that the corresponding
weighted empirical distributions target the intermediate distributions. At
each node, the clouds corresponding to the node’s children are merged using a
product-form estimator. Hence, the algorithm exemplifies the divide-and-
conquer paradigm: carrying computations in the partial order dictated by the
indexing tree, DaC-SMC recursively constructs product-form estimators for each
of the lower-dimensional intermediate targets and, ultimately, for the
original target itself.

While some preliminary work was done in [41], the algorithm’s theoretical
properties remain underexplored. The main aim of this paper is to fill this
gap and firmly establish DaC-SMC as a well-founded algorithm: we show that it
possesses the same fundamental theoretical properties as conventional SMC
(Section 1.1). We also have several other, more minor, associated aims:
(Section 3.3) exploring efficient implementations of the algorithm; (Section
3.4) identifying several interesting variants thereof in addition to those
already pointed out in [41, Section 4]; (Section 4.2) deriving the optimal and
‘locally optimal’ intermediate targets, auxiliary measures, and proposal
kernels and giving guidance on their choice in practice; and (Section 4.3)
formally arguing that DaC-SMC estimators are more statistically efficient than
their direct SMC analogues whenever both algorithms are placed on an equal
footing (this does not, however, necessarily mean that DaC-SMC is more
computationally efficient than SMC as the former’s cost is generally greater
than the latter’s; see Sections 3.3 and 4.3 for more on this).

###  1.1 Relation to the literature and contributions

The most relevant publication is [41] where the DaC-SMC algorithm was
introduced, various variants were considered, several numerical examples were
given, and it was shown that some forms of DaC-SMC estimators for the
normalized targets are consistent and those for the normalizing constants are
unbiased (cf. Section 3.1 for definitions the ‘unnormalized targets’,
‘normalized targets’, and ‘normalizing constants’). Both in the original paper
[41] and also in later work [21], the algorithm was found to have good
empirical properties. Applications of the product-form estimators it builds on
can be found peppered throughout the Monte Carlo literature [55, 41, 1, 53,
37], almost always unnamed and specialized to particular contexts. Their
theoretical properties, at least in any generality, were seemingly only
studied recently in [37] which established their consistency, unbiasedness,
asymptotically normality, and statistical efficiency, examined their
computational efficiency, and considered several variants and extensions.

The contributions of this paper are several:

  * •

_DaC-SMC’s theoretical characterization._ We view the DaC-SMC algorithm
through the lens of product-form estimators (Sections 3.2, 3.3). Using the
insights afforded by this perspective, we obtain a comprehensive description
of the algorithm’s theoretical properties (Section 4). In particular, we show
that DaC-SMC estimators for both the unnormalized and normalized targets
retain the key properties of their SMC analogues: they are consistent
(Theorems 1, 2), asymptotically normal (Theorem 6), and satisfy $L^{p}$
inequalities (Theorem 5). In the case of the unnormalized targets, we also
show that the estimators are unbiased (Theorem 3), while in that of the
normalized ones, we show that the estimators’ bias decays linearly with the
number of particles (Theorem 4). Our proofs here combine well-known arguments
used to analyze SMC algorithms with novel techniques required to deal with the
product-form aspects of the estimators. To be more specific about the former,
for the laws of large numbers and $L^{p}$ inequalities we adapt results from
[14, 44, 17], for the asymptotic normality we borrow from [17, 18], for the
unbiasedness of the unnormalized target estimators we emulate [2], and for the
bias bound of the normalized target estimator we follow [48]; see Section 4.1
for further details.

  * •

_Optimal intermediate targets, auxiliary measures, and proposal kernels._
Using the asymptotic variance expressions obtained in Theorem 6, we derive
(Section 4.2) the intermediate targets, auxiliary measures, and proposal
kernels (cf. Sections 3.1–3.2 for the definitions of these objects) that lead
to zero _asymptotic_ variance DaC-SMC estimators for the normalizing
constants. A bit unexpectedly, and in contrast with the standard case, we find
that these choices do not necessarily lead to these estimators possessing zero
_finite-sample_ variance. We also generalize the well-known characterization
of the ‘locally optimal’ auxiliary measures and proposal kernels of standard
SMC to the DaC-SMC setting (Theorem 8).

  * •

_Statistical efficiency of DaC-SMC._ Also using the asymptotic variances
expressions in Theorem 6, we show that DaC-SMC estimators are, at least
asymptotically, more statistically efficient than their SMC counterparts
(Section 4.3): the asymptotic variances of the former are bounded above by
those of the latter.

###  1.2 Paper structure

Section 2 reviews the two building blocks underpinning the DaC-SMC algorithm:
product-form estimators (Section 2.1) and sequential Monte Carlo (Section
2.2). Section 3 presents the DaC-SMC algorithm: Section 3.1 describes the
target measures the algorithm approximates, Section 3.2 the algorithm itself,
Section 3.3 efficient implementations thereof, and Section 3.4 several
interesting variations (Appendix A contains an example illustrating concepts
discussed in Sections 3.3, 3.4). Section 4 details the algorithm’s theoretical
properties: Section 4.1 describes the methods of proof we use to argue the
properties (to ease the reading, we defer our proofs to Appendices B–I),
Section 4.2 studies the choice of intermediate targets, auxiliary measures,
and proposal kernels, and Section 4.3 compares DaC-SMC with standard SMC
theoretically and computationally. We close the paper with a discussion in
Section 5.

##  2 Background: Product-form estimators and sequential Monte Carlo

The DaC-SMC algorithm combines ideas from two distinct Monte Carlo methods:
product-form estimators and sequential Monte Carlo algorithms. Here, we
provide a brief overview of each.

###  2.1 Product-form estimators

Product-form estimators [37] are Monte Carlo methods that exploit any
conditional independences present in the measure they target, or the proposal
distribution they use, to achieve higher statistical efficiency than normal.
In particular, consider the task of estimating

$\mu(\varphi):=\int\varphi(x)\mu(dx)$

for a given probability distribution $\mu$ on a measurable space
$(E,\mathcal{E})$ and a regular enough function $\varphi$. Suppose that $\mu$
is product-form: it is the product $K$ distributions $\mu_{1},\dots,\mu_{K}$
on respective measurable spaces
$(E_{1},\mathcal{E}_{1}),\ldots,(E_{K},\mathcal{E}_{K})$ satisfying
$(E,\mathcal{E})=\left(\prod_{k=1}^{K}E_{k},\prod_{k=1}^{K}\mathcal{E}_{k}\right)$,
where $\prod_{k=1}^{K}\mathcal{E}_{k}$ denotes the corresponding product
$\sigma$-algebra. To accomplish this task, suppose that we are given $N$
samples $(\bm{X}^{n})_{n=1}^{N}=(X^{n}_{1},\dots,X^{n}_{K})_{n=1}^{N}$
independently drawn from $\mu$. Instead of approximating $\mu(\varphi)$ with
the standard Monte Carlo estimator,

$\mu^{N}(\varphi):=\frac{1}{N}\sum_{n=1}^{N}\varphi(X^{n}),$

consider using the _product-form estimator_ ,

$\mu_{\times}^{N}(\varphi):=\frac{1}{N^{K}}\sum_{\bm{n}\in[N]^{K}}\varphi(\bm{X}^{\bm{n}}),$

obtained by averaging over all of the $N^{K}$ element-wise permutations of
samples $(\bm{X}^{n})_{n=1}^{N}$,

$\bm{X}^{\bm{n}}:=(X^{n_{1}}_{1},\dots,X^{n_{K}}_{K})\quad\forall\bm{n}=(n_{1},\dots,n_{K})\in[N]^{K},\quad\text{where}\quad[N]:=\\{1,\dots,N\\}.$

It is a strongly consistent, unbiased, and asymptotically normal estimator for
$\mu(\varphi)$ ([37, Theorem 1]) that achieves a smaller variance than
$\mu^{N}(\varphi)$ ([37, Corollary 1]),

$\text{Var}(\mu^{N}_{\times}(\varphi))\leq\text{Var}(\mu^{N}(\varphi))\quad\forall
N>0,$

for all square $\mu$-integrable test functions $\varphi$. (For that matter,
$\mu_{\times}^{N}(\varphi)$ yields the best unbiased estimates of
$\mu(\varphi)$ that can be achieved using only the knowledge that $\mu$ is
product-form and $N$ i.i.d. samples drawn from $\mu$, cf. [37, Theorem 2].) In
a nutshell, the standard estimator uses the (non-product-form) empirical
distribution of the unpermuted samples,

$\mu^{N}:=\frac{1}{N}\sum_{n=1}^{N}\delta_{\bm{X}^{n}}$

as an approximation to $\mu$, while the product-form estimator uses the
(product-form) empirical distribution obtained by taking the product of the
marginal empirical distributions,

$\mu^{N}_{\times}:=\left(\frac{1}{N}\sum_{n_{1}=1}^{N}\delta_{X^{n_{1}}_{1}}\right)\times\dots\times\left(\frac{1}{N}\sum_{n_{K}=1}^{N}\delta_{X^{n_{K}}_{K}}\right)=\frac{1}{N^{K}}\sum_{\bm{n}\in[N]^{K}}\delta_{\bm{X}^{\bm{n}}};$

and the latter provides a better approximation to $\mu$ for the reasons
outlined in the introduction (see [37, Section 2.2] for further details).

The above applies similarly to the case where $\mu$ is an importance sampling
(IS) proposal rather than a target. In particular, if our target $\gamma$ is a
change of measure from $\mu$ with Radon-Nikodym derivative $w:=d\gamma/d\mu$,
then

$\gamma_{\times}^{N}(\varphi):=\mu_{\times}^{N}(w\varphi)=\frac{1}{N^{K}}\sum_{\bm{n}\in[N]^{K}}w(\bm{X}^{\bm{n}})\varphi(\bm{X}^{\bm{n}}),$

is a consistent, unbiased, and asymptotically normal estimator for
$\gamma(\varphi)=\mu(w\varphi)$ whose variance is bounded above by that of the
standard IS estimator $\gamma^{N}(\varphi):=\mu^{N}(w\varphi)$; cf. [37,
Section 3.1]. Because this is the case for all sufficiently integrable
$\varphi$, we regard

$\gamma^{N}_{\times}(d\bm{x}):=w(\bm{x})\mu_{\times}^{N}(d\bm{x})=\frac{1}{N^{K}}\sum_{\bm{n}\in[N]^{K}}w(\bm{X}^{\bm{n}})\delta_{\bm{X}^{\bm{n}}}(d\bm{x})$
(1)

as an approximation to $\gamma$. Somewhat abusively, we also refer to this
type of approximation as a ‘product-form estimator’ (e.g.
$\gamma^{N}_{\times}$ is a product-form estimator for $\gamma$).

###  2.2 Sequential Monte Carlo

Sequential Monte Carlo (SMC) algorithms constitute a class of Monte Carlo
methods that, by progressively evolving clouds of weighted particles,
approximate sequences $\rho_{0},\dots,\rho_{T}$ of intractable unnormalized
measures defined on measurable spaces
$(\bm{E}_{0},\bm{\mathcal{E}}_{0}),\dots,(\bm{E}_{T},\mathcal{E}_{T})$ of
increasing dimension (in particular, $\bm{E}_{t}=E_{0}\times\cdots\times
E_{t}$ for some low-dimensional $E_{0},\dots,E_{T}$ and
$\bm{\mathcal{E}_{t}}=\mathcal{E}_{0}\times\cdots\times\mathcal{E}_{t}$ is the
appropriate product $\sigma$-algebra). They also yield weighted particle
approximations of the corresponding normalized distributions
$\mu_{0},\dots,\mu_{T}$ and normalization constants $Z_{0},\dots,Z_{T}$.

The basic SMC approach, often referred to as sequential importance resampling
[23, p. 15], is detailed in Algorithm 1 (we use slightly unconventional
notation here to keep things consistent with Section 3). It is a simple
iterative procedure: given a weighted-particle approximation $\rho_{t-1}^{N}$
to $\rho_{t-1}$, normalize $\rho_{t-1}^{N}$ to obtain an approximation
$\mu^{N}_{t-1}$ of $\mu_{t-1}$, _resample_ $\mu^{N}_{t-1}$ to obtain
unweighted particles $\bm{X}_{t_{-}}^{1,N},\ldots,\bm{X}_{t_{-}}^{N,N}$ (whose
empirical measure also approximates $\mu_{t-1}$), and extend the path of each
resampled particle $\bm{X}_{t_{-}}^{n,N}$ by drawing $X_{t}^{n,N}$ from a
proposal kernel $K_{t}(\bm{X}_{t_{-}}^{n,N},dx_{t})$ and setting
$\bm{X}_{t}^{n,N}:=(X^{n,N}_{t},\bm{X}_{t_{-}}^{n,N})$. Then, re-weight the
particles using the Radon-Nykodim derivative
$w_{t}:=d\rho_{t}/d(\rho_{t-1}\otimes K_{t})$ of $\rho_{t}(d\bm{x}_{t})$
w.r.t. $(\rho_{t-1}\otimes
K_{t})(d\bm{x}_{t}):=\rho_{t-1}(d\bm{x}_{t-1})K_{t}(\bm{x}_{t-1},dx_{t})$ to
obtain an approximation of $\rho_{t}$:

$\displaystyle\rho_{t}^{N}:=\frac{{Z}_{t-1}^{N}}{N}\sum_{n=1}^{N}w_{t}(\bm{X}_{t}^{n,N})\delta_{\bm{X}_{t}^{n,N}},$
(2)

where $Z_{t-1}^{N}:=\rho_{t-1}^{N}(\bm{E}_{t-1})$ denotes the estimate of the
normalizing constant $Z_{t-1}=\rho_{t-1}(\bm{E}_{t-1})$ afforded by
$\rho^{N}_{t-1}$. To initialize the algorithm, we approximate $\rho_{0}$ using
importance sampling with a proposal distribution $K_{0}$.

Resampling is the process of stochastically replicating some particles within
a weighted collection in such a way that the expected number of offspring of
each particle is proportional to its weight. By doing so, we obtain an
unweighted sample suitable for approximating the same target distribution as
the original weighted sample. The simplest way to achieve this, known as
multinomial resampling because the vector of offspring numbers follows a
multinomial distribution, is to independently draw $N$ times from the weighed
empirical measure of the original sample. To not complicate this paper’s
exposition and proofs, and to keep the focus on the differences between
standard SMC and DaC-SMC, we restrict ourselves throughout to the case of
multinomial resampling. However, using more sophisticated resampling schemes
within the DaC-SMC setting is sensible, incurs no further complications than
for SMC, and can in principle be understood using similar theoretical
techniques to those employed in the conventional setting (see, in particular,
[27]).

1: Propose: for $n\leq N$, draw $\bm{X}^{n,N}_{0}$ independently from $K_{0}$.

2: Correct: compute
$\rho^{N}_{0}:=N^{-1}\sum_{n=1}^{N}w_{0}(\bm{X}^{n}_{0})\delta_{\bm{X}^{n}_{0}}$,
where $w_{0}:=d\rho_{0}/dK_{0}$, and
$\mu^{N}_{0}:=\rho^{N}_{0}/\rho^{N}_{0}(\bm{E}_{0})$.

3: for $t=1,\dots,T$ do

4: Resample: for $n\leq N$, draw $\bm{X}_{t_{-}}^{n,N}$ independently from
$\mu_{t-1}^{N}$.

5: Mutate: for $n\leq N$, draw $X^{n,N}_{t}$ independently from
$K_{t}(\bm{X}_{t_{-}}^{n,N},dx_{t})$ and set
$\bm{X}_{t}^{n,N}:=(X^{n,N}_{t},\bm{X}_{t_{-}}^{n,N})$.

6: Correct: compute $\rho_{t}^{N}$ in (2) and
$\mu^{N}_{t}:=\rho_{t}^{N}/\rho_{t}^{N}(\bm{E}_{t})$.

7: end for

Algorithm 1 Sequential importance resampling

SMC algorithms were originally introduced in the filtering literature (see
[23, 25] for reviews in this context) and have been successfully employed to
sample from general distributions [9, 19]. The filtering problem has an
intrinsic sequential structure and the sequence $\rho_{0},\dots,\rho_{T}$
arises naturally. For general distributions $\mu$ lacking this sort of
structure, SMC can instead be applied by building an ‘artificial’ sequence of
distributions such that $\mu_{T}$ (or one of its marginals) coincides with
$\mu$. This can be achieved by incrementally morphing an easy-to-sample-from
distribution into $\mu$, often by tempering [46] or sequentially introducing
observations to transition from a prior to a posterior [9]. At first glance,
Algorithm 1 may not appear applicable in such settings because the artificial
distributions are all defined on the same space. However, it is possible to
construct auxiliary sequences of distributions on spaces of increasing
dimension such that the $t^{\textrm{th}}$ distribution of interest coincides
with the $t^{\textrm{th}}$ ‘coordinate marginal’ of the corresponding
auxiliary distribution, hence allowing _exactly_ this algorithm to be applied
to these problems [19]. The same kind of constructions can be used in the DaC-
SMC framework (see [41, Section 4.2] for a brief discussion) and we focus on
analysing the fundamental algorithm (which, combined with appropriately
extended target distributions, applies directly to these seemingly more
complex settings).

We mention one variant of the SMC algorithm which is of particular relevance
to what follows: the auxiliary SMC (ASMC) scheme originally proposed for
filtering [51]. Although it was initially presented as an alternative proposal
method involving an auxiliary variable construction, this scheme can be
understood as the standard algorithm with an additional re-weighting prior to
the resampling step [7]—the motivation behind the additional re-weighting
being to incorporate as much information about the next distribution in the
sequence as possible _prior_ to resampling. Simply put, it is an SMC algorithm
that targets a slightly different sequence of distributions and employs an
extra importance sampling correction to obtain estimates with respect to the
distributions of interest [31].

In particular, we introduce _auxiliary measures_
$\gamma_{1_{-}},\ldots,\gamma_{T_{-}}$ on $\bm{E}_{0},\dots,\bm{E}_{T-1}$, and
their normalized counterparts $\pi_{1_{-}},\ldots,\pi_{T_{-}}$, and run
Algorithm 1 with these measures as the targets. That is, with the weights used
for resampling at time $t$ set to $w_{t-}:=d\gamma_{t_{-}}/d\gamma_{t-1}$,
where $\gamma_{0}:=K_{0}$ and $\gamma_{t}:=\gamma_{t_{-}}\otimes K_{t}$ for
$t=1,\dots,T$. This gives rise to the approximation

$\displaystyle\gamma_{t_{-}}^{N}:=\frac{\mathcal{Z}_{t-1}^{N}}{N}\sum_{n=1}^{N}w_{t_{-}}(\bm{X}_{t-1}^{n,N})\delta_{\bm{X}_{t-1}^{n,N}}$
(3)

of the auxiliary measure $\gamma_{t_{-}}$, where $\mathcal{Z}_{0}^{N}:=1$ and
$\mathcal{Z}_{t}^{N}:=\gamma_{t_{-}}^{N}(\bm{E}_{t-1})$ for $t=1,\dots,T$. In
the ASMC literature, it is common to find the correction step at the beginning
of iteration $t$ instead of the end of iteration $t-1$ (cf. Algorithm 2). Such
cyclical permutations of the resample, mutate, and correct steps have no
impact beyond labelling; we adopt this ordering because it eases the
presentation of the DaC-SMC algorithm in Section 3.

1: Propose: for $n\leq N$, draw $\bm{X}^{n,N}_{0}$ independently from $K_{0}$.

2: Compute: $\gamma_{0}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}^{n,N}_{0}}$.

3: for $t=1,\dots,T$ do

4: Correct: compute
$\gamma_{t_{-}}^{N}(d\bm{x}_{t-1}):=w_{t_{-}}(\bm{x}_{t-1})\gamma_{t-1}^{N}(d\bm{x}_{t-1})$
and $\pi_{t_{-}}^{N}:=\gamma_{t_{-}}^{N}/\gamma_{t_{-}}^{N}(\bm{E}_{t-1})$.

5: Resample: for $n\leq N$, draw $\bm{X}_{t_{-}}^{n,N}$ independently from
$\pi_{t_{-}}^{N}$.

6: Mutate: for $n\leq N$, draw $X^{n,N}_{t}$ independently from
$K_{t}(\bm{X}_{t_{-}}^{n,N},dx_{t})$ and set
$\bm{X}_{t}^{n,N}:=(X^{n,N}_{t},\bm{X}_{t_{-}}^{n,N})$.

7: Compute:
$\gamma_{t}^{N}:=N^{-1}\mathcal{Z}_{t}^{N}\sum_{n=1}^{N}\delta_{\bm{X}_{t}^{n,N}}$
where $\mathcal{Z}_{t}^{N}:=\gamma_{t_{-}}^{N}(\bm{E}_{t-1})$.

8: end for
Note: At each step $t$, one can use (4) to compute estimates of $\rho_{t}$,
$Z_{t}$, and $\mu_{t}$.

Algorithm 2 Auxiliary sequential importance resampling

The _extended auxiliary measure_ $\gamma_{t}$ is then approximated by

$\gamma_{t}^{N}:=\frac{\mathcal{Z}_{t}^{N}}{N}\sum_{n=1}^{N}\delta_{\bm{X}_{t}^{n,N}},$

where, in the case of multinomial resampling,
$\bm{X}_{t}^{1,N},\dots,\bm{X}_{t}^{N,N}$ denote samples drawn independently
from the normalization of $\gamma_{t_{-}}^{N}\otimes K_{t}$. (We note in
passing that $\gamma_{t-1}$ and $\gamma_{t_{-}}$ can be interpreted as the
‘predictive’ and ‘updated’ time marginals of a Feynman-Kac flow in the sense
of [17].) The above is then corrected using importance sampling and the weight
function $w_{t}:=d\rho_{t}/d\gamma_{t}$ to obtain approximations of the
objects of interest $\rho_{t}$, $Z_{t}$, and $\mu_{t}$:

$\displaystyle\rho_{t}^{N}(d\bm{x}_{t}):=w_{t}(\bm{x}_{t})\gamma_{t}^{N}(d\bm{x}_{t}),\quad{Z}_{t}^{N}:=\rho_{t}^{N}(\bm{E}_{t}),\quad\mu_{t}^{N}(d\bm{x}_{t}):=\frac{1}{{Z}_{t}^{N}}\rho_{t}^{N}(d\bm{x}_{t}),$
(4)

with $\gamma_{0}^{N}$ denoting the empirical distribution of the samples drawn
from $\gamma_{0}=K_{0}$. This construction gives rise to the auxiliary SMC
algorithm described in Algorithm 2.

##  3 The divide-and-conquer sequential Monte Carlo algorithm

In this section, we define our problem setting (Section 3.1), give the DaC-SMC
algorithm (Section 3.2), consider its implementation (Section 3.3), and
discuss various variants thereof (Section 3.4).

###  3.1 Problem definition

We tackle a generalization of the problem in Section 2.2: obtaining weighted-
particle approximations to a collection $(\rho_{u})_{u\in\mathbb{T}}$ of
intractable unnormalized (but finite) measures indexed by the nodes of a
finite rooted tree $\mathbb{T}$. In doing so, we also obtain tractable
approximations to the corresponding collections of normalizing constants
$(Z_{u})_{u\in\mathbb{T}}$ and normalized measures
$(\mu_{u})_{u\in\mathbb{T}}$. Typically, our main interest will lie in
approximating the target $\rho_{\mathfrak{r}}$ indexed by the tree’s root
$\mathfrak{r}$ and all other targets have been introduced to facilitate this
process. Hence, we refer to $\rho_{\mathfrak{r}}$ and $\mu_{\mathfrak{r}}$ as
the ‘final targets’ and to all other $\rho_{u}$ and $\mu_{u}$ as the
‘intermediate targets’.

Although this may seem a somewhat specialized setting at first glance, a great
many problems (including, for example, essentially any in which one wishes to
approximate a probability distribution over $\mathbb{R}^{d}$) may be viewed
within this framework. In some cases, the tree structure is intrinsic to the
model at hand and there are natural candidates for the intermediate targets
(e.g. [30, 49, 26] or Example 1 below). In others, the tree decomposition and
intermediate targets are entirely artificial constructs introduced for the
purpose of providing a computationally tractable approach to approximating the
final targets (e.g. see [41, Section 3.4] or [21]) and DaC-SMC serves as a
generalization of the approaches in [19, 9].

In all of these cases, the intermediate targets transition incrementally from
a collection of measures on low-dimensional spaces indexed by $\mathbb{T}$’s
leaves to the final target. Hence, the measures are defined on spaces whose
dimension grows as we progress up the tree. In particular, the underlying
measurable space $(\bm{E}_{u},\bm{\mathcal{E}}_{u})$ for $\rho_{u}$ is defined
by taking the corresponding partial product of a collection
$(E_{v},\mathcal{E}_{v})_{v\in\mathbb{T}}$ of low-dimensional spaces indexed
by $\mathbb{T}$’s nodes:

$\bm{E}_{u}=\prod_{v\in\mathbb{T}_{u}}E_{v},\quad\bm{\mathcal{E}}_{u}=\prod_{v\in\mathbb{T}_{u}}\mathcal{E}_{v},\quad\forall
u\in\mathbb{T},$

where $\mathbb{T}_{u}$ denotes the sub-tree of $\mathbb{T}$ rooted at $u$
(obtained by removing all nodes from $\mathbb{T}$ except for $u$ and its
descendants).

Similarly to Algorithm 2, we introduce a collection
$(\gamma_{u_{-}})_{u\in\mathbb{T}^{\not\partial}}$ of _auxiliary measures_
indexed by the nodes of our tree that are not leaves. In particular, for any
non-leaf $u$ with set of children $\mathcal{C}_{u}$, $\gamma_{u_{-}}$ is a
finite measure on the corresponding intermediate product space:

$\displaystyle\bm{E}_{\mathcal{C}_{u}}:=\prod_{v\in\mathbb{T}_{u}^{\not
u}}E_{v},\quad\bm{\mathcal{E}}_{\mathcal{C}_{u}}:=\prod_{v\in\mathbb{T}_{u}^{\not
u}}\mathcal{E}_{v},$

where $\mathbb{T}_{u}^{\not u}:=\mathbb{T}_{u}\backslash\\{u\\}$. These
auxiliary distributions fulfil essentially the same role as the auxiliary
distributions in the ASMC algorithm: introducing as much information as
possible about the measure $\rho_{u}$ associated with the current node $u$
_before_ resampling. As we will see in the next section, this mechanism is
particularly important in the divide and conquer context because it allows us
to correlate the components of the resampled particles associated with $u$’s
children.

######  Remark 1 (A connection with SMC).

If the tree is a line (i.e. each node has a single child, except for the first
which has none), then our setting here reduces to that of standard ASMC
(Algorithm 2). If we choose auxiliary distributions
$\gamma_{t_{-}}:=\rho_{t-1}$, then it further reduces to a standard SMC
(Algorithm 1).

For the reasons discussed in Section 4.2, we would like to choose the
intermediate targets and auxiliary measures to be as close as possible to the
corresponding marginals of the final target $\rho_{\mathfrak{r}}$. Of course,
in general, this is more easily said than done and, in practice, $\rho_{u}$
and $\gamma_{u_{-}}$ will often be obtained by ignoring [41, Section 3.4],
rather than ‘integrating-out’, the terms in $\rho_{\mathfrak{r}}$ featuring
variables associated with $u$’s ancestors (and $u$ itself in the case of
$\gamma_{u_{-}}$). To illustrate some of these concepts, consider the
following simplified version of the example in [41, Section 5.2]:

######  Example 1 (Simplified mathematics test dataset).

Suppose we have the following data: the number $M_{sy}$ of students who took a
certain exam in school $s$ and year $y$ and the number $m_{sy}$ which passed
the exam, for all years $y$ that the exam was taken in school $s$ and all
schools $s$ in a given city (we denote the set of the former by
$\mathcal{Y}_{s}$ and that of the latter by $\mathcal{S}$). A standard
hierarchical model [26] for inference on such a dataset is as follows:
$m_{sy}$ is binomially distributed with success probability
$\alpha(\theta_{sy})$ (i.e. $m_{sy}$ has law
$\mathcal{B}(\cdot;M_{sy},\alpha(\theta_{sy}))$), where $\theta_{sy}$ denotes
a latent parameter and $\alpha$ denotes the standard logistic function. The
parameters are related as follows:

$\theta_{sy}=\theta_{s}+\varepsilon_{sy}\enskip\forall
y\in\mathcal{Y}_{s},\enskip
s\in\mathcal{S},\quad\theta_{s}=\theta_{\mathfrak{r}}+\varepsilon_{s}\enskip\forall
s\in\mathcal{S},$

where $(\theta_{s})_{s\in\mathcal{S}}$ and $\theta_{\mathfrak{r}}$ denote
further latent parameters and
$((\varepsilon_{sy})_{y\in\mathcal{Y}_{s}},\varepsilon_{s})_{s\in\mathcal{S}}$
a collection of independent noises with
$\varepsilon_{sy}\sim\mathcal{N}(0,\sigma_{s}^{2})$ and
$\varepsilon_{s}\sim\mathcal{N}(0,\sigma^{2}_{\mathfrak{r}})$ for some unknown
variances $(\sigma_{s}^{2})_{s\in\mathcal{S}}$ and
$\sigma^{2}_{\mathfrak{r}}$. For each parameter variable, we choose as a prior
the pullback $g(d\theta)=\alpha(\theta)(1-\alpha(\theta))d\theta$ of the
uniform distribution on $(0,1)$ via $\alpha$; while for each variance
variable, we choose a unit-mean exponential distribution $f(d\sigma^{2};1)$.
The full (unnormalized) posterior is then given by

$\displaystyle
f(d\sigma^{2}_{\mathfrak{r}};1)g(d\theta_{\mathfrak{r}})\prod_{s\in\mathcal{S}}\Big{[}$
$\displaystyle\mathcal{N}(\theta_{s}-\theta_{\mathfrak{r}};0,\sigma^{2}_{\mathfrak{r}})f(d\sigma^{2}_{s};1)g(d\theta_{s})$
(5)
$\displaystyle\times\prod_{y\in\mathcal{Y}_{s}}\mathcal{N}(\theta_{sy}-\theta_{s};0,\sigma_{s}^{2})\mathcal{B}(m_{sy};M_{sy},\alpha(\theta_{sy}))g(d\theta_{sy})\Big{]}.$

This distribution’s dependence structure is encoded by a tree $\mathbb{T}$
whose paths from root to leaves have the following form: root, school, year.
In other words, $\mathbb{T}$’s root node $\mathfrak{r}$ has a child for each
school ($\mathcal{C}_{\mathfrak{r}}=\mathcal{S}$), each of these has a child
for each year the exam took place in that school
($\mathcal{C}_{s}=\mathcal{Y}_{s}$ for all $s$ in $\mathcal{S}$), and the
latter are all leaves (note that we treat elements in different
$\mathcal{Y}_{s}$s as distinct nodes, even if they represent the same year).
With each leaf node, we associate the corresponding latent variable (so that
$E_{u}:=\mathbb{R}$ if $u$ is a leaf), and, with each non-leaf node, we
associate the corresponding parameter variable and unknown variance (so that
$E_{u}:=\mathbb{R}\times[0,\infty)$ if $u$ is not a leaf).

We set the final target $\rho_{\mathfrak{r}}$ to be our measure of interest
(5). Analytical expressions for its marginals are unavailable, so we instead
pick the intermediate targets as follows: we obtain $\rho_{u}$ by removing
from (5) all terms featuring variables not associated with $u$ or its
descendants:

$\displaystyle\rho_{s}(d\theta_{s},d\sigma^{2}_{s},(d\theta_{sy})_{y\in\mathcal{Y}_{s}})$
$\displaystyle:=f(d\sigma^{2}_{s};1)g(d\theta_{s})$ (6)
$\displaystyle\quad\times\prod_{y\in\mathcal{Y}_{s}}\mathcal{N}(\theta_{sy}-\theta_{s};0,\sigma_{s}^{2})\mathcal{B}(m_{sy};M_{sy},\alpha(\theta_{sy}))g(d\theta_{sy}),$
$\displaystyle\rho_{sy}(d\theta_{sy})$
$\displaystyle:=\mathcal{B}(m_{sy};M_{sy},\alpha(\theta_{sy}))g(d\theta_{sy})\quad\forall
y\in\mathcal{Y}_{s},$

for all $s$ in $\mathcal{S}$. It is also not possible to analytically
integrate over $(\theta_{\mathfrak{r}},\sigma^{2}_{\mathfrak{r}})$ in (5) or
$(\theta_{s},\sigma^{2}_{s})$ in (6), so we instead pick the auxiliary
measures by further doing away with the terms in (5,6) involving the variables
associated with the indexing node:

$\gamma_{\mathfrak{r}_{-}}:=\prod_{s\in\mathcal{S}}\rho_{s},\quad\gamma_{s_{-}}:=\prod_{y\in\mathcal{Y}_{s}}\rho_{sy}\enskip\forall
s\in\mathcal{S}.$

We finish this section with some notation that we use throughout the rest of
the paper:

  * •

(Trees) By a ‘tree’ $\mathbb{T}$ we mean a rooted tree: a connected directed
acyclic graph in which no node has more than one parent and possessing a
single root node, $\mathfrak{r}$, that has no parent. We denote the number of
children a node $u$ has as $c_{u}$, refer to them as $u1,u2,\dots,uc_{u}$, and
denote the set $\\{u1,\dots,uc_{u}\\}$ of them by $\mathcal{C}_{u}$. For any
subset $A$ of the tree $\mathbb{T}$, we use $\left|A\right|$ to denote its
cardinality, $A^{\partial}:=\\{v\in A:c_{v}=0\\}$ to denote the set of all
leaves in $A$, and $A^{\not\partial}:=A\backslash A^{\partial}$ that of all
other nodes in $A$. In particular, $\mathbb{T}^{\partial}$ denotes the set of
$\mathbb{T}$’s leaves and $\mathbb{T}^{\not\partial}$ that of all non-leaf
nodes.

  * •

(Products) We assign the subscript ‘A’ to a symbol, where $A$ is any subset of
$\mathbb{T}$, if the object represented by the symbol depends on the nodes in
$A$. If the symbol is bold, the object also depends on the descendants of the
nodes in $A$. With few exceptions, this dependence takes the form of a
product. For example,

$\displaystyle E_{A}:=$ $\displaystyle\prod_{u\in A}E_{u},$
$\displaystyle\bm{E}_{A}:=$ $\displaystyle E_{\cup_{u\in A}\mathbb{T}_{u}},$
$\displaystyle\mathcal{E}_{A}:=$ $\displaystyle\prod_{u\in A}\mathcal{E}_{u},$
$\displaystyle\bm{\mathcal{E}}_{A}:=$ $\displaystyle\mathcal{E}_{\cup_{u\in
A}\mathbb{T}_{u}},$ $\displaystyle\rho_{A}:=$ $\displaystyle\prod_{u\in
A}\rho_{u},$ $\displaystyle\mu_{A}:=$ $\displaystyle\prod_{u\in A}\mu_{u},$
$\displaystyle Z_{A}:=$ $\displaystyle\prod_{u\in A}Z_{u}.$

for any non-empty subset $A$ of $\mathbb{T}$. The exceptions are $x_{A}$ and
$\bm{x}_{A}$ which denote elements in $E_{A}$ and $\bm{E}_{A}$ rather than any
sort of product over $A$.

###  3.2 The algorithm

To introduce the divide-and-conquer SMC algorithm, we require one last
ingredient: the generalization of the proposal kernels used in SMC/ASMC (cf.
Section 2.2) to extend the path of the resampled particles. In our case, we
have, for each non-leaf node $u$, one Markov kernel
$K_{u}:\bm{E}_{\mathcal{C}_{u}}\times\mathcal{E}_{u}\to[0,1]$ mapping from the
product $\bm{E}_{\mathcal{C}_{u}}=\prod_{v\in\mathcal{C}_{u}}\bm{E}_{v}$ of
the intermediate spaces $\bm{E}_{v}$ indexed by $u$’s children to the low-
dimensional space $E_{u}$ indexed by $u$ itself. For the reasons discussed in
Section 4.2, the _proposal kernel_ $K_{u}$ should be chosen as close as
possible to the conditional distribution over $E_{u}$ under $\mu_{u}$ given
values for its remaining coordinates, see for details.

Algorithm 3 describes DaC-SMC111Note that, in the original paper [41],
Algorithm 3 was referred to as ‘DaC-SMC with mixture resampling’ while the
term ‘DaC-SMC’ was reserved for cases where the auxiliary measures are
product-form and the algorithm’s running time can be sped up through careful
implementation (cf. Section 3.3). For all, at least theoretical, purposes the
latter is a special case of the former and we simply refer to both as the
‘DaC-SMC Algorithm’. recursively: the first thing it does when called on any
given non-leaf node is to call itself on each of that node’s children. Hence,
the algorithm proper begins ‘at the leaves’ by drawing, for each leaf $u$, $N$
samples $\bm{X}_{u}^{1,N},\dots,\bm{X}_{u}^{N,N}$ from a proposal distribution
$K_{u}$ on $(E_{u},\mathcal{E}_{u})$ and setting
$\gamma_{u}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}_{u}^{n,N}}$ to be the
corresponding unweighted sample approximation to $K_{u}$. It then works its
way from leaves to root, at each node combining and augmenting the
approximations indexed by its children.

1: if $u$ is a leaf (i.e. $u\in\mathbb{T}^{\partial}$) then

2: Propose: for $n\leq N$, draw $\bm{X}_{u}^{n,N}$ independently from $K_{u}$.

3: Return: $\gamma_{u}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}_{u}^{n,N}}$.

4: else

5: for $v$ in $\mathcal{C}_{u}$ do

6: Recurse: set $\gamma_{v}^{N}:=\text{dac\\_smc}(v)$.

7: end for

8: Correct: compute $\gamma_{u_{-}}^{N}$ in (8) and
$\pi^{N}_{u_{-}}:=\gamma_{u_{-}}^{N}/\gamma_{u_{-}}^{N}(\bm{E}_{\mathcal{C}_{u}})$.

9: Resample: for $n\leq N$, draw $\bm{X}_{u_{-}}^{n,N}$ independently from
$\pi_{u_{-}}^{N}$.

10: Mutate: for $n\leq N$, draw $X^{n,N}_{u}$ independently from
$K_{u}(\bm{X}_{u_{-}}^{n,N},dx_{u})$ and set
$\bm{X}_{u}^{n,N}:=(X^{n,N}_{u},\bm{X}_{u_{-}}^{n,N})$.

11: Return:
$\gamma_{u}^{N}:=N^{-1}\mathcal{Z}_{u}^{N}\sum_{n=1}^{N}\delta_{\bm{X}_{u}^{n,N}}$
where $\mathcal{Z}_{u}^{N}:=\gamma_{u_{-}}^{N}(\bm{E}_{\mathcal{C}_{u}})$.

12: end if
Note: One can use (11) to compute estimates of $\rho_{u}$, $Z_{u}$, and
$\mu_{u}$.

Algorithm 3 dac_smc$(u)$ for $u$ in $\mathbb{T}$.

More precisely, at any non-leaf node $u$, and for each child $v$ of $u$, the
algorithm first retrieves an unweighted particle approximation
$\gamma_{v}^{N}=\mathcal{Z}_{v}^{N}N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}_{v}^{n,N}}$
(of mass $\mathcal{Z}_{v}^{N}$) to the _extended auxiliary measure_
$\gamma_{v}$ obtained by taking the outer product of the auxiliary measure
$\gamma_{v_{-}}$ and the proposal kernel $K_{v}$:

$\displaystyle\gamma_{v}(d\bm{x}_{v}):=(\gamma_{v_{-}}\otimes
K_{v})(d\bm{x}_{v})=\gamma_{v_{-}}(d\bm{x}_{\mathcal{C}_{v}})K_{v}(\bm{x}_{\mathcal{C}_{v}},dx_{v})\quad\forall
v\in\mathbb{T}^{\not\partial}.$ (7)

(If $v$ is a leaf then the algorithm just retrieves the sample approximation
to $\gamma_{v}:=K_{v}$ mentioned in the preceding paragraph.) It then takes
the product
$\gamma_{\mathcal{C}_{u}}^{N}:=\prod_{v\in\mathcal{C}_{u}}\gamma_{v}^{N}$ of
these and re-weighs it using
$w_{u_{-}}:=d{\gamma_{u_{-}}}/d\gamma_{\mathcal{C}_{u}}$ to obtain a product-
form estimator (cf. (1)) for the auxiliary distribution $\gamma_{u_{-}}$
associated with node $u$:

$\gamma_{u_{-}}^{N}(d\bm{x}_{\mathcal{C}_{u}}):=w_{u_{-}}(\bm{x}_{\mathcal{C}_{u}})\gamma_{\mathcal{C}_{u}}^{N}(d\bm{x}_{\mathcal{C}_{u}})=\frac{\mathcal{Z}^{N}_{\mathcal{C}_{u}}}{N^{c_{u}}}\sum_{\bm{n}\in[N]^{c_{u}}}w_{u_{-}}(\bm{X}_{\mathcal{C}_{u}}^{\bm{n},N})\delta_{\bm{X}_{\mathcal{C}_{u}}^{\bm{n},N}}(d\bm{x}_{\mathcal{C}_{u}})$
(8)

where
$\mathcal{Z}^{N}_{\mathcal{C}_{u}}:=\prod_{v\in\mathcal{C}_{u}}\mathcal{Z}_{v}^{N}$
and, for each $\bm{n}=(n_{1},\dots,n_{c_{u}})$ in $[N]^{c_{u}}$,
$\bm{X}^{\bm{n},N}_{\mathcal{C}_{u}}$ denotes the ‘permuted particle’
$(\bm{X}^{n_{1}}_{u1},\dots,\bm{X}^{n_{c_{u}}}_{uc_{u}})$. Just as in the
standard ASMC case (Algorithm 2), the algorithm proceeds by normalizing
$\gamma_{u_{-}}^{N}$ to obtain an approximation $\pi_{u_{-}}^{N}$ of the
_normalized auxiliary measure_ $\pi_{u_{-}}$,

$\pi_{u_{-}}^{N}:=\frac{\gamma_{u_{-}}^{N}}{\mathcal{Z}^{N}_{u}}\approx\frac{\gamma_{u_{-}}}{\mathcal{Z}_{u}}=:\pi_{u_{-}}\enskip\text{where}\enskip\mathcal{Z}^{N}_{u}:=\gamma_{u_{-}}^{N}(\bm{E}_{\mathcal{C}_{u}}),\enskip\mathcal{Z}_{u}:=\gamma_{u_{-}}(\bm{E}_{\mathcal{C}_{u}}),$
(9)

and it draws $N$ particles $\bm{X}_{u_{-}}^{1,N},\dots,\bm{X}_{u_{-}}^{N,N}$
from $\gamma_{u_{-}}^{N}$ using (multinomial) resampling. The empirical
distribution $N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}_{u_{-}}^{n,N}}$ of these
particles approximates $\pi_{u_{-}}$. Hence, by extending the path of each
resampled particle from $\bm{E}_{\mathcal{C}_{u}}$ to $\bm{E}_{u}$ using the
kernel $K_{u}$,

$\bm{X}_{u}^{n,N}:=(X^{n,N}_{u},\bm{X}_{u_{-}}^{n,N})\enskip\text{where}\enskip
X^{n,N}_{u}\sim K_{u}(\bm{X}_{u_{-}}^{n,N},dx_{u}),$

and defining

$\displaystyle\gamma_{u}^{N}(d\bm{x}_{u}):=$
$\displaystyle\frac{\mathcal{Z}_{u}^{N}}{N}\sum_{n=1}^{N}\delta_{\bm{X}_{u}^{n,N}}\enskip\text{and}\enskip\pi_{u}^{N}(d\bm{x}_{u}):=\frac{1}{N}\sum_{n=1}^{N}\delta_{\bm{X}_{u}^{n,N}},$
(10)

Algorithm 3 returns a finite-sample approximation of the extended auxiliary
measure $\gamma_{u}:=\gamma_{u_{-}}\otimes K_{u}$ indexed by $u$ itself.

Also similarly to the standard ASMC case (Section 2.2), we obtain
approximations of to the targets $\rho_{u}$ and $\mu_{u}$ by running
dac_smc($u$) in Algorithm 3 and applying a simple importance sampling
correction:

$\displaystyle\rho_{u}^{N}(d\bm{x}_{u})=w_{u}(\bm{x}_{u})\gamma_{u}^{N}(d\bm{x}_{u}),\quad
Z_{u}^{N}=\rho_{u}^{N}(\bm{E}_{u}),\quad\mu_{u}^{N}(d\bm{x}_{u})=\frac{\rho_{u}^{N}(d\bm{x}_{u})}{Z_{u}^{N}},$
(11)

where $w_{u}:=d\rho_{u}/d\gamma_{u}$ denotes the appropriate weight function.
As we will see in Section 4, these approximations have appealing theoretical
properties. The computation time of Algorithm 3 can be lowered by
parallelizing the for loop in lines 5–7. Of course, for the approach to work,
we must choose the auxiliary measures
$(\gamma_{u_{-}})_{u\in\mathbb{T}^{\not\partial}}$ and proposal kernels
$(K_{u})_{u\in\mathbb{T}}$ such that the relevant Radon-Nikodym derivatives
exist. To avoid the technical machinery necessary to allow for the possibility
of all particles simultaneously being assigned zero weight (cf. [17, Section
7.4]), we further assume that these derivatives are positive everywhere:

######  Assumption 1.

For all $u$ in $\mathbb{T}$ and $v$ in $\mathbb{T}^{\not\partial}$, $\rho_{u}$
is absolutely continuous w.r.t. $\gamma_{u}$, $\gamma_{v_{-}}$ is absolutely
continuous w.r.t. $\gamma_{\mathcal{C}_{v}}$, and the Radon-Nikodym
derivatives $w_{u}:=d\rho_{u}/d\gamma_{u}$ and
$w_{v_{-}}:=d\gamma_{v_{-}}/d\gamma_{\mathcal{C}_{v}}$ are positive
everywhere.

To finish this section, we illustrate how one may choose the proposal kernels
and obtain the weight functions by revisiting the example of the previous
section:

######  Example 2 (Simplified mathematics test dataset, kernels and weights).

To apply DaC-SMC to Example 1, we need to pick proposal kernels. In the case
of the leaves, the choice is clear: for all $s,y,n$, we generate
$\theta_{sy}^{n,N}$ by drawing a sample from a Beta distribution with
parameters $m_{sy}+1$ and $M_{sy}-m_{sy}+1$ and mapping the sample through the
logit function; so that

$K_{sy}(d\theta_{sy})=(M_{sy}+1)\binom{M_{sy}}{m_{sy}}\mathcal{B}(m_{sy};M_{sy},\alpha(\theta_{sy}))g(d\theta_{sy})$

and the weight functions are constant (for all $s$ in $\mathcal{S}$,
$w_{s_{-}}=\prod_{y\in\mathcal{Y}_{s}}M_{sy}^{-1}$ and, for all $y$ in
$\mathcal{Y}_{s}$, $w_{sy}=M_{sy}^{-1}$). In the case of the other nodes, we
would ideally set the kernel to be the appropriate conditional distribution
(cf. Section 4.2),

$\displaystyle
H_{s}((\theta_{sy})_{y\in\mathcal{Y}_{s}},d\theta_{s},d\sigma^{2}_{s})$
$\displaystyle\propto
f(d\sigma^{2}_{s};1)g(d\theta_{s})\prod_{y\in\mathcal{Y}_{s}}\mathcal{N}(\theta_{sy}-\theta_{s};0,\sigma_{s}^{2})\quad\forall
s\in\mathcal{S};$ $\displaystyle
H_{\mathfrak{r}}((\theta_{s},\sigma^{2}_{s})_{s\in\mathcal{S}},d\theta_{\mathfrak{r}},d\sigma^{2}_{\mathfrak{r}})$
$\displaystyle\propto
f(d\sigma^{2}_{\mathfrak{r}};1)g(d\theta_{\mathfrak{r}})\prod_{s\in\mathcal{S}}\mathcal{N}(\theta_{s}-\theta_{\mathfrak{r}};0,\sigma^{2}_{\mathfrak{r}}).$

These, however, are intractable. So, in an effort to approximate them, we note
that

$\prod_{s\in\mathcal{S}}\mathcal{N}(\theta_{s}-\theta_{\mathfrak{r}};0,\sigma^{2}_{\mathfrak{r}})=\frac{e^{\frac{1}{2\sigma^{2}}\left[\sum_{s\in\mathcal{S}}\theta_{s}\left(\theta_{s}-\frac{1}{\left|\mathcal{S}\right|}\sum_{s^{\prime}\in\mathcal{S}}\theta_{s^{\prime}}\right)\right]}}{\sqrt{\left|\mathcal{S}\right|}\sqrt{2\pi\sigma^{2}_{\mathfrak{r}}}^{\left|S\right|-1}}\mathcal{N}\left(d\theta_{\mathfrak{r}};\frac{1}{\left|\mathcal{S}\right|}\sum_{s\in\mathcal{S}}\theta_{s},\frac{\sigma^{2}_{\mathfrak{r}}}{\left|\mathcal{S}\right|}\right)$

and similarly for
$\prod_{y\in\mathcal{Y}_{s}}\mathcal{N}(\theta_{sy}-\theta_{s};0,\sigma_{s}^{2})$
and all $s$ in $\mathcal{S}$; and we instead set

$\displaystyle
K_{s}((\theta_{sy})_{y\in\mathcal{Y}_{s}},d\theta_{s},d\sigma^{2}_{s})$
$\displaystyle:=f(d\sigma^{2}_{s};1)\mathcal{N}\left(d\theta_{s};\frac{1}{\left|\mathcal{Y}_{s}\right|}\sum_{y\in\mathcal{Y}_{s}}\theta_{sy},\frac{\sigma^{2}_{s}}{\left|\mathcal{Y}_{s}\right|}\right),$
$\displaystyle
K_{\mathfrak{r}}((\theta_{s},\sigma^{2}_{s})_{s\in\mathcal{S}},d\theta_{\mathfrak{r}},d\sigma^{2}_{\mathfrak{r}})$
$\displaystyle:=f(d\sigma^{2}_{\mathfrak{r}};1)\mathcal{N}\left(d\theta_{\mathfrak{r}};\frac{1}{\left|\mathcal{S}\right|}\sum_{s\in\mathcal{S}}\theta_{s},\frac{\sigma^{2}_{\mathfrak{r}}}{\left|\mathcal{S}\right|}\right);$

in which case $w_{\mathfrak{r}_{-}}=\prod_{s\in\mathcal{S}}w_{s}$,

$\displaystyle
w_{s}(\theta_{s},\sigma^{2}_{s},(\theta_{sy})_{y\in\mathcal{Y}_{s}})=\frac{g(\theta_{s})e^{\frac{1}{2\sigma^{2}_{s}}\left[\sum_{y\in\mathcal{Y}_{s}}\theta_{sy}\left(\theta_{sy}-\frac{1}{\left|\mathcal{Y}_{s}\right|}\sum_{y^{\prime}\in\mathcal{Y}_{s}}\theta_{sy^{\prime}}\right)\right]}}{\sqrt{\left|\mathcal{Y}_{s}\right|}\sqrt{2\pi\sigma^{2}_{s}}^{\left|\mathcal{Y}_{s}\right|-1}}\quad\forall
s\in\mathcal{S},$ $\displaystyle
w_{\mathfrak{r}}(\theta_{\mathfrak{r}},\sigma^{2}_{\mathfrak{r}},(\theta_{s},\sigma^{2}_{s},(\theta_{sy})_{y\in\mathcal{Y}_{s}})_{s\in\mathcal{S}})=\frac{g(\theta_{\mathfrak{r}})e^{\frac{1}{2\sigma^{2}_{\mathfrak{r}}}\left[\sum_{s\in\mathcal{S}}\theta_{s}\left(\theta_{s}-\frac{1}{\left|\mathcal{S}\right|}\sum_{s^{\prime}\in\mathcal{S}}\theta_{s^{\prime}}\right)\right]}}{\sqrt{\left|\mathcal{S}\right|}\sqrt{2\pi\sigma^{2}_{\mathfrak{r}}}^{\left|\mathcal{S}\right|-1}}.$

To sample $(\theta_{s},\sigma^{2}_{s})$ from
$K_{s}((\theta_{sy})_{y\in\mathcal{Y}_{s}},d\theta_{s},d\sigma^{2}_{s})$,
first sample $\sigma^{2}_{s}$ from $f(d\sigma^{2}_{s};1)$ and then sample
$\theta_{s}$ from
$\mathcal{N}(d\theta_{s};\left|\mathcal{Y}_{s}\right|^{-1}\sum_{y\in\mathcal{Y}_{s}}\theta_{sy},\left|\mathcal{Y}_{s}\right|^{-1}\sigma^{2}_{s})$.
Similarly for $K_{\mathfrak{r}}$.

###  3.3 Efficient computation

The main drawback of Algorithm 3 is the $\mathcal{O}(N^{c_{u}})$ computational
cost of resampling from the product-form $\pi_{u_{-}}^{N}$ in (9) for a
general auxiliary measure $\gamma_{u_{-}}$ (to draw $N$ samples from
$\pi_{u_{-}}^{N}$ in $\mathcal{O}(N^{c_{u}})$ operations one can, for example,
use the alias method [57]). This results in a total algorithmic cost of
$\mathcal{O}(N^{d})$, where $d$ denotes the tree’s degree (i.e. the largest
number of children that any of the tree’s nodes posses). In some cases, $d$ is
small enough that the cost is manageable. For instance, when the collection of
target distributions $(\rho_{u})_{u\in\mathbb{T}_{u}}$ is an artificial
construct introduced for computational purposes and the interest lies in the
final target, it is common [21, 20, 13] to choose the indexing tree to be
binary, in which case $d=2$ and the algorithm’s cost is $\mathcal{O}(N^{2})$.

Otherwise, one can avoid the $\mathcal{O}(N^{c_{u}})$ cost by choosing
auxiliary measures that factorize. For instance, suppose that, as in Example
1, they factorize fully:

$\gamma_{u_{-}}:=\prod_{v\in\mathcal{C}_{u}}\bar{\gamma}_{v}\quad\forall
u\in\mathbb{T}^{\not\partial},$

where $(\bar{\gamma}_{v})_{v\in\mathbb{T}}$ denotes a collection of measures
on the respective spaces $(\bm{E}_{v},\bm{\mathcal{E}}_{v})_{v\in\mathbb{T}}$.
In this case, the weight function $w_{u_{-}}$ decomposes into the product

$\displaystyle
w_{u_{-}}=\frac{d\gamma_{u_{-}}}{d\gamma_{\mathcal{C}_{u}}}=\prod_{v\in\mathcal{C}_{u}}\frac{d\bar{\gamma}_{v}}{d\gamma_{v}}=:\prod_{v\in\mathcal{C}_{u}}\bar{w}_{v},$

and the correction step (line 8 in Algorithm 3) breaks down into the
following: compute

$\gamma_{u_{-}}^{N}:=\prod_{v\in\mathcal{C}_{u}}\bar{\gamma}_{v}^{N},\quad\pi_{u_{-}}^{N}:=\prod_{v\in\mathcal{C}_{u}}\bar{\pi}_{v}^{N},$

where, for all children $v$ of $u$,

$\bar{\gamma}_{v}^{N}(d\bm{x}_{v}):=\bar{w}_{v}(\bm{x}_{v})\gamma_{v}(d\bm{x}_{v})=\frac{\mathcal{Z}_{v}^{N}}{N}\sum_{n=1}^{N}\bar{w}_{v}(\bm{x}_{v})\delta_{\bm{X}_{v}^{n,N}}(d\bm{x}_{v}).$
(12)

Resampling from
$\pi_{u_{-}}^{N}:=\bar{\gamma}_{v}^{N}/\bar{\gamma}_{v}^{N}(\bm{E}_{v})$ can
then be achieved in $\mathcal{O}(N)$ operations by independently drawing
samples $N$ from $\bar{\pi}_{u1}^{N},\dots,\bar{\pi}_{uc_{u}}^{N}$, e.g. by
applying the alias method [57] to each approximation separately, and
concatenating them. In other words, we can replace lines 5–9 in Algorithm 3
with those in Algorithm 4 (and obtain the algorithm referred to as ‘DaC-SMC’
in [41], cf. Footnote 1). The running time can then be lowered further by
parallelizing these operations across $u$’s children.

1: for $v$ in $\mathcal{C}_{v}$ do

2: Recurse: set $\gamma_{v}^{N}:=\text{dac\\_smc}(v)$.

3: Correct: compute $\bar{\gamma}_{v}^{N}$ in (12) and
$\bar{\pi}^{N}_{v}:=\bar{\gamma}_{v}^{N}/\bar{\gamma}_{v}^{N}(\bm{E}_{v})$.

4: Resample: for $n\leq N$, $\bm{\bar{X}}_{v}^{n,N}$ independently from
$\bar{\pi}_{v}^{N}$.

5: end for

6: Concatenate: for $n\leq N$, set
$\bm{X}_{u_{-}}^{n,N}:=(\bm{\bar{X}}_{v}^{n,N})_{v\in\mathcal{C}_{u}}$.

Algorithm 4 $\mathcal{O}(N)$-cost replacement of lines 5–9 in Algorithm 3 for
fully-factorized $\gamma_{u_{-}}$.

These choices result in DaC-SMC employing extended auxiliary measures with
partial product structure:

$\gamma_{u}(d\bm{x}_{u})=K_{u}(\bm{x}_{\mathcal{C}_{u}},dx_{u})\prod_{v\in\mathcal{C}_{v}}\bar{\gamma}_{v}(d\bm{x}_{v})\quad\forall
u\in\mathbb{T}^{\not\partial}.$

(I.e. if $(Y_{u},\bm{Y}_{u1},\dots,\bm{Y}_{uc_{u}})\sim\pi_{u}$, then
$\bm{Y}_{u1},\dots,\bm{Y}_{uc_{u}}$ are independent.) If the targets
$\rho_{u}$ and $\mu_{u}$ do not share the above structure, then this
$\mathcal{O}(N)$ approach will generally require large
$\gamma_{u}^{N}\mapsto\rho_{u}^{N}$ corrections and, consequently, result in
greater estimator variance (which may, or may not, be compensated by the
computational gains). As suggested in [41, Section 4.2], one way to mitigate
this issue is via tempering [19].

Another way to incorporate non-product structure into the auxiliary measures
without sacrificing the $\mathcal{O}(N)$ cost is setting them to be sums of
fully-factorized measures:

$\gamma_{u_{-}}:=\sum_{i=1}^{I}\prod_{v\in\mathcal{C}_{u}}\bar{\gamma}_{v}^{i}=:\sum_{i=1}^{I}\gamma_{u_{-}}^{i}\quad\forall
u\in\mathbb{T}^{\not\partial},$

where $(\bar{\gamma}_{v}^{1},\dots,\bar{\gamma}_{v}^{I})_{v\in\mathbb{T}}$
denote collections of measures on the respective spaces
$(\bm{E}_{v},\bm{\mathcal{E}}_{v})_{v\in\mathbb{T}}$. In this case, the
correction step becomes: return

$\gamma_{u_{-}}^{N}:=\sum_{i=1}^{I}\gamma_{u_{-}}^{i,N},\quad\pi_{u_{-}}^{N}:=\sum_{i=1}^{I}\frac{\gamma_{u_{-}}^{i,N}(\bm{E}_{\mathcal{C}_{u}})}{\gamma_{u_{-}}^{N}(\bm{E}_{\mathcal{C}_{u}})}\pi_{u_{-}}^{i,N}=:\sum_{i=1}^{I}\omega_{u_{-}}^{i,N}\pi_{u_{-}}^{i,N},$

where
$\gamma_{u_{-}}^{i,N}:=\prod_{v\in\mathcal{C}_{u}}\bar{\gamma}_{v}^{i,N}$ and
$\pi_{u_{-}}^{i,N}:=\prod_{v\in\mathcal{C}_{u}}\bar{\pi}_{v}^{i,N}$, with
$\bar{\gamma}_{v}^{i,N},\bar{\pi}_{v}^{i,N}$ defined analogously to
$\bar{\gamma}_{v}^{N},\bar{\pi}_{v}^{N}$ in (12). We can then draw $N$ samples
from $\pi_{u_{-}}^{N}$ in $\mathcal{O}(N)$ operations by, for instance,
drawing $(m_{1},\dots,m_{I})$ from a multinomial distribution with weights
$(\omega_{u_{-}}^{i,N})_{i=1}^{I}$ and, for each $i$, generating $m_{i}$
samples from $\pi_{u_{-}}^{i,N}$ as described above for the $I=1$ case.

Similar considerations apply if the auxiliary measures are chosen to be sums
of partially factorized functions rather than fully factorized ones. The main
difference is that resampling from $\pi_{u_{-}}^{N}$ has a ‘conditional
aspect’ to it and costs $\mathcal{O}(N^{a})$ where the exponent $1<a<c_{u}$
depends on the amount of factorization, see Appendix A.1 for an example.

###  3.4 Variants

To simplify the algorithm’s presentation and analysis, we focus throughout on
the particular version given in Algorithm 3. There are, however, many variants
that one might wish to consider in different settings and whose analysis may
be tackled using straightforward modifications of the arguments given in this
paper.

First off, we have the low-cost variant referred to as ‘lightweight mixture
resampling’ in [41, Section 4.1] (see ‘multiple matching’ in [40, Section 2.2]
for similar ideas) where the product-form estimator
$\gamma_{u_{-}}^{N}(d\bm{x}_{\mathcal{C}_{u}})$ in (8) is replaced with its
‘incomplete’ version,

$\frac{\mathcal{Z}^{N}_{\mathcal{C}_{u}}}{\left|\mathcal{M}\right|}\sum_{\bm{n}\in\mathcal{M}}w_{u_{-}}(\bm{X}_{\mathcal{C}_{u}}^{\bm{n},N})\delta_{\bm{X}_{\mathcal{C}_{u}}^{\bm{n},N}}(d\bm{x}_{\mathcal{C}_{u}}),$
(13)

with $\mathcal{M}\subseteq[N]^{c_{u}}$ indexing a user-chosen subset of
permuted particles $\bm{X}_{\mathcal{C}_{u}}^{\bm{n},N}$. In particular, by
picking $\mathcal{M}$ carefully, it might be possible to substantially reduce
the cost without sacrificing too much of the variance reduction (cf. [35] and
references therein for similar feats in the U-statistics literature). While
one may opt for more sophisticated approaches (e.g. ones along the lines of
those in [35]), we can offer two simple rules of thumb for choosing
$\mathcal{M}$ that are motivated by the considerations in [37, Section 4]: (a)
set $\mathcal{M}$’s size to be such that the cost of resampling from (13)’s
normalization is comparable to that of generating the unpermuted particles,
$(\bm{X}_{u1}^{\bm{n},N},\dots,\bm{X}_{uc_{u}}^{\bm{n},N})_{n=1}^{N}$; and (b)
minimize the number of over-lapping components in the elements of
$\mathcal{M}$. In terms of theoretical analysis for this variant, one can
combine the methods in Appendices C–G with the techniques used to study
incomplete U-statistics [39, Chapter 4.3] and their generalizations.

We also need not resample at every step and may instead choose to do so
adaptively (e.g. resample only when the effective sample size drops underneath
a predetermined threshold [34]). To prevent the number of particles from
blowing up, any product-form approximation $\pi_{u_{-}}^{N}$ to a normalized
auxiliary distribution $\pi_{u_{-}}$ that is not resampled must be replaced
with an $N$ sample one, e.g.

$\frac{\sum_{n=1}^{N}w_{u_{-}}(\bm{X}_{u1}^{n,N},\dots,\bm{X}_{uc_{u}}^{n,N})\delta_{(\bm{X}_{u1}^{n,N},\dots,\bm{X}_{uc_{u}}^{n,N})}}{\sum_{n=1}^{N}w_{u_{-}}(\bm{X}_{u1}^{n,N},\dots,\bm{X}_{uc_{u}}^{n,N})}$
(14)

as in standard ASMC. We should point out here that, for DaC-SMC, and in
contrast with SMC/ASMC, there are benefits to resampling beyond just
mitigating weight degeneracy. In particular, by permuting and resampling the
particles, we explore areas of the target space that would be missed were we
to only use the $N$ unpermuted particles [37, Section 3.1].

There is also flexibility in the resampling scheme employed: the theoretical
analysis in Section 4 applies to Algorithm 3 in which multinomial resampling
is explicitly encoded. But one could extend our analysis along the lines of
[27] to cover a much broader class of resampling schemes, including
substantially lower-variance ones. An interesting question of practical
importance is how best to efficiently implement low-variance resampling
schemes when the auxiliary measures decompose into sums of partially-
factorized measures (e.g. consider the example in Appendix A.1).

One may also wish to implement minor variations to Algorithm 3 dictated by the
particular target of interest. For instance, in some cases, it might be much
easier to specify the auxiliary measure $\gamma_{u-}$ over the variables
taking values in $\bm{E}_{u}$ than over those taking values in
$\bm{E}_{\mathcal{C}_{u}}$. In this case, mutation would need to be carried
out before correction and resampling (line 10 in Algorithm 3 before lines 8–9)
leading to a slightly different algorithm with essentially the same
properties. This approach allows one to specify an analogue of the extended
auxiliary distribution independently of the proposal kernel and directly
exploit the information it encodes in the resampling step. The cost of doing
so is having to perform mutation $N^{c_{u}}$ times, rather than $N$ times,
although this does not increase the overall complexity of the basic algorithm
(at least for unfactorized auxiliary meaures).

Along similar lines, for targets with particular dependence structures, one
may use generalizations of product-form estimators (e.g. [37, Section 3.2 and
Appendix F]) that directly account for these structures; see Appendix A.2 for
an example. The idea is that, by doing so, one can ease the burden placed on
the correction steps which, in Algorithm 3, are the sole responsibles for
introducing dependencies among children into the approximations.

Lastly, one need not generate the same number of particles at each node and
may instead opt to allocate more computational power to the nodes that prove
most problematic (e.g. the $u$s whose associated space $E_{u}$ have greatest
dimension), or even take a more sophisticated approach of the sort in [38].
One would just have to replace $\gamma_{u}^{N}$ in (8) with

$\left(\prod_{v\in\mathcal{C}}\frac{\mathcal{Z}^{N_{v}}_{v}}{N_{v}}\right)\sum_{n_{1}=1}^{N_{u1}}\dots\sum_{n_{1}=1}^{N_{uc_{u}}}w_{u_{-}}(\bm{X}_{u1}^{n_{1},N_{u1}},\dots,\bm{X}_{uc_{u}}^{n_{c_{u}},N_{uc_{u}}})\delta_{(\bm{X}_{u1}^{n_{1},N_{u1}},\dots,\bm{X}_{uc_{u}}^{n_{c_{u}},N_{uc_{u}}})}(d\bm{x}_{\mathcal{C}_{u}})$

where, for each child $v$ of $u$, $N_{v}$ denotes the amount of particles
resampled at node $v$.

##  4 Theoretical characterization

We now turn to the main results of the paper showing that Algorithm 3 is well-
founded. We give an overview of the proofs for these results in Section 4.1,
postponing the full details until Appendices C–G. To simplify the exposition,
we only state results for the particle approximations to the targets; however,
analogous statements for the approximations to the auxiliary measures and
their extensions can be extracted from Appendices C–F. Also to keep the
exposition simple, we focus throughout on the case of bounded test functions,
writing $\mathcal{B}_{b}(S)$ for the space of bounded measurable real-valued
test functions on a measurable space $(S,\mathcal{S})$ and
$\left|\left|{\varphi}\right|\right|:=\sup_{x\in
S}\left|\varphi(x)\right|<\infty$ for the supremum norm on
$\mathcal{B}_{b}(S)$, and we assume that the weight functions are bounded:

######  Assumption 2.

For all $u$ in $\mathbb{T}^{\not\partial}$ and $v$ in $\mathbb{T}$,
$w_{u_{-}}=d\gamma_{u_{-}}/d\gamma_{\mathcal{C}_{u}}$ and
$w_{v}=d\rho_{v}/d\gamma_{v}$ are bounded: $||w_{u_{-}}||<\infty$ and
$\left|\left|{w_{v}}\right|\right|<\infty$.

To begin with, Algorithm 3 produces strongly consistent estimators for the
targets:

######  Theorem 1 (Strong laws of large numbers).

If Assumptions 1–2 are satisfied, $u$ belongs to $\mathbb{T}$, and $\varphi$
belongs to $\mathcal{B}_{b}(\bm{E}_{u})$, then

$\displaystyle\lim_{N\to\infty}\rho^{N}_{u}(\varphi)=\rho_{u}(\varphi),\quad\lim_{N\to\infty}\mu^{N}_{u}(\varphi)=\mu(\varphi),\quad\lim_{N\to\infty}Z^{N}_{u}=Z_{u},\quad\text{almost
surely.}$

###### Proof.

See Appendix C. ∎

If the underlying spaces posses nice enough topological properties, the above
almost sure pointwise convergence can be strengthened to almost sure weak
convergence:

######  Theorem 2 (Almost sure weak convergence).

If, in addition to Assumptions 1–2, the spaces $(E_{u})_{u\in\mathbb{T}}$ are
Polish and $(\mathcal{E}_{u})_{u\in\mathbb{T}}$ are the corresponding Borel
sigma algebras, then

$\displaystyle\rho_{u}^{N}\rightharpoonup\rho_{u},\quad\mu_{u}^{N}\rightharpoonup\mu_{u},\quad\text{almost
surely},$

for each $u$ in $\mathbb{T}$, where $\rightharpoonup$ denotes weak convergence
as $N\to\infty$.

###### Proof.

See Appendix C.∎

The estimators for the unnormalized targets are unbiased:

######  Theorem 3 (Unbiasedness of $(\rho_{u}^{N})_{u\in\mathbb{T}}$).

If Assumptions 1–2 hold, then

$\displaystyle\mathbb{E}\left[\rho^{N}_{u}(\varphi)\right]=\rho_{u}(\varphi),\quad\mathbb{E}\left[Z^{N}_{u}\right]=Z_{u},\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}),\enskip u\in\mathbb{T}.$

###### Proof.

See Appendix D. ∎

The non-linearity in the correction step introduces a bias in the estimators
for the normalized targets. However, the bias decays linearly with the number
of particles $N$:

######  Theorem 4 (Bias estimates for $(\mu_{u}^{N})_{u\in\mathbb{T}}$).

If Assumptions 1–2 hold, the weight functions are bounded below (i.e., for all
$u$ in $\mathbb{T}^{\not\partial}$ and $v$ in $\mathbb{T}$,
$w_{u_{-}}\geq\beta_{u_{-}}$ and $w_{v}\geq\beta_{v}$ for some constants
$\beta_{u_{-}},\beta_{v}>0$), and $u$ belongs to $\mathbb{T}$, then there
exists a constant $C_{u}<\infty$ such that

$\left|\mathbb{E}\left[\mu^{N}_{u}(\varphi)\right]-\mu_{u}(\varphi)\right|\leq\frac{C_{u}\left|\left|{\varphi}\right|\right|}{N}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}).$

###### Proof.

See Appendix E. ∎

All estimators converge at a rate proportional the square root of the number
of particles:

######  Theorem 5 ($L^{p}$ inequalities).

If Assumptions 1–2 hold, then, for each $p\geq 1$ and $u$ in $\mathbb{T}$,
then there exists constants $C_{u}^{\rho},C_{u}^{\mu}<\infty$ such that

$\displaystyle\mathbb{E}\left[\left|\rho^{N}_{u}(\varphi)-\rho(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C_{u}^{\rho}\left|\left|{\varphi}\right|\right|}{N^{1/2}},\quad\mathbb{E}\left[\left|\mu^{N}_{u}(\varphi)-\mu_{u}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C_{u}^{\mu}\left|\left|{\varphi}\right|\right|}{N^{1/2}},$

for all $N>0$ and $\varphi$ in $\mathcal{B}_{b}(\bm{E}_{u})$. In particular,
$\mathbb{E}\left[\left|Z^{N}_{u}-Z_{u}\right|^{p}\right]^{1/p}\leq
C_{u}^{\rho}/N^{1/2}$ for all $N>0$.

###### Proof.

See Appendix C. ∎

To further characterize the $\sqrt{N}$ rate of convergence of $\rho_{u}^{N}$
and $\mu^{N}_{u}$, we obtain a central limit theorem for each (given in
Theorem 6 below). To state these, we first need to introduce some notation:
given two measurable spaces $(A,\mathcal{A})$ and $(B,\mathcal{B})$, a kernel
$M:A\times\mathcal{B}\to[0,\infty)$ from the first to the second, and a test
function $\varphi:B\rightarrow\mathbb{R}$, we use $M\varphi$ to denote the
real-valued function on $A$ defined by

$\displaystyle(M\varphi)(a):=\int M(a,db)\varphi(b)\quad\forall a\in A,$

assuming that the above integrals are well-defined.

As we will see in Theorem 6, the asymptotic variances of our estimators
decompose into sums of $\left|\mathbb{T}_{u}\right|$ terms, each one
accounting for the variance introduced by the computations carried out in a
different descendant $v$ of $u$. The term corresponding to a particular
descendant $v$ involves a kernel
$\Gamma_{v,u}:\bm{E}_{v}\times\bm{\mathcal{E}}_{u}\to[0,\infty)$ mapping
$\gamma_{v}$ to $\gamma_{u}$:

$\displaystyle\gamma_{v}(\Gamma_{v,u}\varphi)=\gamma_{u}(\varphi)\quad\forall\varphi\in\mathcal{B}_{b}(\bm{E}_{u}).$
(15)

The kernel is defined recursively:
$\Gamma_{u,u}(\bm{x}_{u},d\bm{y}_{u})=\delta_{\bm{x}_{u}}(d\bm{y}_{u})$,

$\Gamma_{v,u}(\bm{x}_{v},d\bm{y}_{u}):=\delta_{\bm{x}_{v}}(d\bm{y}_{v})\gamma_{\mathcal{C}_{u}}^{\not
v}(d\bm{y}_{\mathcal{C}_{u}}^{\not
v})w_{u_{-}}(\bm{y}_{\mathcal{C}_{u}})K_{u}(\bm{x}_{\mathcal{C}_{u}},dy_{u})\quad\forall
v\in\mathcal{C}_{u}$ (16)

where $\gamma_{\mathcal{C}_{u}}^{\not v}:=\prod_{r\in\mathcal{C}_{u}^{\not
v}}\gamma_{r}$ with $\mathcal{C}_{u}^{\not
v}:=\mathcal{C}_{u}\setminus\\{v\\}$. For all other descendants $v\neq u$, we
set

$\Gamma_{v,u}=\Gamma_{v,r_{1}}\Gamma_{r_{1},r_{2}}\dots\Gamma_{r_{l},u},$ (17)

where $v,r_{1},\dots,r_{l},u$ denotes the branch of $\mathbb{T}_{u}$
connecting $v$ and $u$. Because $\Gamma_{v,u}$ trivially satisfies (15) if $v$
is a child of $u$, it follows from (17) that (15) also holds for all other
descendants $v$ of $u$. Note that the collection
$(\Gamma_{u,v})_{u\in\mathbb{T},v\in\mathbb{T}_{u}}$ of all these kernels
amounts to a generalization of the semigroup describing the propagation of
local sampling errors in standard SMC that is central to much of its analysis
(e.g. see [17]).

######  Theorem 6 (Central limit theorems).

If Assumptions 1–2 hold, then, as $N\to\infty$,

$\displaystyle N^{1/2}\left(\rho^{N}_{u}(\varphi)-\rho_{u}(\varphi)\right)$
$\displaystyle\Rightarrow\mathcal{N}(0,\sigma^{2}_{\rho_{u}}(\varphi)),\quad
N^{1/2}\left(\mu^{N}_{u}(\varphi)-\mu_{u}(\varphi)\right)\Rightarrow\mathcal{N}(0,\sigma^{2}_{\mu_{u}}(\varphi)),$

for any given $u$ in $\mathbb{T}$ and $\varphi$ in
$\mathcal{B}_{b}(\bm{E}_{u})$, where $\Rightarrow$ denotes convergence in
distribution,

$\displaystyle\sigma^{2}_{\rho_{u}}(\varphi)$
$\displaystyle:=\sum_{v\in\mathbb{T}_{u}}\pi_{v}([\mathcal{Z}_{v}\Gamma_{v,u}[w_{u}\varphi]-\rho_{u}(\varphi)]^{2}),$
$\displaystyle\sigma^{2}_{\mu_{u}}(\varphi)$
$\displaystyle:=\sum_{v\in\mathbb{T}_{u}}\pi_{v}([\mathcal{Z}_{v}\Gamma_{v,u}[w_{u}Z_{u}^{-1}[\varphi-\mu_{u}(\varphi)]]]^{2}).$

In particular,
$N^{1/2}\left(Z_{u}^{N}-Z_{u}\right)\Rightarrow\mathcal{N}(0,\sigma^{2}_{Z_{u}})$
as $N\to\infty$ with

$\sigma^{2}_{Z_{u}}:=Z_{u}^{2}\sum_{v\in\mathbb{T}_{u}}\pi_{v}\left(\left[\frac{d\mu_{u}^{v}}{d\pi_{v}}-1\right]^{2}\right),$
(18)

where $\mu_{u}^{v}$ denotes the $\bm{E}_{v}$-marginal of $\mu_{u}$ (i.e.
$\mu_{u}^{v}(A):=\mu_{u}(A\times E_{\mathbb{T}_{u}\backslash\mathbb{T}_{v}})$
for all $A$ in $\bm{\mathcal{E}}_{v}$).

###### Proof.

See Appendix F. ∎

###  4.1 Methods of proof

The approximations to the targets are obtained from those to the extended
auxiliary measures via a single importance sampling step, cf. (11). Hence,
just as most results for ASMC are easily extracted from those for SMC [31],
the fundamental objects we need to analyze here are the approximations to the
extended auxiliary measures.

To this end, note that the main difference between DaC-SMC (Algorithm 3) and
standard ASMC (Algorithm 2) is the correction step: in the former case, we
employ the product-form estimator (8) while, in the latter, we instead use the
usual estimator (3). Consequently, our proofs for the above results combine
well-known methods previously used to establish analogous results for standard
SMC and ASMC with novel techniques that control the errors introduced by the
product-form estimators embedded within DaC-SMC.

At a conceptual level, we expand the approximation error into a sum of
products of ‘local errors’, each accounting for the approximations made at a
different node, and we control them separately (in particular, we show that,
just as in the standard case, these errors are $\mathcal{O}(N^{-1/2})$).
However, in contrast with standard SMC and ASMC whose local errors only get
‘propagated forward in time’, those of DaC-SMC get multiplied together as we
move up the tree. Consequently, we end up with a far greater number of terms
in our expansion than normal. Herein, lies the main novelty of our analysis:
dispatching these extra products-of-local-errors terms. Because the products
are taken over nodes on separate branches, the errors in the products are
independent of each other. Hence, the product of $k$ of these
$\mathcal{O}(N^{-1/2})$ errors should intuitively be $\mathcal{O}(N^{-k/2})$;
indeed, this was already shown in [37, Lemma 1] for the simple case of a tree
composed of a single parent with an arbitrary number of children.

More specifically, we begin by proving the $L^{p}$ inequalities (Theorem 5).
To do so, we follow the iterative approach taken in [14, 44]: we call on the
Marcinkiewicz-Zygmund inequality ([17, Lemma 7.3.3]) to obtain $L^{p}$
inequalities for the empirical distributions of the particles indexed by the
tree’s leaves and we show, step-by-step, that the algorithm preserves these
inequalities. The resampling and mutation steps then follow from arguments
analogous to those in [14, 44]. For the correction step, we need to do some
extra work and show that the product
$\gamma_{\mathcal{C}_{u}}^{N}=\prod_{v\in\mathcal{C}_{v}}\gamma_{v}^{N}$ also
maintains the inequalities. That is, we need to prove that the $L^{p}$ norm of
the product of the errors is $\mathcal{O}(N^{-1/2})$; something we do roughly
in Lemma 3 by using the boundedness of the test and weight functions to
uniformly bound all but one of the errors in the product, and exploiting the
$\mathcal{O}(N^{-1/2})$ size of the remaining error. To complete the argument
for the correction step, we then resume with the approach of [14, 44] and show
that the re-weighting in (8) preserves the inequalities. Next, we extract the
laws of large numbers (Theorem 1) from the $L^{p}$ inequalities using the
usual approach involving Chebyshev’s inequality and the Borel-Cantelli lemma
(e.g. [6, p. 17]). The weak convergence (Theorem 2) then follows using
standard techniques described in detail for SMC in [53, Section S1].

The proof of the unbiasedness in Theorem 3 is a straightforward generalization
of the arguments used in [11, Section 16.4.1] to establish the analogous
result for standard SMC. It too works recursively: we show that
$\gamma^{N}_{u_{-}}$ is an unbiased estimator for $\gamma_{u_{-}}$ if
$\gamma^{N}_{v_{-}}$ is a unbiased estimator for $\gamma_{v_{-}}$ for each
non-leaf child $v$ of $u$. The unbiasedness of $\rho_{u}^{N}$ then follows
with a simple use of the tower property.

To obtain the bias bounds in Theorem 4, we follow the approach of [48]. This
argument also works its way through the algorithm step by step, this time
showing that each step increases the bias by at most $\mathcal{O}(N^{-1})$.
Just as with the $L^{p}$ inequalities, our innovation here (Lemma 8) deals
with the product
$\gamma_{\mathcal{C}_{u}}^{N}=\prod_{v\in\mathcal{C}_{v}}\gamma_{v}^{N}$. As
in that case, we obtain the necessary bias bound for the product by uniformly
bounding the bias of all but one of the approximations in the product and
using a previously derived $\mathcal{O}(N^{-1})$ bound for the remaining
approximation.

Our proof for the central limit theorems (Theorem 6) follows the conceptual
approach outlined above and extends the arguments in [17, Chapter 9].
Similarly to in that chapter, we obtain an expression for the approximation
error in terms of local errors (and their ‘propagations up the tree’), the
difference being that the expression is not just a sum of the local errors,
but of their products too. An analogous argument to [17, Corollary 9.3.1] then
shows that the local errors are asymptotically normal and independent. Hence,
if we are able to demonstrate that their products are $o(N^{-1/2})$, then the
CLTs follow as in [17, p. 301] by applying Slutsky’s theorem and the
continuous mapping theorem, and exploiting the fact that a linear combination
of independent normal random variables is normal. To this end, we prove the
following product version of the $L^{2}$ inequality in Theorem 5:

######  Theorem 7 (Product $L^{2}$ inequality).

If, in addition to Assumptions 1–2, $u$ and $v$ lie in separate branches (i.e.
$u\not\in\mathbb{T}_{v}$ and $v\not\in\mathbb{T}_{u}$), then there exists a
constant $C_{u,v}<\infty$ such that

$\mathbb{E}\left[(\gamma^{N}_{u}-\gamma_{u})\times(\gamma_{v}^{N}-\gamma_{v})(\varphi)^{2}\right]^{\frac{1}{2}}\leq\frac{C_{u,v}\left|\left|{\varphi}\right|\right|}{N}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}\times\bm{E}_{v}).$

###### Proof.

See Appendix G.∎

To argue the above, we expand [37, Lemma 1]’s proof using an inductive
approach similar to those we take for the $L^{p}$ inequalities, unbiasedness,
and bias bounds: one by one, we show that each step of algorithm preserves the
above inequality. For the base case dealing with empirical distributions of
the particles indexed by the tree’s leaves, we call on [37, Lemma 1].

###  4.2 Optimal intermediate targets, auxiliary measures, and proposal
kernels

While the optimal choice of intermediate targets, auxiliary measures, and
proposal kernels generally depends on the particular average we are interested
in estimating, there is one choice that leads to zero variance estimates of
the final target’s normalizing constant $Z_{\mathfrak{r}}$. In particular,
suppose that the final target $\rho_{\mathfrak{r}}$ and the underlying space
$(\bm{E}_{\mathfrak{r}},\bm{\mathcal{E}}_{\mathfrak{r}})$ are nice enough
(e.g. $(\bm{E}_{\mathfrak{r}},\bm{\mathcal{E}}_{\mathfrak{r}})$ is Borel [32,
Theorem 8.5], see also [24]) that, for each $u$ in
$\mathbb{T}^{\not\partial}$, there exists a regular conditional probability
distribution mapping the $\bm{E}_{\mathcal{C}_{u}}$-marginal
$\mu_{\mathfrak{r}}^{u_{-}}$ (i.e.
$\mu_{\mathfrak{r}}^{u_{-}}(A):=\mu_{\mathfrak{r}}(A\times\bm{E}_{\mathbb{T}\backslash(\mathbb{T}_{u}^{\not
u})})$ for all $A$ in $\bm{\mathcal{E}}_{\mathcal{C}_{u}}$) to its
$\bm{E}_{u}$-marginal $\mu_{\mathfrak{r}}^{u}$. That is, a Markov kernel
$M_{u}:\bm{E}_{\mathcal{C}_{u}}\times\mathcal{E}_{u}\to[0,1]$ satisfying
$\mu_{\mathfrak{r}}^{u_{-}}\otimes M_{u}=\mu_{\mathfrak{r}}^{u}$. In this
case, setting

$\mu_{u}:=\mu_{\mathfrak{r}}^{u},\quad K_{u}:=M_{u}\enskip\forall
u\in\mathbb{T},\quad\pi_{u_{-}}:=\mu_{\mathfrak{r}}^{u_{-}}\enskip\forall
u\in\mathbb{T}^{\not\partial},$

we have that

$\pi_{v}=\pi_{v_{-}}\otimes M_{v}=\mu_{\mathfrak{r}}^{v_{-}}\otimes
M_{v}=\mu_{\mathfrak{r}}^{v}=\mu_{u}^{v}\quad\forall
v\in\mathbb{T}_{u},\enskip u\in\mathbb{T}.$

It follows from (18) that our estimator $Z_{\mathfrak{r}}^{N}$ for
$Z_{\mathfrak{r}}$ achieves zero (asymptotic) variance (and similarly for
$Z_{u}^{N}$ for all other nodes $u$). The above choices generalize those for
which ASMC yields zero-variance estimates of the marginal likelihood in the
context of filtering, cf. [28, Proposition 2]. Note also that, just as in the
ASMC case, the normalizing constants in the intermediate targets auxiliary
measures are immaterial and we are free to choose them as we find convenient:
something unsurprising given that they do not influence the resampling
operation, and ultimately cancel in the computation of the normalizing
constant estimates. However, perhaps a bit unexpectedly, and in contrast with
the ASMC case, these choices do not necessarily lead to zero _finite-sample_
variance estimates for the normalizing constants. In particular, given the
above,

$\displaystyle
w_{u}=\frac{d\rho_{u}}{d\gamma_{u}}=\frac{Z_{u}}{\mathcal{Z}_{u}}\frac{d\mu_{u}}{d\pi_{u}}=\frac{Z_{u}}{\mathcal{Z}_{u}},\quad
w_{u_{-}}=\frac{d\gamma_{u_{-}}}{d\gamma_{\mathcal{C}_{u}}}=\frac{\mathcal{Z}_{u}}{\mathcal{Z}_{\mathcal{C}_{u}}}\frac{d\pi_{u_{-}}}{d\pi_{\mathcal{C}_{u}}}=\frac{\mathcal{Z}_{u}}{\mathcal{Z}_{\mathcal{C}_{u}}}\frac{d\mu^{u_{-}}_{\mathfrak{r}}}{d\prod_{v\in\mathcal{C}_{u}}\mu^{v}_{\mathfrak{r}}},$
$\displaystyle\Rightarrow
Z^{N}_{u}=\gamma_{u}^{N}(w_{u})=\frac{Z_{u}}{\mathcal{Z}_{u}}\mathcal{Z}_{u}^{N}=\frac{Z_{u}}{\mathcal{Z}_{u}}\gamma_{\mathcal{C}_{u}}^{N}(w_{u_{-}})=\frac{Z_{u}}{\mathcal{Z}_{\mathcal{C}_{u}}}\gamma_{\mathcal{C}_{u}}^{N}\left(\frac{d\mu^{u_{-}}_{\mathfrak{r}}}{d\prod_{v\in\mathcal{C}_{u}}\mu^{v}_{\mathfrak{r}}}\right),$

and the rightmost term will not be constant unless
$\mu^{u_{-}}_{\mathfrak{r}}=\prod_{v\in\mathcal{C}_{u}}\mu^{v}_{\mathfrak{r}}$.
On the other hand, were this to be the case for all nodes $u$, we could
iterate the above down the tree to find that $Z^{N}_{u}=Z_{u}$ for all $u$.

In short, the above choices lead to zero-finite-sample-variance estimates for
the normalizing constants if and only if $\mathbb{T}$ precisely matches the
target’s dependence structure. Some insight as to why this might be the case
can be gleaned by comparing the asymptotic and finite variance expressions
given in [37, Theorem 1] for product-form estimators. In essence, these
choices ensure that the ‘one-dimensional’ marginals of the sampling
distribution $\pi_{u}$ at node $u$ coincide with those of $\mu_{u}$ and negate
the corresponding $\mathcal{O}(N^{-1/2})$ terms in the finite-sample variance
expansions (i.e. those also featuring in the asymptotic variance expansions).
However, the sampling distribution factorizes over $u$’s children and, unless
the target does so too, the higher-dimensional marginals are mismatched and
the corresponding $o(N^{-1/2})$ terms in the expansion persevere.

Just as is the case in standard ASMC [28], the above choices can rarely be
implemented in practice: with few exceptions, $\mu_{\mathfrak{r}}$’s marginals
cannot be computed explicitly and the conditional distributions are unknown
(much less can be sampled from). A more common situation in practice is that
where the intermediate targets are fixed and we seek to choose only the
proposal kernels and/or auxiliary measures. In standard ASMC (cf. [11, Chapter
10] and references therein), this is typically done by deriving the ‘locally
optimal’ proposal kernel and auxiliary measure for each step of the algorithm
(i.e. those that minimize the finite-sample variance of the normalizing
constant estimate obtained at that step) and, if these prove intractable,
approximating them as best as one can. To extend this approach to DaC-SMC, we
provide the following generalization of the characterization for the locally
optimal kernels and measures previously obtained for SMC and ASMC (e.g. see
[11, Theorems 10.1–10.2]).

######  Theorem 8.

Fix any $u$ in $\mathbb{T}^{\not\partial}$ and $N>0$, and let $\mathcal{S}$
denote the set of all $(\gamma_{u_{-}},K_{u})$ pairs satisfying Assumption 1
and $\mathcal{S}(K_{u})$ and $\mathcal{S}(\gamma_{u_{-}})$ be its slices:

$\displaystyle\mathcal{S}:=\\{(\gamma_{u_{-}},K_{u}):\gamma_{u_{-}}\sim\gamma_{\mathcal{C}_{u}},\enskip\rho_{u}\sim\gamma_{u_{-}}\otimes
K_{u}\\},$
$\displaystyle\mathcal{S}(K_{u}):=\\{\gamma:(\gamma,K_{u})\in\mathcal{S}\\},\quad\mathcal{S}(\gamma_{u_{-}}):=\\{K:(\gamma_{u_{-}},K)\in\mathcal{S}\\},\quad\forall(\gamma_{u_{-}},K_{u})\in\mathcal{S}.$

If $f(\gamma_{u_{-}},K_{u})$ denotes $\textrm{Var}(Z_{u}^{N})$ written
explicitly as a function of $(\gamma_{u_{-}},K_{u})$, then, for all
$(\gamma_{u_{-}},K_{u})$ in $\mathcal{S}$ and $C>0$,

$\displaystyle C\sqrt{K_{u}\omega_{u}(K_{u})^{2}}\rho_{\mathcal{C}_{u}}$
$\displaystyle\in\operatorname*{arginf}_{\gamma\in\mathcal{S}(K_{u})}f(\gamma,K_{u}),$
(19) $\displaystyle M_{u}$
$\displaystyle\in\operatorname*{arginf}_{K\in\mathcal{S}(\gamma_{u_{-}})}f(\gamma_{u_{-}},K),$
(20) $\displaystyle(C\mu_{u}^{u_{-}},M_{u})$
$\displaystyle\in\operatorname*{arginf}_{(\gamma,K)\in\mathcal{S}}f(\gamma,K),$
(21)

where $\omega_{u}(K_{u}):=d\rho_{u}/d\rho_{\mathcal{C}_{u}}\otimes K_{u}$ and
$M_{u}$ denotes the regular conditional probability distribution mapping
$\mu_{u}^{u_{-}}$ to $\mu_{u}$ (i.e. $\mu_{u}^{u_{-}}\otimes M_{u}=\mu_{u}$),
assuming that it exists.

###### Proof.

See Appendix H. ∎

###  4.3 Comparison with ASMC

Theorems 1–6 demonstrate that DaC-SMC is a well-founded algorithm: it leads to
estimators for the unnormalized and normalized targets with the same basic
properties as those possessed by SMC and ASMC estimators [17]. It is only
natural to now ask whether it performs better than its ASMC analogue
(Algorithm 2).

To answer this question, suppose we are interested in approximating the final
target $\rho_{\mathfrak{r}}$ and its normalization $\mu_{\mathfrak{r}}$. We
first need to figure out how to approximate $\rho_{\mathfrak{r}}$ and
$\mu_{\mathfrak{r}}$ using ASMC. This requires somehow lumping the nodes in
$\mathbb{T}$ to obtain a line $\mathbb{L}$ (Remark 1). While there are many
ways that it could be done, we focus on perhaps the most obvious: merge
together all nodes in each level of the tree to obtain a line
$\mathcal{L}_{0}\to\mathcal{L}_{1}\to\dots\to\mathcal{L}_{T}=\\{\mathfrak{r}\\}$
whose length equals $\mathbb{T}$’s depth. Taking the corresponding products of
the auxiliary measures,

$\gamma_{t_{-}}:=\prod_{u\in\mathcal{L}_{t}^{\not\partial}}\gamma_{u_{-}},\quad\pi_{t_{-}}:=\prod_{u\in\mathcal{L}_{t}^{\not\partial}}\pi_{u_{-}},\quad\forall
t=1,\dots,T,$

we obtain sequences $(\gamma_{t_{-}})_{t=1}^{T}$ and $(\pi_{t_{-}})_{t=1}^{T}$
that can be used in a standard ASMC set-up. In particular, replacing
$(\gamma_{u_{-}})_{u\in\mathbb{T}^{\not\partial}}$,
$(\pi_{u_{-}})_{u\in\mathbb{T}^{\not\partial}}$, and
$(K_{u})_{u\in\mathbb{T}}$ with $(\gamma_{t_{-}})_{t=1}^{T}$,
$(\pi_{t_{-}})_{t=1}^{T}$, and $(K_{t})_{t=0}^{T}$, for some proposal kernels
$K_{t}:\bm{E}_{t_{-}}\times\mathcal{E}_{t}\to[0,1]$ mapping from
$\bm{E}_{t_{-}}:=\prod_{u\in\mathcal{L}_{t}^{\not\partial}}\bm{E}_{\mathcal{C}_{u}}$
to $\mathcal{E}_{t}:=\mathcal{E}_{\mathcal{L}_{t}}$, Algorithm 3 reduces to
Algorithm 2. To keep both algorithms on as equal a footing as possible, we set
the ASMC proposal kernels to be the corresponding products of the DaC-SMC
ones:

$K_{t}(\bm{x}_{t_{-}},dx_{t}):=\left(\prod_{u\in\mathcal{L}_{t}^{\partial}}K_{u}(dx_{u})\right)\left(\prod_{u\in\mathcal{L}_{t}^{\not\partial}}K_{u}(\bm{x}_{\mathcal{C}_{u}},dx_{u})\right)\quad\forall
t=0,\dots,T.$ (22)

Theoretically, Theorems 1–6 tell us that both ASMC and DaC-SMC produce
consistent estimators of $\rho_{\mathfrak{r}}$ and $\mu_{\mathfrak{r}}$,
unbiased in the case of $\rho_{\mathfrak{r}}$ and with an
$\mathcal{O}(N^{-1})$ bias in that of $\mu_{\mathfrak{r}}$, whose errors are
$\mathcal{O}(N^{-1/2})$ and asymptotically normal. Comparing the asymptotic
variance expressions in Theorem 6 for both algorithms, we find that the DaC-
SMC estimators are more statistically efficient than the ASMC ones:

######  Theorem 9.

If Assumptions 1–2 are satisfied, $\varphi$ belongs to
$\mathcal{B}_{b}(\bm{E}_{\mathfrak{r}})$, and
$\sigma^{2}_{\rho_{\mathfrak{r}},smc}(\varphi)$ and
$\sigma^{2}_{\mu_{\mathfrak{r}},smc}(\varphi)$ denote the respective
asymptotic variances of the ASMC estimators for $\rho_{\mathfrak{r}}(\varphi)$
and $\mu_{\mathfrak{r}}(\varphi)$, then

$\sigma^{2}_{\rho_{\mathfrak{r}}}(\varphi)\leq\sigma^{2}_{\rho_{\mathfrak{r}},smc}(\varphi),\qquad\sigma^{2}_{\mu_{\mathfrak{r}}}(\varphi)\leq\sigma^{2}_{\mu_{\mathfrak{r}},smc}(\varphi).$

In particular, the asymptotic variance of the DaC-SMC normalizing constant
estimate is bounded above by that of the SMC estimate.

###### Proof.

See Appendix I. ∎

Of course, the above theorem, and this entire comparison, applies only if we
are using factorized proposal kernels in (22) for ASMC, a needless restriction
for this algorithm and one that will often be detrimental to its performance.
Indeed, unless the intermediate targets exhibit the same kind of conditional
independence encoded in the factorized kernels, standard results for ASMC
(e.g. [11, Theorems 10.1,10.2] or [28, Proposition 2]) tell us that these
kernels will not be optimal in any usual sense.

However, for such factorized proposals, the gains in statistical efficiency
can sometimes be drastic (e.g. exponential in the tree’s degree), see Examples
2 and 4 and Section 3.4 in [37]. In exchange, DaC-SMC generally has a higher
computational cost. More specifically, the mutation steps of DaC-SMC require
drawing $N$ samples from $K_{u}$ for each $u$ in $\mathbb{T}$, while ASMC
requires drawing $N$ samples from $K_{t}$ for each $t=0,\dots,T$. However,
$K_{t}$’s definition implies that drawing a sample from it is equivalent to
drawing one from each $K_{u}$ with $u$ in $\mathcal{L}_{t}$ and concatenating
them. Hence, the cost incurred by the mutation steps of both algorithms is the
same. The cost of the correction and resampling steps, however, is not the
same. In the case of ASMC, we need to sample $N$ particles from
$\pi_{t_{-}}^{N}$ for each $t$ in $0,\dots,T$ which requires only $N$
evaluations of $w_{t_{-}}=\prod_{u\in\mathcal{L}_{t}^{\not\partial}}w_{u_{-}}$
per $t$. For DaC-SMC, we must instead sample $N$ particles from
$\pi_{u_{-}}^{N}$ for each $u$ in $\mathbb{T}^{\not\partial}$. In the absence
of any sort of special structure in the auxiliary measures, this involves
$N^{c_{u}}$ evaluations of $w_{u_{-}}$ per $u$ in $\mathbb{T}^{\not\partial}$.
Hence, for general auxiliary measures, DaC-SMC’s cost is $\mathcal{O}(N^{d})$,
where $d$ denotes $\mathbb{T}$’s degree, while ASMC’s is just
$\mathcal{O}(N)$.

However, this extra cost materializes only in nodes with more than one child
and concentrates in those with $d$ children. If these nodes only feature high
up in the tree, then the high cost incurred by the generation of the particles
indexed by the nodes’ children can outweigh the extra overhead incurred by the
resampling and correction steps. Similarly if the proposal kernels for those
nodes are particularly expensive to sample from. It is worth noting that in
many settings of practical interest (e.g. [13, 20, 21]), $d=2$ and this cost
can be borne. Otherwise, as described in Section 3.3, DaC-SMCs cost may be
brought down by choosing auxiliary measures that are partially or fully
factorized (or sums of such measures)—and when these differ too much from the
optimal auxiliary measures the impact of this discrepancy upon the variance
can be mitigated via tempering [19]. Lastly, in cases where $d$ is not small
and the preferred auxiliary measures do not factorize sufficiently, one can
instead turn to the ‘incomplete’ estimators (13) similar to those successfully
employed in the U-statistic literature to tackle analogous issues (cf. [35]
and references therein).

##  5 Discussion

In this paper, we study DaC-SMC algorithm introduced in [41] and show that it
is theoretically well-founded: it possesses the same basic properties that
standard SMC algorithms do (Section 4). To achieve this, we combine well-known
methods previously used to study standard SMC algorithms with novel techniques
that control the errors introduced by the product-form estimators (8) embedded
within DaC-SMC (see Section 4.1 for an overview and Appendices C–G for the
details). Our analysis here can be sharpened and refined in various ways. The
positivity requirements in Assumptions 1 can be circumvented by emulating the
use of stopping times in [17]. The boundedness ones in Assumptions 2 can be
avoided by introducing families of appropriately integrable functions at each
node, e.g. similarly as in [10, 22]. While the additional requirement in
Theorem 4 that the weight functions are lower bounded is common in the SMC
literature (e.g. [18, Proposition 9.5.6] and [48]), we anticipate that it
could be relaxed by considering a more restricted class of test functions
(than that of all bounded and measurable ones) and dealing with the possible
‘extinction’ of the particle cloud as done in [17, Theorem 7.4.3]. Lastly, we
do not believe that the particular lumping construction in Section 4.3 is
crucial for the variance bounds in Theorem 9 to hold, only that the ASMC
auxiliary measures and proposals are obtained by taking appropriate products
of the DaC-SMC ones. In fact, in Appendix H we actually show that lumping the
children of a node never reduces the variance of DaC-SMC estimators. However,
to argue the bounds in full generality one would likely have to further show
that lumping ‘cousins’ does not lower the variance, and we opted to focus on
that particular lumping construct so not overly complicate Section 4.3 and
Theorem 8’s proof.

Whenever the auxiliary measures and proposal kernels used for DaC-SMC are the
corresponding products of those used for standard ASMC, the former proves more
statistically efficient than the latter (Section 4.3): using the same number
of particles, DaC-SMC estimators achieve lower variances than their standard
counterparts. In exchange, DaC-SMC generally carries a higher computational
cost. How much higher depends on the amount of product structure present in
the auxiliary measures $(\gamma_{u_{-}})_{u\in\mathbb{T}^{\not\partial}}$
employed (Section 3.3). If they are fully-factorized (or sums of such
measures), then the DaC-SMC’s cost is $\mathcal{O}(N)$ where $N$ denotes the
number of particles (the same as ASMC). If they do not factorize at all (nor
break down into sums of partially-factorized functions), then DaC-SMC’s cost
is $\mathcal{O}(N^{d})$ where $d$ denotes the tree’s degree. For anything
between these two extremes, the cost is $\mathcal{O}(N^{k})$ where $1<k<d$
depends on amount of factorization (the greater it is, the lower $k$ is).

At each node $u$, the algorithm approximates the auxiliary measure
$\gamma_{u_{-}}$ with the product-form $\gamma^{N}_{\mathcal{C}_{u}}$ and then
corrects this approximation using the weight function $w_{u_{-}}$ (8). The
corrected approximation, $\gamma_{u_{-}}^{N}$, is then extended using a
mutation step, and the resulting extended approximation, $\gamma_{u}^{N}$, is
a corrected using the weight function $w_{u}$ (11) to produce an
approximation, $\rho_{u}^{N}$, of the intermediate target $\rho_{u}$. This
means that any non-product structure in $\rho_{u}$ must be introduced into our
approximations via these two re-weighting steps (aside from structure linking
the components indexed by $u$’s children with $u$ itself, which can be
introduced using the mutation kernels). Consequently, if the intermediate
target $\rho_{u}$ posses pronounced non-product structure, we end up with
competing interests: (a) we would like to choose the auxiliary measure
$\gamma_{u_{-}}$ as close to factorized as possible so that the algorithm’s
cost is low (Section 3.3) and the first re-weighting does not necessitate
large corrections and lead to high-variance approximations of
$\gamma_{u_{-}}$; and (b) we would like to choose the auxiliary measure
$\gamma_{u_{-}}$ close $\rho_{u}$. In fact, the analysis in Section 4.2
suggests that, purely in terms of variance and disregarding any considerations
of the algorithm’s cost, $\gamma_{u_{-}}$ should be chosen as close as
possible to the corresponding marginal of $\rho_{u}$. How to balance these two
competing considerations is an interesting question beyond the scope of this
paper and likely best dealt with on a case-by-case basis.

Alternatively, it might also be possible to introduce non-product structure
without blowing up the computational cost or the variance in the correction
steps using tempering [19], ‘incomplete’ versions of product-form estimators
(13), or generalizations of these estimators that account for more complicated
dependence structures (see Appendix A.2 for an example). We find each of these
variants, and the question of which one should be used under what
circumstances, worthy of further investigation; and we believe this can be
done using extensions of the techniques we employed in Appendices C–G.
Similarly for the other, perhaps more familiar, refinements that will also
likely improve the algorithm’s performance in practice (e.g. adaptive
resampling, more efficient resampling schemes, and varying the computational
power spent at each node; see Section 3.4).

Lastly, we anticipate that a variety of interesting applications of the DaC-
SMC algorithm will involve its merger with other well-known Monte Carlo
methodology and our results pave the way for such combinations. For instance,
its merger with pseudo-marginal MCMC [3] to tackle targets with intractable
densities and that with particle-marginal Metropolis-Hastings [2] to allow for
Metropolis-Hastings proposals with intractable densities; two approaches
justified by the unbiasedness in Theorem 3. It can also be used within general
particle MCMC [2], including particle Gibbs (e.g. pertinent for models such as
Markov random fields with collections of correlated parameter variables
lacking time-series structure), the formal justification of which requires the
characterization of the joint law of all random variables generated during the
running of the Algorithm 3. This characterization can be found in [15, Chapter
8] which generalizes the argument given in [41, Appendix A.1] for balanced
binary trees.

In summary, DaC-SMC is a theoretically-sound algorithm and a promising
extension of SMC. Its development, however, still lies in its infancy and
there remain many open questions regarding its use. For instance, ‘how exactly
should we carry out the resampling and which nodes should we allocate more
computational power to?’, ‘how do we introduce non-product structure into the
approximations while balancing the algorithm’s cost and the variance of the
correction steps?’, and ‘with which other Monte algorithms would it prove
fruitful to merge DaC-SMC and for what type of targets should which merger be
used for?’. We look forward to their resolution.

## References

  * [1] L. Aitchison.  Tensor Monte Carlo: Particle methods for the GPU era.  In Adv Neural Inf Process Syst, volume 32, pages 7148–7157, 2019\.
  * [2] C. Andrieu, A. Doucet, and R. Holenstein.  Particle Markov chain Monte Carlo methods.  J R Stat Soc B, 72(3):269–342, 2010.  doi:10.1111/j.1467-9868.2009.00736.x.
  * [3] C. Andrieu and G. O. Roberts.  The pseudo-marginal approach for efficient Monte Carlo computations.  Ann Stat, 37(2):697–725, 2009.  doi:10.1214/07-AOS574.
  * [4] S. Asmussen and W. Glynn.  Stochastic Simulation: Algorithms and Analysis.  Springer-Verlag New York, 2007.  doi:10.1007/978-0-387-69033-9.
  * [5] P. Berti, L. Pratelli, and P. Rigo.  Almost sure weak convergence of random probability measures.  Stochastics and Stochastics Reports, 78(2):91–97, 2006.  doi:10.1080/17442500600745359.
  * [6] A. Boustati, O. D. Akyildiz, T. Damoulas, and A. M. Johansen.  Generalised Bayesian Filtering via Sequential Monte Carlo.  Adv Neural Inf Process Syst, 33:418–429, 2020.
  * [7] J. Carpenter, P. Clifford, and P. Fearnhead.  An improved particle filter for non-linear problems.  IEE P-Radar Son Nav, 146(1):2–7, 1999.
  * [8] R. Chan, M. Pollock, A. M. Johansen, and G. O. Roberts.  Divide-and-conquer Monte Carlo fusion.  arXiv preprint arXiv:2110.07265, 2021.
  * [9] N. Chopin.  A sequential particle filter method for static models.  Biometrika, 89(3):539–551, 2002.  doi:10.1093/biomet/89.3.539.
  * [10] N. Chopin.  Central limit theorem for sequential Monte Carlo methods and its application to Bayesian inference.  Ann Stat, 32(6):2385–2411, 2004.  doi:10.1214/009053604000000698.
  * [11] N. Chopin and O. Papaspiliopoulos.  An Introduction to Sequential Monte Carlo.  Springer, 2020.  doi:10.1007/978-3-030-47845-2.
  * [12] M. E. Conway.  A multiprocessor system design.  In Proceedings of the November 12-14, 1963, Fall Joint Computer Conference, pages 139–146, 1963.  doi:10.1145/1463822.1463838.
  * [13] A. Corneflos, N. Chopin, and S. Särkkä.  A parallel-in-time particle smoother.  In preparation, 2021.
  * [14] D. Crisan and A. Doucet.  A survey of convergence results on particle filtering methods for practitioners.  IEEE T Signal Proces, 50(3):736–746, 2002.  doi:10.1109/78.984773.
  * [15] F. R. Crucinio.  Some Interacting Particle Methods with Non-Standard Interactions.  Ph.D. thesis, University of Warwick, 2021.
  * [16] H. Dai, M. Pollock, and G. Roberts.  Monte Carlo fusion.  J Appl Probab, 56(1):174–191, 2019.  doi:10.1017/jpr.2019.12.
  * [17] P. Del Moral.  Feynman-Kac Formulae.  Springer, 2004.  doi:10.1007/978-1-4684-9393-1.
  * [18] P. Del Moral.  Mean field simulation for Monte Carlo integration.  Monographs on Statistics and Applied Probability. CRC Press, 1 edition, 2013.  doi:10.1201/b14924.
  * [19] P. Del Moral, A. Doucet, and A. Jasra.  Sequential monte carlo samplers.  J R Stat Soc B, 68(3):411–436, 2006.  doi:10.1111/j.1467-9868.2006.00553.x.
  * [20] D. Ding.  Monte Carlo Algorithms for Hypothesis Testing and for Hidden Markov Models.  Ph.D. thesis, Imperial College, 2019.
  * [21] D. Ding and A. Gandy.  Tree-based particle smoothing algorithms in a hidden markov model.  arXiv preprint arXiv:1808.08400, 2018.
  * [22] R. Douc and É. Moulines.  Limit theorems for weighted samples with applications to sequential monte carlo methods.  Ann Stat, 36(5):2344–2376, 2008.  doi:10.1214/07-AOS514.
  * [23] A. Doucet and A. M. Johansen.  A tutorial on particle filtering and smoothing: Fiteen years later.  In D. Crisan and B. Rozovsky, editors, The Oxford Handbook of Nonlinear Filtering, pages 656–704. Oxford University Press, 2011.
  * [24] A. M. Faden.  The existence of regular conditional probabilities: Necessary and sufficient conditions.  Ann Probab, 13(1):288–298, 1985.  doi:10.1214/aop/1176993081.
  * [25] P. Fearnhead and H. R. Künsch.  Particle filters and data assimilation.  Annu Rev Stat Appl, 5(1):421–449, 2018.  doi:10.1146/annurev-statistics-031017-100232.
  * [26] A. Gelman and J. Hill.  Data Analysis Using Regression and Multilevel/Hierarchical Models.  Cambridge University Press, 2006.  doi:10.1017/CBO9780511790942.
  * [27] M. Gerber, N. Chopin, and N. Whiteley.  Negative association, ordering and convergence of resampling methods.  Ann Stat, 47(4):2236–2260, 2019.  doi:10.1214/18-AOS1746.
  * [28] P. Guarniero, A. M. Johansen, and A. Lee.  The iterated auxiliary particle filter.  J Am Stat Assoc, 112(520):1636–1647, 2017.  doi:10.1080/01621459.2016.1222291.
  * [29] M. Gürbüzbalaban, X. Gao, Y. Hu, and L. Zhu.  Decentralized stochastic gradient langevin dynamics and hamiltonian monte carlo.  J Mach Learn Res, 22(239):1–69, 2021.
  * [30] S. W. Jewell.  Divide and conquer sequential Monte Carlo for phylogenetics.  Master thesis, University of British Columbia, 2015.
  * [31] A. M. Johansen and A. Doucet.  A note on the auxiliary particle filter.  Stat Probab Lett, 78(12):1498–1504, 2008.  doi:10.1016/j.spl.2008.01.032.
  * [32] O. Kallenberg.  Foundations of Modern Probability.  Springer, Cham, 2021.  doi:10.1007/978-3-030-61871-1.
  * [33] N. Kantas, A. Doucet, S. S. Singh, J. Maciejowski, and N. Chopin.  On Particle Methods for Parameter Estimation in State-Space Models.  Stat Sci, 30(3):328–351, 2015.  doi:10.1214/14-STS511.
  * [34] A. Kong, J.S. Liu, and W.H. Wong.  Sequential imputations and bayesian missing data problems.  J Am Stat Assoc, 89(425):278–288, 1994.  doi:10.1080/01621459.1994.10476469.
  * [35] X. Kong and W. Zheng.  Design based incomplete u-statistics.  Stat Sinica, 31:1593–1618, 2021.  doi:10.5705/ss.202019.0098.
  * [36] V. S. Korolyuk and Y. V. Borovskich.  Theory of U-statistics.  Springer Science & Business Media, 1994.
  * [37] J. Kuntz, F. R. Crucinio, and A. M. Johansen.  Product-form estimators: exploiting independence to scale up Monte Carlo.  arXiv preprint arXiv:2102.11575, 2021.
  * [38] A. Lee and N. Whiteley.  Variance estimation in the particle filter.  Biometrika, 105(1):609–625, 2018.  doi:10.1093/biomet/asy028.
  * [39] A. J. Lee.  U-Statistics: Theory and Practice.  CRC Press, 1990.
  * [40] M. T. Lin, J. L. Zhang, Q. Cheng, and R. Chen.  Independent particle filters.  J Am Stat Assoc, 100(472):1412–1421, 2005.  doi:10.1198/016214505000000349.
  * [41] F. Lindsten, A. M. Johansen, C. A. Naesseth, B. Kirkpatrick, T. B. Schön, J. A. D. Aston, and A. Bouchard-Côté.  Divide-and-Conquer with sequential Monte Carlo.  J Comput Graph Stat, 26(2):445–458, 2017.  doi:10.1080/10618600.2016.1237363.
  * [42] J. S. Liu.  Monte Carlo Strategies in Scientific Computing.  Springer-Verlag New York, 2001.  doi:10.1007/978-0-387-76371-2.
  * [43] T. G. Mattson, B. Sanders, and B. Massingill.  Patterns for parallel programming.  Pearson Education, 2004.
  * [44] J. Míguez, D. Crisan, and P. M. Djurić.  On the convergence of two sequential Monte Carlo methods for maximum a posteriori sequence estimation and stochastic global optimization.  Stat Comput, 23(1):91–107, 2013.  doi:10.1007/s11222-011-9294-4.
  * [45] S. Minsker, S. Srivastava, L. Lin, and D. B. Dunson.  Scalable and robust Bayesian inference via the median posterior.  In 31st Int Conf Mach Learn, volume 32, pages 1656–1664, 2014.
  * [46] R. M. Neal.  Annealed importance sampling.  Stat Comput, 11(2):125–139, 2001.  doi:10.1023/A:1008923215028.
  * [47] W. Neiswanger, C. Wang, and E. Xing.  Asymptotically Exact, Embarrassingly Parallel MCMC.  In 13th Conf Uncertain Artif Intell, pages 623–632, 2014.
  * [48] J. Olsson and T. Rydén.  The bootstrap particle filtering bias.  Lund University, Technical Report 929081, (2004:24), 2004.
  * [49] B. Paige and F. Wood.  Inference networks for sequential monte carlo in graphical models.  In 33rd Int Conf Mach Learn, volume 48, pages 3040–3049, 2016.
  * [50] A. Parayil, H. Bai, J. George, and P. Gurram.  Decentralized langevin dynamics for bayesian learning.  Adv Neural Inf Process Syst, 33, 2020.
  * [51] M. K. Pitt and N. Shephard.  Filtering via simulation: Auxiliary particle filters.  J Am Stat Assoc, 94(446):590–599, 1999.
  * [52] L. J. Rendell, A. M. Johansen, A. Lee, and N. Whiteley.  Global consensus Monte Carlo.  J Comput Graph Stat, 30(2):249–259, 2021.  doi:10.1080/10618600.2020.1811105.
  * [53] S. M. Schmon, G. Deligiannidis, A. Doucet, and M. K. Pitt.  Large-sample asymptotics of the pseudo-marginal method.  Biometrika, 108(1):37–51, 2020.  doi:10.1093/biomet/asaa044.
  * [54] S. L. Scott, A. W. Blocker, F. V. Bonassi, H. A. Chipman, E. I. George, and R. E. McCulloch.  Bayes and Big Data: The Consensus Monte Carlo Algorithm.  Int J Manag Sci Eng Manag, 11(2):78–88, 2016.  doi:10.1080/17509653.2016.1142191.
  * [55] M.-H. Tran, M. Scharth, M. K. Pitt, and R. Kohn.  Importance sampling squared for Bayesian inference in latent variable models.  arXiv preprint arXiv:1309.3339, 2013.
  * [56] M. Vono, N. Dobigeon, and P. Chainais.  Asymptotically exact data augmentation: Models, properties, and algorithms.  J Comput Graph Stat, 30(2):335–348, 2021.  doi:10.1080/10618600.2020.1826954.
  * [57] A. J. Walker.  An efficient method for generating discrete random variables with general distributions.  ACM T Math Software, 3(3):253–256, 1977.  doi:10.1145/355744.355749.
  * [58] X. Wang and D. B. Dunson.  Parallelizing MCMC via Weierstrass Sampler.  arXiv preprint arXiv:1312.4605, 2013.
  * [59] D.-M. Yuan and S.-J. Li.  Extensions of several classical results for independent and identically distributed random variables to conditional cases.  J Korea Math Soc, 52(2):431–445, 2015.  doi:10.4134/JKMS.2015.52.2.431.

Appendices

Appendix A contains an example illustrating several concepts discussed in
Sections 3.3–3.4. The other appendices contain the proofs of Theorems 1–9
given in the main text. In particular, Appendix B contains notation used
throughout Appendices C–I, Appendix C contains the proofs for the DaC-SMC
estimators’ consistency and $L^{p}$ inequalities (Theorems 1, 2, and 5),
Appendix D that for the unbiasedness of the unnormalized target estimators
(Theorem 3), Appendix E that for the bias bound of the normalized target ones
(Theorem 4), Appendix F that for the asymptotic normality (Theorem 6),
Appendix G that for the product $L^{2}$ inequality (Theorem 7), Appendix H the
derivation of the locally optimal proposals and auxiliary measures (Theorem
8), and Appendix I the argument showing that the asymptotic variances of DaC-
SMC estimators are bounded above by those of their SMC counterparts (Theorem
9).

##  Appendix A A further example

Consider the following time-varying version of the toy hierarchical model in
[37, Section 3.1]:

$Y_{t,l}\sim\mathcal{N}(X_{t,l},1),\enskip
X_{t,l}\sim\mathcal{N}(0,\theta_{t}),\enskip\forall
l\in[L],\quad\theta_{t}\sim f(d\theta;\theta_{t-1}),\enskip\forall
t\in[T],\enskip\theta_{0}\sim f(d\theta;1),$

where $T,L>0$ denote integers, $Y_{0:T,1:L}:=((Y_{t,l})_{l=1}^{L})_{t=0}^{T}$
some observed variables, $X_{0:T,1:L}$ $:=((X_{t,l})_{l=1}^{L})_{t=0}^{T}$
some latent variables, $\theta_{0:T}=(\theta_{t})_{t=0}^{T}$ some unknown
parameters, and $f(d\theta;\eta)$ the exponential distribution with mean
$\eta$. Suppose we observe data $y_{0:T,1:L}$ and wish to draw inferences from
the smoothing distribution proportional to

$\displaystyle\rho(d\theta_{0:T},dx_{0:T,1:L}):=$ $\displaystyle
f(d\theta_{0};1)\left(\prod_{t=1}^{T}f(d\theta_{t};\theta_{t-1})\right)$
$\displaystyle\times\left(\prod_{t=0}^{T}\prod_{l=1}^{L}\mathcal{N}(y_{t,k};x_{t,l},1)\mathcal{N}(dx_{t,l};0,\theta_{t})\right).$

We consider two approaches with which we could tackle this problem:

###  A.1 The DaC-SMC approach

One way to obtain tractable approximations of the posterior $\rho$ is to apply
DaC-SMC with the following tree. Each parameter variable $\theta_{t}$ has a
node $t$ assigned to it (with $E_{t}=[0,\infty)$), each latent variable
$X_{t,l}$ has a node $(t,l)$ assigned to it (with $E_{t,l}=\mathbb{R}$), and
the tree has an extra ‘dummy’ root node $T+1$ (which we introduce purely for
notational convenience). The nodes assigned to $\theta_{0}$ and $X_{1:T,1:L}$
are the tree’s leaves and, for each $t$ in $[T+1]$, the children of $t$ are
$t-1$ and $(t-1,1),\dots,(t-1,L)$. For the auxiliary measure $\gamma_{t_{-}}$
associated with the node $t$ assigned to $\theta_{t}$, we simply remove all
terms in $\rho$ that involve variables $\theta_{s}$ and
$x_{s,1},\dots,x_{s,L}$ with $s\geq t$:

$\displaystyle\gamma_{(t+1)_{-}}(d\theta_{0:t},dx_{0:t,1:L}):=$ $\displaystyle
f(d\theta_{0};1)\left(\prod_{t^{\prime}=1}^{t}f(d\theta_{t^{\prime}};\theta_{t^{\prime}-1})\right)$
$\displaystyle\times\left(\prod_{t^{\prime}=0}^{t}\prod_{l=1}^{L}\mathcal{N}(y_{t^{\prime},l};x_{t^{\prime},l},1)\mathcal{N}(dx_{t^{\prime},l};0,\theta_{t^{\prime}})\right)\quad\forall
t=0,\dots,T.$

(Note that $\gamma_{(T+1)_{-}}$ coincides with the unnormalized smoothing
distribution $\rho$, so $\pi_{(T+1)_{-}}^{N}$ approximates the smoothing
distribution.) For the proposal kernels we pick

$\displaystyle K_{0}(d\theta_{0})$ $\displaystyle:=f(d\theta_{0};1),\quad
K_{t}(\theta_{t-1},d\theta_{t}):=f(d\theta_{t};\theta_{t-1})\enskip\forall
t=1,\dots,T,$ (23) $\displaystyle K_{t,k}(dx_{t,l})$
$\displaystyle:=\mathcal{N}(dx_{t,l};y_{t,l},1)\enskip\forall
t=0,\dots,T,\enskip l=1,\dots,L;$ (24)

in which case

$w_{(t+1)_{-}}(\theta_{t},x_{t,1:L})=\prod_{l=1}^{L}\mathcal{N}(x_{t,l};0,\theta_{t-1})\quad\forall
t=0,\dots,T.$

Hence, with
$g_{t,l}(\theta):=\sum_{n_{l}=1}^{N}\mathcal{N}(X_{t,l}^{n_{l}};0,\theta)$ and
$g_{t}(\theta):=\prod_{l=1}^{L}g_{t,l}(\theta)$,

$\displaystyle\frac{N^{L+1}}{\mathcal{Z}_{\mathcal{C}_{t}}^{N}}\gamma_{(t+1)_{-}}^{N}$
$\displaystyle=\sum_{n=1}^{N}\left(\prod_{l=1}^{L}\sum_{n_{l}=1}^{N}\mathcal{N}(X_{t,l}^{n_{l}};0,\theta_{t}^{n})\delta_{X_{t,l}^{n_{l}}}\right)\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})}$
(25)
$\displaystyle=\sum_{n=1}^{N}\left(\prod_{l=1}^{L}\sum_{n_{l}=1}^{N}\frac{\mathcal{N}(X_{t,l}^{n_{l}};0,\theta_{t}^{n})}{g_{t,l}(\theta_{t}^{n})}\delta_{X_{t,l}^{n_{l}}}\right)g_{t}(\theta_{t}^{n})\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})}.$

(Note that we are omitting the $N$ superscripts from $\theta^{n}_{t}$ and
$X^{n}_{t,l}$ to simplify the notation.) Setting
$\omega_{t,l}^{n,n_{l}}:=\frac{\mathcal{N}(X_{t,l}^{n_{l}};0,\theta_{t}^{n})}{g_{t,l}(\theta_{t}^{n})}$
and
$\omega_{t}^{n}:=\frac{g_{t}(\theta_{t}^{n})}{\sum_{m=1}^{N}g_{t}(\theta_{t}^{m})}$,
and normalizing, we find that

$\pi_{(t+1)_{-}}^{N}=\sum_{n=1}^{N}\left(\prod_{l=1}^{L}\sum_{n_{l}=1}^{N}\omega_{t,l}^{n,n_{l}}\delta_{X_{t,l}^{n_{l}}}\right)\omega_{t}^{n}\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})}.$

Hence, we are able to draw $N$ samples from $\pi_{(t+1)_{-}}^{N}$ in
$\mathcal{O}(N^{2})$ operations by: $(a)$ evaluating
$(\omega_{t}^{n})_{n=1}^{N}$; $(b)$ drawing $N$ indices $(m_{n})_{n=1}^{N}$
from $(\omega_{t}^{n}\delta_{n})_{n=1}^{N}$ using the alias method; $(c)$ for
each $n$ and $l$, evaluating $(\omega_{t,l}^{m_{n},k})_{k=1}^{N}$; $(d)$ for
each $n$ and $l$, drawing an index $k_{n,l}$ from
$(\omega_{t}^{m_{n},k}\delta_{k})_{k=1}^{N}$; and $(e)$ setting

$\bm{X}_{(t+1)_{-}}^{n}:=(\theta_{1:t}^{m_{n}},X_{1:t-1,1:L}^{m_{n}},X_{t,1}^{k_{n,1}},\dots,X_{t,L}^{k_{n,L}})\quad\forall
n\in[N].$

This process can be trivially sped up by storing and re-using the weights
$(\omega_{t}^{m_{n},k})_{k=1}^{N}$ for indices $m_{n}$ drawn more than once.
If there are $m_{n}$s drawn many times, one may be better off computing
probability and alias tables for the corresponding weights
$(\omega_{t}^{m_{n},k})_{k=1}^{N}$ and using these to generate the
corresponding $k_{n}$s.

###  A.2 An alternative approach

This model has two types of direct dependencies linking its variables:
$\theta_{t+1}$ is drawn from an exponential distribution with mean
$\theta_{t}$ and $X_{t,l}$ is drawn from a zero-mean Gaussian distribution
with variance $\theta_{t}$. The first of these two is encoded into the kernels
in (23,24): we propose a value for $\theta_{t+1}^{n}$ by sampling
$f(d\theta_{t+1};\theta_{t}^{n})$. The second type, however, is not:
$X_{t,l}^{n}$ are drawn independently of $\theta_{t}^{n}$ by sampling from
$\mathcal{N}(y_{t,l},1)$.

To account for the second type of dependency, we could instead try drawing

$\theta^{n}_{t}\sim f(d\theta_{t};\theta_{t-1}^{n})\enskip\text{(or
$f(d\theta_{0};1)$ if $t=0$)},\quad
X_{t,1:L}^{n,m}\sim\prod_{l=1}^{L}\mathcal{N}\left(\frac{\theta_{t}^{n}}{1+\theta_{t}^{n}}y_{t,l},\frac{\theta_{t}^{n}}{1+\theta_{t}^{n}}\right),$

for each for all $n$ and $m$ in $[N]$; resulting in a total cost of
$\mathcal{O}(N^{2})$. We would then have to replace the product-form
approximation to $\gamma_{(t+1)_{-}}$ in (25) with the following ‘partially
product-form’ one [37, Section 3.2]:

$\displaystyle\gamma_{(t+1)_{-}}^{N}:$
$\displaystyle=\frac{\mathcal{Z}_{t}^{N}}{N^{L+1}}\sum_{n=1}^{N}\left[\prod_{l=1}^{L}\sum_{n_{l}=1}^{N}\mathcal{N}(y_{t,l};0,1+\theta_{t}^{n})\delta_{X^{n,n_{l}}_{t,l}}\right]\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})}$
$\displaystyle=\frac{\mathcal{Z}_{t}^{N}}{N^{L+1}}\sum_{n=1}^{N}\left[\prod_{l=1}^{L}\frac{1}{N}\sum_{n_{l}=1}^{N}\delta_{X^{n,n_{l}}_{t,l}}\right]g_{t}(\theta_{t}^{n})\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})},$

where $g_{t}(\theta):=\prod_{l=1}^{L}\mathcal{N}(y_{t,l};0,1+\theta)$ and
$\mathcal{Z}_{t}^{N}$ denotes the mass of the previously computed
$\gamma_{t_{-}}^{N}$ (with $\mathcal{Z}_{0}^{N}:=1$). Drawing $N$ samples from
the approximation’s normalization,

$\pi_{(t+1)_{-}}^{N}:=\frac{\gamma_{(t+1)_{-}}^{N}}{\mathcal{Z}_{(t+1)_{-}}^{N}}=\sum_{n=1}^{N}\left[\prod_{l=1}^{L}\frac{1}{N}\sum_{n_{l}=1}^{N}\delta_{X^{n,n_{l}}_{t,l}}\right]\omega_{t}^{n}\delta_{(\theta_{1:t}^{n},X^{n}_{1:t-1,1:L})}$

where
$\omega_{t}^{n}:=g_{t}(\theta_{t}^{n})/\sum_{m=1}^{N}g_{t}(\theta_{t}^{m})$,
can then be done in $\mathcal{O}(N)$ operations by: $(a)$ evaluating
$(\omega_{t}^{n})_{n=1}^{N}$; $(b)$ drawing $N$ indices $(m_{n})_{n=1}^{N}$
from $(\omega_{t}^{n}\delta_{n})_{n=1}^{N}$ using the alias method; $(c)$ for
each $n$ and $l$, drawing an index $k_{n,l}$ from
$(N^{-1}\delta_{k})_{k=1}^{N}$ ; and $(d)$ setting

$\bm{X}_{(t+1)_{-}}^{n}:=(\theta_{1:t}^{m_{n}},X_{1:t-1,1:L}^{m_{n}},X_{t-1,1}^{m_{n},k_{n,1}},\dots,X_{t-1,L}^{m_{n},k_{n,L}})\quad\forall
n\in[N].$

##  Appendix B Important notation for the proofs

Throughout Appendices C–I, we employ the notation introduced in Sections
3.1–3.2 and the following:

  * •

Unless specified otherwise, $(\Omega,\mathcal{F},\mathbb{P})$ denotes the
underlying probability space on which all random variables introduced in
Algorithm 3 (including those indexed by different $N$s) are jointly defined,
and $\mathbb{E}$ denotes expectations with respect to $\mathbb{P}$.

  * •

$C$ denotes a generic positive constant dependent on $u$, $v$, and $p$ but not
on $\varphi$ and $N$ that may change from one line to the next. We use this
notation whenever we are interested only in the existence of such a constant
and its precise value is unimportant. To not overly repeat ourselves,
throughout the appendices, we often introduce $C$ without mention (e.g. Lemma
1 below) and leave the corresponding ‘there exists some $C<\infty$ such that’
or ‘for some $C<\infty$’ statements implicit.

  * •

Given any two subsets $A\subseteq B\subseteq\mathbb{T}$, a measure $\nu$ on
$(E_{A},\mathcal{E}_{A})$, and measurable function $\psi$ on
$(E_{B},\mathcal{E}_{B})$, we use $\nu(\psi)$ to denote the measurable
function on $(E_{B\backslash A},\mathcal{E}_{B\backslash A})$ obtained by
integrating the arguments of $\psi$ indexed by $A$ with respect to $\nu$:

$\nu(\psi)(x_{B\backslash A}):=\int\psi(x_{A},x_{B\backslash
A})\nu(dx_{A})\quad\forall x_{B\backslash A}\in E_{B\backslash A},$

under the assumption that the integral is well-defined for all $x_{B\backslash
A}$ in $E_{B\backslash A}$.

  * •

For any $u$ in $\mathbb{T}^{\not\partial}$ and $N>0$,
$\epsilon_{u}^{N}:=N^{-1}\sum_{n=1}^{N}\delta_{\bm{X}^{n,N}_{u_{-}}}$ denotes
the empirical distribution of the resampled particles (not to be confused with
$\pi_{u}^{N}$ in (10) which is empirical distribution of these particles after
mutation) and
$\pi_{\mathcal{C}_{u}}^{N}:=\prod_{v\in\mathcal{C}_{u}}\pi_{v}^{N}$ denotes
the product of the particle approximations to the extended normalized
auxiliary measures indexed by $u$’s children. Additionally,
$\mathcal{F}_{\mathcal{C}_{u}}^{N}$ denotes the sigma algebra generated by the
mutated particles
$(\bm{X}_{\mathcal{C}_{u}}^{n,N})_{n=1}^{N}:=(\bm{X}_{u1}^{n,N},\dots,\bm{X}_{uc_{u}}^{n,N})_{n=1}^{N}$
indexed by $u$’s children, and that $\mathcal{F}_{u_{-}}^{N}$ that generated
by the resampled particles $(\bm{X}_{u_{-}}^{n,N})_{n=1}^{N}$ indexed by $u$.

##  Appendix C Proofs of Theorems 1, 2, and 5

Throughout this appendix, we use the notation described in Appendix B. The aim
of the appendix is to prove $L^{p}$ inequalities in Theorem 5. The laws of
large numbers in Theorem 1 follow from the inequalities and the ensuing
consequence of the Borel-Cantelli lemma:

######  Lemma 1 (e.g. [6, p. 17]).

Let $(Z_{N})_{N=1}^{\infty}$ be a sequence of real-valued random variables. If
there exists a $p>2$ such that

$\mathbb{E}\left[\left|Z_{N}\right|^{p}\right]\leq\frac{C}{N^{p/2}}\quad\forall
N>0,$

then $Z_{N}\to 0$ with probability one as $N\to\infty$.

The weak convergence in Theorem 2 follows directly:

###### Proof of Theorem 2.

Because the product of (finitely-many) Borel $\sigma$-algebras on Polish
spaces coincides with the Borel $\sigma$-algebra on the product space, as the
former are separable cf. [32, Lemma 1.2], and the latter is also Polish, this
follows immediately from Theorem 1 and the equivalence established in [5,
Theorem 2.2]. See [53, Section S1, Theorem 1] for details in an SMC context. ∎

To establish the $L^{p}$ inequalities in Theorem 5, we first obtain such
results for $(\pi_{u}^{N})_{u\in\mathbb{T}}$:

######  Theorem 10 ($L^{p}$ inequalities for
$(\pi_{u}^{N})_{u\in\mathbb{T}}$).

If Assumptions 1–2 hold, then

$\displaystyle\mathbb{E}\left[\left|\pi^{N}_{u}(\varphi)-\pi_{u}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}},\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}),\enskip p\geq 1,\enskip
u\in\mathbb{T}.$

To prove the above we follow the approach for standard SMC in [14] adapted to
cover the case of general $p$ (as in, for example, [44]). In particular, we
derive an $L^{p}$ inequality for the empirical distributions of the particles
indexed by the tree’s leaves and show that each step of the algorithm
preserves this inequality. Let’s start: the Marcinkiewicz-Zygmund-type
inequality in [17, Lemma 7.3.3] shows that

$\mathbb{E}\left[\left|\pi^{N}_{u}(\varphi)-\pi_{u}(\varphi)\right|^{m}\right]^{\frac{1}{m}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}),$

for all $u$ in $\mathbb{T}^{\partial}$ and positive integers $m$. We extend
the above to all $p\geq 1$ using Jensen’s inequality: if $m$ is the smallest
integer no smaller than $p$, then

$\displaystyle\mathbb{E}\left[\left|\pi^{N}_{u}(\varphi)-\pi_{u}(\varphi)\right|^{p}\right]^{\frac{1}{p}}$
$\displaystyle=\mathbb{E}\left[\left(\left|\pi^{N}_{u}(\varphi)-\pi_{u}(\varphi)\right|^{m}\right)^{\frac{p}{m}}\right]^{\frac{1}{p}}\leq\mathbb{E}\left[\left|\pi^{N}_{u}(\varphi)-\pi_{u}(\varphi)\right|^{m}\right]^{\frac{1}{m}}$
$\displaystyle\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{u}).$ (26)

Taking the product of $\pi_{v}^{N}$ over all children $v$ of a node $u$
preserves (26):

######  Lemma 2 (Product step).

If, in addition to Assumptions 1–2, (26) is satisfied for each child $v$ (i.e.
it holds with $v$ replacing $u$ therein) of a node $u$ in
$\mathbb{T}^{\not\partial}$ and some $p\geq 1$, then

$\mathbb{E}\left[\left|\pi^{N}_{\mathcal{C}_{u}}(\varphi)-\pi_{\mathcal{C}_{u}}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{\mathcal{C}_{u}}).$ (27)

The key to proving the above is following lemma giving an $L^{p}$ inequality
for products of independent random probability measures satisfying their own
individual $L^{p}$ inequalities:

######  Lemma 3.

Suppose that $(\eta_{1}^{N})_{N=1}^{\infty}$ and
$(\eta_{2}^{N})_{N=1}^{\infty}$ are independent sequences of random
probability measures defined on some common probability triplet
$(\Omega,\mathcal{F},\mathbb{P})$, respectively taking values in some
measurable spaces $(S_{1},\mathcal{S}_{1})$ and $(S_{2},\mathcal{S}_{2})$, and
satisfying

$\mathbb{E}\left[\left|\eta^{N}_{k}(\varphi)-\eta_{k}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(S_{k}),\enskip k\in\\{1,2\\},$ (28)

for some $p\geq 1$, with limits $\eta_{1}$ and $\eta_{2}$ that are also
probability measures. Then,

$\mathbb{E}\left[\left|(\eta^{N}_{1}\times\eta^{N}_{2})(\varphi)-(\eta_{1}\times\eta_{2})(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(S_{1}\times S_{2}).$ (29)

###### Proof.

Fix any $N>0$ and $\varphi$ in $\mathcal{B}_{b}(S_{1}\times S_{2})$ and note
that

$\displaystyle\eta^{N}_{1}\times\eta^{N}_{2}(\varphi)-\eta_{1}\times\eta_{2}(\varphi)=$
$\displaystyle(\eta^{N}_{1}-\eta_{1})\times(\eta^{N}_{2}-\eta_{2})(\varphi)$
(30)
$\displaystyle+(\eta^{N}_{1}-\eta_{1})\times\eta_{2}(\varphi)+\eta_{1}\times(\eta^{N}_{2}-\eta_{2})(\varphi).$

Hence, Minkowski’s inequality implies that

$\displaystyle\mathbb{E}\left[\left|\eta^{N}_{1}\times\eta^{N}_{2}(\varphi)-\eta_{1}\times\eta_{2}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq$
$\displaystyle\mathbb{E}\left[\left|(\eta^{N}_{1}-\eta_{1})\times(\eta^{N}_{2}-\eta_{2})(\varphi)\right|^{p}\right]^{\frac{1}{p}}$
(31)
$\displaystyle+\mathbb{E}\left[\left|\eta^{N}_{1}(\eta_{2}(\varphi))-\eta_{1}(\eta_{2}(\varphi))\right|^{p}\right]^{\frac{1}{p}}$
$\displaystyle+\mathbb{E}\left[\left|\eta^{N}_{2}(\eta_{1}(\varphi))-\eta_{2}(\eta_{1}(\varphi))\right|^{p}\right]^{\frac{1}{p}}.$

Because $\eta_{1}$ and $\eta_{2}$ are probability measures,
$\eta_{1}(\varphi)$ and $\eta_{2}(\varphi)$ are bounded above by
$\left|\left|{\varphi}\right|\right|$. For this reason, the $L^{p}$
inequalities in (28) imply that

$\displaystyle\mathbb{E}\left[\left|\eta^{N}_{1}(\eta_{2}(\varphi))-\eta_{1}(\eta_{2}(\varphi))\right|^{p}\right]^{\frac{1}{p}}$
$\displaystyle\leq\frac{C\left|\left|{\eta_{2}(\varphi)}\right|\right|}{N^{1/2}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}},$
(32)
$\displaystyle\mathbb{E}\left[\left|\eta^{N}_{2}(\eta_{1}(\varphi))-\eta_{2}(\eta_{1}(\varphi))\right|^{p}\right]^{\frac{1}{p}}$
$\displaystyle\leq\frac{C\left|\left|{\eta_{1}(\varphi)}\right|\right|}{N^{1/2}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}.$
(33)

To control the remaining term in (31), let $\mathcal{F}_{2}$ denote the
$\sigma$-algebra generated by the $\eta_{2}^{N}$s. For each $\omega$ in
$\Omega$,
$x_{1}\mapsto\int\varphi(x_{1},x_{2})\eta_{2}^{N}(\omega,dx_{2})-\eta_{2}(\varphi)(x_{1})$
is a bounded function on $S_{1}$. Hence, (28) and the independence of
$(\eta_{1}^{N})_{N=1}^{\infty}$ and $(\eta_{2}^{N})_{N=1}^{\infty}$ imply that

$\displaystyle\mathbb{E}\left[\left|(\eta^{N}_{1}-\eta_{1})\times(\eta^{N}_{2}-\eta_{2})(\varphi)\right|^{p}\right]^{\frac{1}{p}}$
(34)
$\displaystyle\qquad\qquad=\mathbb{E}\left[\left|\eta^{N}_{1}(\eta^{N}_{2}(\varphi)-\eta_{2}(\varphi))-\eta_{1}(\eta^{N}_{2}(\varphi)-\eta_{2}(\varphi))\right|^{p}\right]^{\frac{1}{p}}$
$\displaystyle\qquad\qquad=\mathbb{E}\left[\mathbb{E}\left[\left|\eta^{N}_{1}(\eta^{N}_{2}(\varphi)-\eta_{2}(\varphi))-\eta_{1}(\eta^{N}_{2}(\varphi)-\eta_{2}(\varphi))\right|^{p}|\mathcal{F}_{2}\right]\right]^{\frac{1}{p}}$
$\displaystyle\qquad\qquad\leq\frac{C\mathbb{E}\left[\left|\left|{\eta^{N}_{2}(\varphi)-\eta_{2}(\varphi)}\right|\right|^{p}\right]^{\frac{1}{p}}}{N^{1/2}}\leq\frac{C\mathbb{E}\left[(\left|\left|{\eta^{N}_{2}(\varphi)}\right|\right|+\left|\left|{\eta_{2}(\varphi)}\right|\right|)^{p}\right]^{\frac{1}{p}}}{N^{1/2}}$
$\displaystyle\qquad\qquad\leq\frac{C\mathbb{E}\left[(2\max\\{\left|\left|{\eta^{N}_{2}(\varphi)}\right|\right|,\left|\left|{\eta_{2}(\varphi)}\right|\right|\\})^{p}\right]^{\frac{1}{p}}}{N^{1/2}}\leq\frac{2C\left|\left|{\varphi}\right|\right|}{N^{1/2}}.$

Putting (31–34) together, we obtain the $L^{p}$ inequality (29) for the
product. ∎

###### Proof of Lemma 2.

Because
$(\pi^{N}_{u1})_{N=1}^{\infty},\dots,(\pi^{N}_{uc_{u}})_{N=1}^{\infty}$ are
independent sequences of probability measures by construction, so are
$(\pi^{N}_{u[k]})_{N=1}^{\infty}$ and $(\pi^{N}_{u(k+1)})_{N=1}^{\infty}$ for
all $k<c_{u}$, where $u[k]:=\\{u1,\dots,uk\\}$ denotes the set containing the
first $k$ children of $u$ and $\pi^{N}_{u[k]}:=\prod_{v\in u[k]}\pi^{N}_{v}$
the corresponding product of $\pi_{v}^{N}$s. Hence, starting from the premise
and repeatedly applying Lemma 3, we obtain the $L^{p}$ inequality for
$\pi^{N}_{u[c_{u}]}=\pi_{\mathcal{C}_{u}}^{N}$. ∎

Emulating the approach of [14, Lemma 4] and [44, Lemma 1], we find that the
correction step also respects the inequality:

######  Lemma 4 (Correction step).

If, in addition to Assumptions 1–2, (27) is satisfied for some $u$ in
$\mathbb{T}^{\not\partial}$ and $p\geq 1$, then

$\mathbb{E}\left[\left|\pi_{u_{-}}^{N}(\varphi)-\pi_{u_{-}}(\varphi)\right|^{p}\right]^{\frac{1}{p}}\leq\frac{C\left|\left|{\varphi}\right|\right|}{N^{1/2}}\quad\forall
N>0,\enskip\varphi\in\mathcal{B}_{b}(\bm{E}_{\mathcal{C}_{u}}).$ (35)

###### Proof.

Fix any $p\geq 1$, $N>0$, and $\varphi$ in
$\mathcal{B}_{b}(\bm{E}_{\mathcal{C}_{u}})$. Recall the definitions in Section
3.2:

$\displaystyle\gamma_{u_{-}}^{N}(d\bm{x}_{\mathcal{C}_{u}})=w_{u_{-}}(\bm{x}_{\mathcal{C}_{u}})\gamma_{\mathcal{C}_{u}}^{N}(d\bm{x}_{\mathcal{C}_{u}}),\quad\gamma_{\mathcal{C}_{u}}^{N}(d\bm{x}_{\mathcal{C}_{u}})$
$\displaystyle=\mathcal{Z}_{\mathcal{C}_{u}}^{N}\pi_{\mathcal{C}_{u}}^{N}(d\bm{x}_{\mathcal{C}_{u}}).$

It follows that

$\pi^{N}_{u_{-}}(\varphi)=\frac{\gamma_{u_{-}}^{N}(\varphi)}{\gamma_{u_{-}}^{N}(\bm{E}_{\mathcal{C}_{u}})}=\frac{\gamma^{N}_{\mathcal{C}_{u}}(w_{u_{-}}\varphi)}{\gamma^{N}_{\mathcal{C}_{u}}(w_{u_{-}})}=\frac{\pi^{N}_{\mathcal{C}_{u}}(w_{u_{-}}\varphi)}{\pi^{N}_{\mathcal{C}_{u}}(w_{u_{-}})}.$

Because $w_{u_{-}}=d\gamma_{u_{-}}/d\gamma_{\mathcal{C}_{u}}$, we similarly
have that

$\pi_{u_{-}}(\varphi)=\frac{\gamma_{u_{-}}(\varphi)}{\gamma_{u_{-}}(\bm{E}_{\mathcal{C}_{u}})}=\frac{\gamma_{\mathcal{C}_{u}}(w_{u_{-}}\varphi)}{\gamma_{\mathcal{C}_{u}}(w_{u_{-}})}=\frac{\pi_{\mathcal{C}_{u}}(w_{u_{-}}\varphi)}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}.$

Hence,

$\displaystyle\left|\pi_{u_{-}}^{N}(\varphi)-\pi_{u_{-}}(\varphi)\right|\leq\left|\pi_{u_{-}}^{N}(\varphi)-\frac{\pi_{\mathcal{C}_{u}}^{N}(w_{u_{-}}\varphi)}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}\right|+\left|\frac{\pi_{\mathcal{C}_{u}}^{N}(w_{u_{-}}\varphi)}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}-\frac{\pi_{\mathcal{C}_{u}}(w_{u_{-}}\varphi)}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}\right|.$
(36)

To control the first term on the right-hand side, we use

$\displaystyle\left|\pi_{u_{-}}^{N}(\varphi)-\frac{\pi_{\mathcal{C}_{u}}^{N}(w_{u_{-}}\varphi)}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}\right|$
$\displaystyle\leq\frac{\left|\pi_{u_{-}}^{N}(\varphi)\right|\left|\pi_{\mathcal{C}_{u}}(w_{u_{-}})-\pi_{\mathcal{C}_{u}}^{N}(w_{u_{-}})\right|}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}$
(37)
$\displaystyle\leq\frac{\left|\left|{\varphi}\right|\right|\left|\pi_{\mathcal{C}_{u}}(w_{u_{-}})-\pi_{\mathcal{C}_{u}}^{N}(w_{u_{-}})\right|}{\pi_{\mathcal{C}_{u}}(w_{u_{-}})}.$

Because,