# A divide and conquer sequential Monte Carlo approach to high dimensional
filtering

Francesca R. Crucinio Department of Statistics, University of Warwick,
Coventry, CV4 7AL, UK. Email: $\\{$francesca.crucinio,
<EMAIL_ADDRESS>Adam M. Johansen11footnotemark: 1

###### Abstract

We propose a divide-and-conquer approach to filtering which decomposes the
state variable into low-dimensional components to which standard particle
filtering tools can be successfully applied and recursively merges them to
recover the full filtering distribution. It is less dependent upon
factorization of transition densities and observation likelihoods than
competing approaches and can be applied to a broader class of models.
Performance is compared with state-of-the-art methods on a benchmark problem
and it is demonstrated that the proposed method is broadly comparable in
settings in which those methods are applicable, and that it can be applied in
settings in which they cannot.

##  1 Introduction

Particle filters (PFs), an instance of sequential Monte Carlo (SMC) methods,
are a popular class of algorithms to perform state estimation for state space
models (SSM) — or general-state-space hidden Markov models as they are
sometimes known. We consider the class of SSMs with a latent d-valued process
$(X_{t})_{t\geq 1}$ and conditionally independent p-valued observations
$(Y_{t})_{t\geq 1}$. Such a SSM $(X_{t},Y_{t})_{t\geq 1}$ is defined by the
transition density $f_{t}(x_{t-1},x_{t})$ of the latent process, with the
convention that $f_{1}(x_{0},x_{1})\equiv f_{1}(x_{1})$, and the observation
likelihood $g_{t}(y_{t}|x_{t})$. In this work, we are interested in
approximating the sequence of filtering distributions,
$(p(x_{t}|y_{1:t}))_{t\geq 1}$, i.e. at each time $t$ the distribution of the
latent state at that time given the observations obtained by that time.

Basic PF algorithms are known to suffer from the curse of dimensionality,
requiring an exponential increase in computational requirements as the
dimension $d$ grows, limiting their applicability to large systems (Rebeschini
and Van Handel, 2015; Bengtsson et al., 2008). While the ensemble Kalman
filter (Evensen, 2009) can tackle high dimensional problems, it involves
approximations which do not disappear even asymptotically and does not perform
well if the model is far from linear and Gaussian (Lei et al., 2010).

To extend the use of particle filters to high dimensional problems it is
natural to attempt to exploit the fact that dependencies in high dimensional
SSMs encountered in practice are often local in space in order to decompose
the filtering problem into a collection of local low-dimensional problems
which can somehow be combined; examples of this strategy are the block
particle filter (Rebeschini and Van Handel, 2015), space-time particle filters
(STPF; Beskos et al. (2017)) and nested sequential Monte Carlo (NSMC) methods
(Næsseth et al., 2015, 2019).

We propose a divide-and-conquer approach in which the state space is divided
into smaller subsets over which standard particle filtering ideas can be
applied, these smaller subsets are then recursively merged in a principled
manner to obtain approximations over the full state space. Our method is an
extension of the divide-and-conquer sequential Monte Carlo (DaC-SMC) algorithm
introduced by Lindsten et al. (2017) to the filtering context, which exploits
ideas akin those in Klaas et al. (2005); Lin et al. (2005) to marginalize out
the past $x_{1:t-1}$ at a given time $t$.

In order to apply DaC-SMC to the filtering problem, we define a non-standard
sequence of targets evolving both in space and in time: at a given time $t$ we
define $d$ univariate targets serving as proxies for the marginals of the
filtering distribution, $p(x_{t}(i)|y_{1:t})$ for $i=1,\dots,d$, we then
iteratively combine these lower dimensional targets to obtain approximations
of higher-dimensional marginals (e.g. $p(x_{t}(i:i+1)|y_{1:t})$) until we
recover the full filtering distribution, $p(x_{t}|y_{1:t})\equiv
p(x_{t}(1:d)|y_{1:t})$.

Unlike NSMC and STPF, this approach does not require analytic expressions for
marginals of the transition density or observation likelihood, but only point-
wise evaluations of $f_{t}$ and $g_{t}$, making it suitable to tackle high
dimensional SSM in which the observations are correlated in non-trivial ways
(cf. the model in Section 4.2), as is common in real applications (see, e.g.,
Chib et al. (2009, Section 2)).

We review the basic ideas of particle filtering, its marginal variant and
divide-and-conquer SMC in Section 2; we then show how to extend the ideas
underlying DaC-SMC to the filtering problem in Section 3, where we also
discuss strategies to improve computational cost and accuracy. Finally, in
Section 4 we compare the performances of the divide-and-conquer approach with
NSMC and STPF on a simple linear Gaussian SSM for which the Kalman filter
provides the exact filtering distribution. Our experiments show that the
errors in approximating the true filtering distribution obtained with the
divide-and-conquer approach are comparable to those of NSMC and STPF. We then
consider a spatial model whose correlation structure in the observation model
makes it impossible to apply NSMC or STPFs (at least without additional
approximations) and empirically show that the proposed approach can recover
stable estimates of the filtering distribution which, in small dimensional
settings, coincide with those obtained with a bootstrap PF with a large number
of particles.

##  2 Background

###  2.1 Particle Filtering

We describe here the basic SMC approach, often referred to as sequential
importance resampling (Doucet and Johansen, 2011, p. 15), and refer to Liu
(2001); Chopin and Papaspiliopoulos (2020) for a more extensive treatment.

Given the sequence of unnormalized target densities $(\gamma_{t})_{t\geq 1}$,
with

$\displaystyle\gamma_{t}(x_{1:t})=p(x_{1:t},y_{1:t})=\prod_{k=1}^{t}f_{k}(x_{k-1},x_{k})g_{k}(y_{k}|x_{k}),$
(1)

defined on $({}^{d})^{t}$, PFs proceed iteratively, and, at time $t-1$
approximate

$\pi_{t-1}:=\gamma_{t-1}/\int\gamma_{t-1}(x_{1:t-1})\mathrm{d}x_{1:t-1}$

with a cloud of particles $\\{x_{1:t-1}^{n}\\}_{n=1}^{N}$. The particles are
propagated forward in time using a Markov kernel $K_{t}(x_{1:t-1},\cdot)$,
reweighted using the weight function $w_{t}:=\gamma_{t}/\gamma_{t-1}\otimes
K_{t}$ and resampled to obtain a new particle population
$\\{x_{1:t}^{n}\\}_{n=1}^{N}$ approximating $\pi_{t}$.

Standard PFs formally target distributions (1) whose dimension increases at
every time step $t$, although one is often interested only in the final time
marginal of the distributions so approximated, in this case the filtering
distribution. An alternative to this approach is given by marginal particle
filters (MPFs; Klaas et al. (2005)) and the closely related ideas of Lin et
al. (2005). MPFs target the filtering distribution directly

$\displaystyle\gamma_{t}(x_{t})=p(x_{t}|y_{1:t})=g_{t}(y_{t}|x_{t})\int
f_{t}(x_{t-1},x_{t})p(x_{t-1}|y_{1:t-1})\mathrm{d}x_{t-1};$ (2)

and, since the integral w.r.t. $x_{t-1}$ is intractable, replace
$p(x_{t-1}|y_{1:t-1})$ with its particle approximation obtained at time $t-1$.
Given the new sequence of targets, MPFs proceed as standard PFs, with the only
difference being that, whenever we need to compute an integral w.r.t.
$x_{t-1}$, this is approximated using $\pi_{t-1}^{N}$, obtained by normalizing
$\gamma_{t-1}^{N}$, the particle approximation of $p(x_{t-1}|y_{1:t-1})$.
Basic MPFs incur an $O(N^{2})$ cost for each time step, because of the
presence of the integral w.r.t. $x_{t-1}$ in the weight computations, although
lower cost strategies might be employed in some cases (Lin et al. (2005); see
also Klaas et al. (2006)).

###  2.2 Particle Filters for High Dimensional Problems

We briefly summarize three classes of particle filters which make use of space
decompositions to tackle the filtering problem which we believe to be the
state-of-the-art in Monte Carlo approximation of high dimensional filtering
distributions.

The block particle filter (BPF; Rebeschini and Van Handel (2015)) algorithm
relies on a decomposition of the state space d into lower dimensional blocks
on which, at each $t$, one step of a standard PF is run. The approximation of
the filtering distribution over the whole state space is obtained by taking
the product of the lower dimensional approximations on each block. BPFs are
inherently biased because of the decomposition into blocks; although this bias
can be eliminated asymptotically by allowing the blocks to grow at an
appropriate rate with computational effort.

Space-time particle filters (STPF; Beskos et al. (2017)) exploit local
dependence structures in the observation $y_{t}$ to gradually introduce the
likelihood term by decomposing the space dimension into smaller subsets and
running $N$ independent particle filters on each of the subsets (also called
islands) which are then combined using an importance resampling step, the
presence of the latter, guarantees asymptotically consistent approximations,
contrary to BPFs. Crucial to the implementation of STPF is that the joint law
(1) at time $t$ can be factorized so that the marginal of $x_{t}(i)$ given the
observations and the past only depends on a neighbourhood of $x_{t}(i)$,
$\\{x_{t}(j):\ j\in\mathcal{A}\\}$ for some
$\mathcal{A}\subset\\{1,\dots,d\\}$, for all $i=1,\dots,d$. STPFs are
particularly amenable to SSM which are time discretizations of SDEs, since in
this case one can build time discretization schemes which guarantee analytical
forms for the marginals (Akyildìz et al., 2022). A marginal version of STPFs
also exists (Beskos et al., 2017; Xu and Jasra, 2019).

Nested sequential Monte Carlo (NSMC; Næsseth et al. (2015)) treats the problem
of recovering the filtering distribution as a smoothing problem, where the
time variable is replaced by the dimension $d$, approximates the fully adapted
proposal of Pitt and Shephard (1999) with an inner SMC iteration and then uses
the result in the outer level which corresponds to a standard forward
filtering backward simulation algorithm. NSMC is particularly well-suited for
Markov random fields in which the temporal and the spatial components can be
separated, as this makes the backward simulation straightforward.

###  2.3 Divide and Conquer SMC

Divide-and-Conquer SMC (DaC-SMC; Lindsten et al. (2017)) is an extension of
standard SMC in which a collection of (unnormalized) target distributions
$(\gamma_{u})_{u\in\mathbb{T}}$ is indexed by the nodes of a rooted tree,
$\mathbb{T}$, and particles evolve from the leaves to the root,
$\mathfrak{R}$, rather than along a sequence of distributions indexed by
(algorithmic) time. It shares many of the convergence properties of standard
SMC (Kuntz et al., 2021).

The target distributions are defined on spaces whose dimension grows as we
progress up the tree: for each $u$, $\pi_{u}\propto\gamma_{u}$ is a density
over ${}^{|\mathbb{T}_{u}|}$ 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) and $|\mathbb{T}_{u}|$ denotes its
cardinality. We focus here on the case in which the state space is d, however,
essentially the same construction allows for much more general spaces,
including those with discrete components.

As in standard SMC, each distribution $\gamma_{u}$ is approximated by a
particle population $\\{x_{u}^{n}\\}_{n=1}^{N}$. However, these distributions
do not evolve ‘linearly’ but are merged whenever the corresponding branches of
$\mathbb{T}$ merge. For simplicity, we describe here the case in which
$\mathbb{T}$ is a binary tree and each non-leaf node $u$ has two children, a
left child $\ell(u)$ and a right child $r(u)$.

If $u$ is a leaf node, the algorithm performs a simple importance sampling
step with proposal $K_{u}$ and importance weight $w_{u}:=\gamma_{u}/K_{u}$ to
obtain a weighted particle population $\\{x_{u}^{n},w_{u}^{n}\\}_{n=1}^{N}$
approximating $\gamma_{u}$. Otherwise, to obtain a particle population
approximating $\gamma_{u}$, we gather the particle populations associated with
each of $u$’s children and compute the (weighted) product form estimator
(Kuntz et al., 2022)

$\displaystyle\gamma_{\mathcal{C}_{u}}^{N}=\frac{1}{N^{2}}\sum_{n_{1}=1}^{N}\sum_{n_{2}=1}^{N}w_{\ell(u)}(x_{\ell(u)}^{n_{1}})w_{r(u)}(x_{r(u)}^{n_{2}})\delta_{(x_{\ell(u)}^{n_{1}},x_{r(u)}^{n_{2}})}$
(3)

to approximate the product of the marginal distributions
$\gamma_{\mathcal{C}_{u}}:=\gamma_{\ell(u)}\times\gamma_{r(u)}$. The
$O(N^{2})$ cost of evaluating $\gamma_{\mathcal{C}_{u}}^{N}$ can be
prohibitively large, we discuss lower cost alternatives in Section 3.3.

We reweight the particle approximation $\gamma_{\mathcal{C}_{u}}^{N}$ of
$\gamma_{\mathcal{C}_{u}}$ to target $\gamma_{u}$; the resulting _mixture_
(importance) weights

$\displaystyle
m_{u}(x_{\ell(u)},x_{r(u)}):=\frac{\gamma_{u}(x_{\ell(u)},x_{r(u)})}{\gamma_{\ell(u)}(x_{\ell(u)})\gamma_{r(u)}(x_{r(u)})}$
(4)

capture the mismatch between $\gamma_{u}$ and
$\gamma_{\ell(u)}\times\gamma_{r(u)}$, and are incorporated prior to
resampling similarly to the auxiliary “twisting” function in the auxiliary PF
(see, e.g., Chopin and Papaspiliopoulos (2020, Chapter 10)). This leads to
weights of the form

$\displaystyle\tilde{w}_{u}(x_{\ell(u)},x_{r(u)}):=w_{\ell(u)}(x_{\ell(u)})w_{r(u)}(x_{r(u)})m_{u}(x_{\ell(u)},x_{r(u)}).$
(5)

Resampling $N$ times from $\tilde{w}_{u}\gamma_{\mathcal{C}_{u}}^{N}$, using
any unbiased resampling scheme (cf. Gerber et al. (2019)), we obtain an
equally weighted particle population
$\\{\tilde{x}_{u}^{n},w_{u}^{n}=1\\}_{n=1}^{N}$ approximating $\gamma_{u}$. If
necessary, we can then apply a $\pi_{u}$-invariant Markov kernel $K_{u}$, to
rejuvenate the particles. Algorithm 1, which is applied to the root node to
carry out the sampling process, summarizes this.

The DaC approach in Algorithm 1 is a special case of that considered in
Lindsten et al. (2017); Kuntz et al. (2021), in which the target at each non-
leaf node is defined on the product of the spaces on which each of its child
targets are defined. DaC-SMC is particularly amenable to distributed
implementation (Lindsten et al., 2017, Section 5.3).

1: if $u$ is a leaf then

2: Initialize: draw $x_{u}^{n}\sim K_{u}$ and compute
$w_{u}^{n}=\gamma_{u}/K_{u}$ for all $n\leq N$.

3: else

4: Recurse: set $(\\{x_{v}^{n},w_{v}^{n}\\}_{n=1}^{N}):=\text{dac\\_smc}(v)$
for $v$ in $\\{\ell(u),r(u)\\}$ and obtain $\gamma_{\mathcal{C}_{u}}^{N}$ in
(3).

5: Merge: compute $\tilde{w}_{u}^{(n_{1},n_{2})}$ in (5) for all
$n_{1},n_{2}\leq N$.

6: Resample: draw $\\{\tilde{x}_{u}^{n}\\}_{n=1}^{N}$ using weights
$\tilde{w}_{u}^{(n_{1},n_{2})}$ and set $w_{u}^{n}=1$ for all $n\leq N$.

7: (Optionally): draw $x_{u}^{n}\sim K_{u}(\tilde{x}_{u}^{n},\cdot)$ for all
$n\leq N$.
(Otherwise): set $x_{u}^{n}=\tilde{x}_{u}^{n}$ for all $n\leq N$.

8: end if

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

##  3 Divide and Conquer within Marginal SMC for Filtering

To apply the DaC-SMC algorithm described above to the filtering problem, we
need to identify a suitable collection
$(\widetilde{\gamma}_{t,u})_{u\in\mathbb{T}}$ indexed by the nodes of the tree
$\mathbb{T}$, for each time $t$. Graphically, this corresponds to a path graph
(corresponding to time) in which each node has associated with it a copy of
the tree $\mathbb{T}$ (corresponding to space). In this case Algorithm 1 takes
as input at the leaves a particle population approximating the filtering
distribution at time $t-1$ and outputs at the root a particle population
approximating the filtering distribution at time $t$.

At a given time $t$, to build the collection
$(\widetilde{\gamma}_{t,u})_{u\in\mathbb{T}}$, we consider spatial
decompositions of $x_{t}$ into low dimensional (often univariate) elements.
Here, we consider a simple decomposition obtained by identifying the $d$
components $(x_{t}(1),\dots,x_{t}(d))$ with the leaves of a tree $\mathbb{T}$.
As we move up the tree, the components are merged pairwise until
$x_{t}=x_{t}(1:d)$ is recovered at the root node $\mathfrak{R}$. For
simplicity, we assume that $d=2^{D}$ for some $D\in\mathbb{N}$ so that
$\mathbb{T}$ is a perfect binary tree; essentially the same construction
applies to general $d$.

We denote the set of components associated with node $u$ by $\mathcal{V}_{u}$,
its cardinality increases from leaves to root: at the level of the leaves
$|\mathcal{V}_{u}|=1$, while $|\mathcal{V}_{\mathfrak{R}}|=d$. Figure 1 shows
the space decomposition for $d=8$.

$x_{t}(1:8)$$x_{t}(1:4)$$x_{t}(1:2)$$x_{t}(1)$ Figure 1: Space decomposition
for $d=8$.

As the filtering problem has an inherent (temporal) sequential structure, the
collection $(\widetilde{\gamma}_{t,u})_{u\in\mathbb{T}}$ at time $t$ is most
easily specified in terms of the filtering distribution at time $t-1$,
$\widetilde{\gamma}_{t-1,\mathfrak{R}}$ — as shown below in (6). Similarly to
MPFs, we deal with this dependence by approximately marginalizing out the
previous timestep using the existing sample approximation. We introduce
auxiliary functions $f_{t,u}:{}^{d}\times{}^{|\mathcal{V}_{u}|}\to\real$ and
$g_{t,u}:{}^{|\mathcal{V}_{u}|}\times{}^{|\mathcal{V}_{u}|}\to\real$ for
$t\geq 1$ and $u\in\mathbb{T}$, such that
$f_{t,\mathfrak{R}}=f_{t},g_{t,\mathfrak{R}}=g_{t}$ and for
$u\in\mathbb{T}\setminus\mathfrak{R}$, $f_{t,u}$ and $g_{t,u}$ serve as
proxies for marginals of the transition density and observation likelihood,
respectively. These auxiliary functions are used to define our collection of
target densities $(\widetilde{\gamma}_{t,u})_{t\geq 1,u\in\mathbb{T}}$ over
${}^{|\mathcal{V}_{u}|}$:

$\displaystyle\widetilde{\gamma}_{t,u}(z_{t,u})$
$\displaystyle=g_{t,u}(z_{t,u},(y_{t}(i))_{i\in\mathcal{V}_{u}})\int
f_{t,u}(x_{t-1},z_{t,u})\widetilde{\gamma}_{t-1,\mathfrak{R}}(x_{t-1})\mathrm{d}x_{t-1}$
(6)

where $z_{t,u}=(x_{t}(i))_{i\in\mathcal{V}_{u}}$ are the components of $x_{t}$
associated with node $u$ and $x_{t-1}$ denotes the previous state of the
system. The requirement that
$g_{t,\mathfrak{R}}=g_{t},f_{t,\mathfrak{R}}=f_{t}$ ensures that at the root
we obtain the distribution in (2). The integral w.r.t.
$\widetilde{\gamma}_{t-1,\mathfrak{R}}$ in (6) cannot be computed
analytically, hence, as in MPFs, we use a sample approximation of this
integral, as described in the next section.

If the marginals of $f_{t},g_{t}$ are available, one could use them to define
$g_{t,u},f_{t,u}$, but this is not essential: these intermediate distributions
can be essentially arbitrary up to the absolute continuity required to justify
the importance sampling steps although, of course, the variance of the
estimator will be influenced by this choice. The issue of specifying these
distributions is closely related to that of choosing the sequence of
artificial targets in a standard SMC sampler when one is only interested in
the final target distribution (see Del Moral et al. (2006), where they note
that optimizing this sequence is a very difficult problem). See Kuntz et al.
(2021, Section 4.2) for a theoretical perspective in the divide-and-conquer
context, which suggests that the optimal choice of intermediate target
distribution is the appropriate marginal of the root target, hence the
suggestion to use the marginals of $f_{t}$ and $g_{t}$ where these are
available. In practice, when doing this exactly is not feasible, this
perspective suggests that we should seek to approximate these distributions
with distributions with comparable tail behaviour to keep all importance
weights well controlled. If there is a substantial mismatch between the
children of a node and itself then choosing a parametric path between the two,
adaptively specifying a tempering sequence along that path and using SMC
techniques to approximate each distribution in turn may mitigate some
difficulties (See Appendix B).

While we believe that the choice of the best $f_{t,u},g_{t,u}$ will be model
dependent, we found that, in some cases, some choices are preferable. As an
example, take $f_{t}$ to be a Gaussian distribution with mean $\mu$ and
covariance $\Sigma$. Then, a possible choice for $\gamma_{t,u}$ are Gaussian
distributions with mean $\mu_{u}$ equal to the restriction of $\mu$ to the
components corresponding to node $u$ and covariance $\Sigma_{u}$ obtained by
subsetting $\Sigma$ and discarding all components not in $u$. In this context,
we found that setting $\Sigma_{u}=\Sigma(\mathcal{V}_{u})^{-1}$, i.e.
selecting the components of $\Sigma$ corresponding to node $u$ and then
inverting this matrix, instead of $\Sigma_{u}=\Sigma^{-1}(\mathcal{V}_{u})$
has lower computational cost and leads to more diffuse distributions and thus
better behaved mixture weights (10).

###  3.1 The Algorithm

For each time $t$, having identified the space decomposition over $\mathbb{T}$
and the collection of distributions
$(\widetilde{\gamma}_{t,u})_{u\in\mathbb{T}}$, we can apply Algorithm 1 to the
root $\mathfrak{R}$ of $\mathbb{T}$. However, since the integral in (6) is not
analytically tractable, we replace $\widetilde{\gamma}_{t-1,\mathfrak{R}}$
with an approximation provided by the particle population at the root of the
tree corresponding to $t-1$, $\\{z_{t-1,\mathfrak{R}}^{n}\\}_{n=1}^{N}$, i.e.
its particle approximation obtained at the previous time step, as is normally
done in MPFs, and define

$\displaystyle\gamma_{t,u}(z_{t,u})$
$\displaystyle:=g_{t,u}(z_{t,u},(y_{t}(i))_{i\in\mathcal{V}_{u}})\frac{1}{N}\sum_{n=1}^{N}f_{t,u}(z_{t-1,\mathfrak{R}}^{n},z_{t,u}).$
(7)

Given $\\{z_{t-1,\mathfrak{R}}^{n}\\}_{n=1}^{N}$, at each leaf node of the
tree we sample one component of $x_{t}$ per node from
$N^{-1}\sum_{n=1}^{N}K_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)$, the importance
weights are then given by

$\displaystyle w_{t,u}(z_{t,u},x_{1:t-1,u})$
$\displaystyle=\frac{g_{t,u}(z_{t,u},(y_{t}(i))_{i\in\mathcal{V}_{u}})\sum_{n=1}^{N}f_{t,u}(z_{t-1,\mathfrak{R}}^{n},z_{t,u})}{\sum_{n=1}^{N}K_{t,u}(z_{t-1,\mathfrak{R}}^{n},z_{t,u})}.$
(8)

As in the MPF, if we choose $K_{t,u}=f_{t,u}$, (8) simplifies dramatically to
$w_{t,y}(z_{t,u},x_{1:t-1,u})=g_{t,u}(z_{t,u},(y_{t}(i))_{i\in\mathcal{V}_{u}})$,
considerably reducing the cost of evaluating the weights at the leaves.

For any non-leaf node $u$ we gather the particle populations
$\\{z_{t,\ell(u)}^{n}\\}_{n=1}^{N}$ and $\\{z_{t,r(u)}^{n}\\}_{n=1}^{N}$ on
its left and right child, respectively, and obtain an approximation of the
product measure
$\gamma_{t,\mathcal{C}_{u}}:=\gamma_{t,\ell(u)}\times\gamma_{t,r(u)}$ using
the weighted product form estimator (3)

$\displaystyle\gamma_{t,\mathcal{C}_{u}}^{N}=N^{-2}\sum_{n_{1}=1}^{N}\sum_{n_{2}=1}^{N}w_{t,\ell(u)}(z_{t,\ell(u)}^{n_{1}})w_{t,r(u)}(z_{t,r(u)}^{n_{2}})\delta_{(z_{t,\ell(u)}^{n_{1}},z_{t,r(u)}^{n_{2}})},$
(9)

or one of the lower cost alternatives discussed in Section 3.3.

As in the DaC-SMC setting, we reweight the particle approximation of
$\gamma_{t,\mathcal{C}_{u}}$ to target $\gamma_{t,u}$. In this case, the
mixture weights are given by

$\displaystyle m_{t,u}(z_{t,\mathcal{C}_{u}})=$
$\displaystyle\frac{\gamma_{t,u}(z_{t,\mathcal{C}_{u}})}{\gamma_{t,\mathcal{C}_{u}}(z_{t,\mathcal{C}_{u}})}$
(10) $\displaystyle=$
$\displaystyle\frac{g_{t,u}(z_{t,\mathcal{C}_{u}},(y_{t}(i))_{i\in\mathcal{V}_{u}})}{g_{t,\ell(u)}(z_{t,\ell(u)},(y_{t}(i))_{i\in\mathcal{V}_{\ell(u)}})g_{t,r(u)}(z_{t,r(u)},(y_{t}(i))_{i\in\mathcal{V}_{r(u)}})}\times$
$\displaystyle\frac{N^{-1}\sum_{n=1}^{N}f_{t,u}(z_{t-1,\mathfrak{R}}^{n},z_{t,\mathcal{C}_{u}})}{N^{-1}\sum_{n=1}^{N}f_{t,\ell(u)}(z_{t-1,\mathfrak{R}}^{n},z_{t,\ell(u)})N^{-1}\sum_{n=1}^{N}f_{t,r(u)}(z_{t-1,\mathfrak{R}}^{n},z_{t,r(u)})},$

where we defined $z_{t,\mathcal{C}_{u}}:=(z_{t,\ell(u)},z_{t,r(u)})$ the
vector obtained by merging the components on the left and on the right child
of $u$.

For each pair in (9) we obtain the incremental mixture weights
$m_{t,u}^{(n_{1},n_{2})}:=m_{t,u}((z_{t,\ell(u)}^{n_{1}},z_{t,r(u)}^{n_{2}}))$
in (10) and the updated weights

$\displaystyle\tilde{w}_{t,u}^{(n_{1},n_{2})}=\tilde{w}_{t,u}((z_{t,\ell(u)}^{n_{1}},z_{t,r(u)}^{n_{2}})):=w_{t,\ell(u)}(z_{t,\ell(u)}^{n_{1}})w_{t,r(u)}(z_{t,r(u)}^{n_{2}})m_{t,u}((z_{t,\ell(u)}^{n_{1}},z_{t,r(u)}^{n_{2}})),$

for $n_{1},n_{2}=1,\dots,N$. To avoid unbounded growth in the number of
particles, the weights $\tilde{w}_{t,u}$ are then used to resample a
population of $N$ particles approximating $\gamma_{t,u}$,
$\\{\tilde{z}_{t,u}^{n},w_{t,u}=1\\}_{n=1}^{N}$ — although one could allow the
number of particles retained to grow as the simulation approaches the root of
the tree to accommodate the growing dimension of the space. Algorithm 2
summarizes the procedure described above with the natural modification at
$t=1$, where it is not necessary to “marginalise” over previous states and
simple importance sampling can be used at the leaf nodes. Contrary to
Algorithm 1, we do not include an additional MCMC step in this statement of
the algorithm, but one could easily be added as discussed in Section 3.2. In
this case, given a $\pi_{t,u}$-invariant kernel $Q_{t,u}$, with
$\pi_{t,u}\propto\gamma_{t,u}$, line $8$ should be replaced by: draw
$z_{t,u}^{n}\sim Q_{t,u}(\tilde{z}_{t,u}^{n},\cdot)$ for all $n\leq N$.

1: for $u$ leaf node do

2: Initialize: draw $z_{t,u}^{n}\sim
N^{-1}\sum_{n=1}^{N}K_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)$ and compute
$w_{t,u}^{n}$ as in (8) for all $n\leq N$.

3: end for

4: for $u$ non-leaf node do

5: Recurse: set
$(\\{z_{t,v}^{n},w_{t,v}^{n}\\}_{n=1}^{N}):=\text{dac\\_smc}(t,v)$ for $v$ in
$\\{\ell(u),r(u)\\}$ and obtain $\gamma_{\mathcal{C}_{u}}^{N}$ in (9).

6: Merge: compute the mixture weights $m_{t,u}^{(n_{1},n_{2})}$ in (10) and
$\tilde{w}_{t,u}^{(n_{1},n_{2})}$ for all $n_{1},n_{2}\leq N$.

7: Resample: draw $\\{\tilde{z}_{t,u}^{n}\\}_{n=1}^{N}$ using weights
$\tilde{w}_{t,u}^{(n_{1},n_{2})}$ and set $w_{u}^{n}=1$ for all $n\leq N$.

8: Update: set $z_{t,u}^{n}=\tilde{z}_{t,u}^{n}$ for all $n\leq N$.

9: end for

10: Output $(\\{z_{t,\mathfrak{R}}^{n}\\}_{n=1}^{N})$.

Algorithm 2 dac_smc$(t)$ for $t\geq 1$. Given
$(\\{z_{t-1,\mathfrak{R}}^{n}\\}_{n=1}^{N}):=\text{dac\\_smc}(t-1)$.

The mixture resampling strategy described above requires evaluating the
mixture weights (10) for each of the $N^{2}$ pairs in (9). The $O(N^{2})$ cost
of this operation is often prohibitive for large $N$; we describe smaller cost
alternatives in Section 3.3 and demonstrate that these approaches have good
performance in Section 4.

###  3.2 Choice of Proposals

Algorithm 2 describes a general strategy to perform filtering using DaC-SMC;
as in the case of standard SMC the performances of the algorithm are heavily
influenced by the choice of the proposals at the leaf nodes. We discuss here a
simple strategy to select $K_{t,u}$.

We assume that we can sample from $f_{t,u}$ in (6) and set
$K_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)=f_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)$
so that the importance weights (8) reduce to $w_{t,u}=g_{t,u}$. This choice
corresponds to the proposal used in the bootstrap PF of Gordon et al. (1993),
(locally) optimal proposals also exist (see e.g., Chopin and Papaspiliopoulos
(2020, Chapter 10)) and are expected to lead to better performances but incur
a higher computational cost.

While picking
$K_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)=f_{t,u}(z_{t-1,\mathfrak{R}}^{n},\cdot)$
causes standard marginal PFs to reduce to the bootstrap PF (as described in
Klaas et al. (2005, Section 3.3)), this is not true for our marginal DaC-SMC,
since the integral w.r.t. $z_{t-1,\mathfrak{R}}$ still appears in the mixture
weights (10).

If needed, to avoid particle impoverishment, one might consider applying a
Markov kernel $Q_{t,u}$ which leaves $\pi_{t,u}\propto\gamma_{t,u}$ invariant
after the resampling step in line 7 of Algorithm 2. These
$\pi_{t,u}$-invariant kernels can be selected exploiting the vast literature
on sequential MCMC methods (e.g., Gilks and Berzuini (2001); Septier et al.
(2009); Carmi et al. (2012); Septier and Peters (2016); Pal and Coates (2018);
Han and Nakamura (2021)) to employ proposals whose cost is not $O(N)$ as it
would be for some naïve choices.

###  3.3 Adaptive Lightweight Mixture Resampling

The mixture resampling in line 6–7 of Algorithm 2 becomes computationally
impractical for large $N$, since it requires evaluating the mixture weights
for $N^{2}$ particles (Lindsten et al., 2017; Kuntz et al., 2021; Corneflos et
al., 2022). Several strategies have been proposed to alleviate this cost by
only constructing a subset of the $N^{2}$ combinations in (9), e.g. the
multiple matching strategy of Lin et al. (2005) which gives rise to the
lightweight mixture resampling of Lindsten et al. (2017) in this context,
strategies borrowed from the literature on incomplete U-statistics (Kuntz et
al., 2021) and lazy resampling schemes (Corneflos et al., 2022).

We consider the lightweight version of mixture resampling proposed in Lindsten
et al. (2017) which only considers a subset $\theta N$ with $\theta\ll N$ of
the $N^{2}$ possible pairs. However, instead of setting $\theta$ to some pre-
specified value (e.g. $\theta=\lceil\sqrt{N}\rceil$), we propose a simple
strategy to select $\theta$ adaptively based on the effective sample size
(ESS; Kong et al. (1994)),

$\displaystyle\textrm{ESS}:=\left(\sum_{n}\tilde{w}_{t,u}^{n}\right)^{2}/\sum_{n}(\tilde{w}_{t,u}^{n})^{2},$
(11)

where the sum is over all pairs $n=(n_{1},n_{2})$ obtained from (9), which is
similar in spirit to the adaptive tempering strategies commonly encountered in
the SMC literature (see, e.g., Jasra et al. (2010); Johansen (2015)) and
aiming to do just enough computation to obtain a good $N$-sample
approximation.

The merge step in lines 6–7 of Algorithm 2 is replaced by Algorithm 3: after
building all the $N$ pairs obtained by concatenating the particles associated
with each of the two children, further permutations are added until the ESS
achieves a pre-specified value $\textrm{ESS}^{\star}$ (e.g.
$\textrm{ESS}^{\star}=N$). To avoid $\theta$ getting too large, we stop adding
permutations when $\theta=\lceil\sqrt{N}\rceil$ thereby allowing us to bound
the worst-case computational cost by $O(N^{3/2})$. We empirically compare
several mixture resampling approaches in Appendix A.

1: Correct: compute the mixture weights
$m_{t,u}(x_{\ell(u)}^{n},x_{r(u)}^{n})$ as in (10) and
$\tilde{w}_{t,u}(x_{\ell(u)}^{n},x_{r(u)}^{n})$ for all $n\leq N$ and the ESS
(11).

2: Set: $\theta\leftarrow 1$ and
$\tilde{x}_{u}^{n}=(x_{\ell(u)}^{n},x_{r(u)}^{n})$ for $n\leq N$.

3: while $\textrm{ESS}((\tilde{x}^{n})_{n=1}^{\theta N})<\textrm{ESS}^{\star}$
and $\theta<\left\lceil\sqrt{N}\right\rceil$ do

4: Set: $\theta\leftarrow\theta+1$.

5: Permute: draw one permutation of $N$, $\pi(N)$, set
$\tilde{x}_{u}^{N(\theta-1)+n}=(x_{\ell(u)}^{n},x_{r(u)}^{\pi(n)})$, compute
the mixture weights $m_{t,u}(\tilde{x}_{u}^{N(\theta-1)+n})$ in (10) and the
updated weights $\tilde{w}_{t,u}(\tilde{x}_{u}^{N(\theta-1)+n})$ for $n\leq N$
and update the ESS.

6: end while

7: Resample: draw $\\{z_{t,u}^{n}\\}_{n=1}^{N}$ from
$\tilde{x}_{u}^{1:N\theta}$ with weights
$\tilde{w}_{t,u}(\tilde{x}_{u}^{1:N\theta})$ and set $w_{t,u}^{n}=1$ for all
$n\leq N$.

Algorithm 3 Adaptive lightweight mixture resampling.

Algorithm 3 can be implemented in a space-efficient manner by storing the
permutations $\pi(N)$ corresponding to each value of $\theta$ rather than
building the $\theta N$ pairs $\tilde{x}_{u}^{1:N\theta}$.

We report in Figure 2 the distribution of the number of permutation $\theta$
selected by Algorithm 3 for the linear Gaussian model in Section 4. We observe
that the highest values of $\theta$ are chosen at the level above the leaves
(panel $1$), indeed, at the leaf level the observation $y_{t}$ is
incorporated, causing a larger adjustment to the distribution at the first
mixture resampling step than that needed as we move up the tree. The spike at
$\theta=\lceil\sqrt{N}\rceil$ at level 1 shows that sometimes the target ESS
is not reached, which suggests that the product of the proposal distributions
over the children might be a poor proposal for $\gamma_{t,u}$ — this is in
part likely to be due to the use of a “bootstrap proposal” and could be
mitigated in the same way as standard particle filters by seeking to design
(marginal) proposal distributions which incorporate the influence of
observations. In addition, this phenomenon can be mitigated with the use of
(adaptive) tempering, as shown in, e.g., Jasra et al. (2010); Johansen (2015);
Wang et al. (2020); Zhou et al. (2016) for standard SMC and Lindsten et al.
(2017, Section 4.2) for DaC-SMC. However, naïve implementation of this
approach will bear a substantial computational cost in the marginal context;
designing MCMC kernels that are efficient in this context has been explored in
the sequential MCMC context — see, e.g., Septier and Peters (2016, Section
III-B) — and in marginal STPFs (Xu and Jasra, 2019). We give details of a
tempering strategy for Algorithm 3 in Appendix B.

Figure 2: Distribution of the number of permutations $\theta$ selected by
Algorithm 3 for the simple linear Gaussian model in Section 4.1 with $d=32$,
$N=10^{3}$ and 10 time steps; the 5 panels correspond to the 5 levels of the
tree from the level above the leaves (level 1) to the root (level 5).

###  3.4 Computational Cost

At a given node, $u$, the runtime cost of Algorithm 2 with lightweight mixture
resampling given the particle approximations for its children, is
$O(h(N,d)\theta N)$ where $h(N,d)$ denotes the cost of obtaining the mixture
weights (10).

In general, the dependence of $h(N,d)$ on the number of particles $N$ is
$O(N)$ as for MPFs. Similarly to MPFs, one could try to reduce the cost of
computing (10) making use of techniques from $N$-body learning (e.g. Gray and
Moore (2000); Lang et al. (2005)) as shown in Klaas et al. (2005).
Alternatively, one could consider efficient implementations using GPUs
(Charlier et al., 2021) as shown in (Clarté et al., 2019, Section 4) for sums
of the form of those in (10).

For some models, it might be possible to pick the auxiliary functions
$f_{t,u}$ so that the dependence on the past (i.e. $z_{t-1,\mathfrak{R}}$ in
(7)) vanishes and obtain an $O(1)$ cost w.r.t. $N$. However, we expect that
this type of decompositions will require larger corrections at the root, where
$f_{t,\mathfrak{R}}=f_{t}$, which might offset the cost savings.

In the adaptive case, worst case costs are given by $O(h(N,d)\theta N)$ with
$\theta$ replaced by the upper bound imposed upon the number of permutations
considered. For instance, for the examples in Section 4 we have an
$h(N,d)=O(N)$ cost to obtain the mixture weights and we set the upper bound
for $\theta$ to be $N^{1/2}$ leading to a cost of order $O(N^{5/2})$ w.r.t.
the number of particles.

Denoting $C_{u}(d,\theta_{u},N)$ the cost of running Algorithm 2 at node $u$,
we can then bound the total cost of serial implementations Algorithm 2 applied
at the root node $\mathfrak{R}$ by $O(dt\sup_{u}C_{u})$, where the supremum is
taken over all nodes in $\mathbb{T}$; a lower running cost of
$O(t\sup_{u}C_{u}\log_{2}d)$, can be achieved parallelizing the computations
over each level of the tree (Lindsten et al., 2017, Section 5.3).

In the adaptive case, this upper bound is far from being tight, since, as
shown in the histogram in Figure 2, $\theta$ tends to be high when the
observation is incorporated (level 1) but stabilizes as we move up the tree.
Additionally, for large $N$ one expects the number of permutations required to
obtain a good $N$-particle approximation to converge to some fixed integer and
hence the cost for sufficiently large $N$ will with high probability be of
smaller order than these bounds. Furthermore the constants multiplying the
$N^{5/2}$ contribution arising from the level above the leaves are
sufficiently small that this is not the dominant cost in our experiments — and
is likely to be typical in high-dimensional settings in which it is rarely
feasible to employ very large numbers of particles and the objective is to
obtain a good approximation at acceptable time _and space_ costs.

##  4 Experiments

We compare the results obtained with DaC with those of NSMC and STPF; we do
not include simpler strategies because both standard PF and BPF have been
shown to have worse performances than NSMC and STPF for the model considered
here (Næsseth et al., 2015, 2019; Beskos et al., 2017), and the marginal
version of STPF because of the higher cost for large $d$.

The functions $g_{t,u},f_{t,u}$ in (6) are obtained from $g_{t},f_{t}$,
respectively, by discarding all the terms in those functions involving
components $i\not\in\mathcal{V}_{u}$, further details are given in Appendix C.
For Algorithm 2 we use the proposals discussed in Section 3.2 and the
lightweight mixture resampling strategies described in Section 3.3. All
resampling steps are performed using stratified resampling (Kitagawa, 1996).

First, we consider a simple linear Gaussian SSM, and compare the results
obtained by the three algorithms with the exact filtering distribution given
by the Kalman filter. Then, we consider a spatial model with simple latent
dynamics but non-trivial spatial correlations structure between observations,
moving away from the assumption of i.i.d. observations which is convenient
from a computational perspective, but rarely satisfied in practice (Chib et
al., 2009).

All the experiments have been executed in serial using a single core of a
Intel(R) Xeon(R) CPU E5-2440 0 @ 2.40GHz using R 4.1.0.

###  4.1 Simple Linear Gaussian Model

We start by considering a simple linear Gaussian SSM taken from Næsseth et al.
(2015), for which the filtering distributions can be computed exactly with the
Kalman filter. The model is given by
$f_{t}(x_{t-1},x_{t})=\mathcal{N}(x_{t};Ax_{t-1},\Sigma)$, and
$g_{t}(x_{t},y_{t})=\mathcal{N}(y_{t};x_{t},\sigma^{2}_{y}\textrm{Id}_{d})$,
with $A\in{}^{d\times d}$, $\sigma_{y}^{2}>0$, $\Sigma\in{}^{d\times d}$ a
tridiagonal covariance matrix and $\textrm{Id}_{d}$ the $d$-dimensional
identity matrix (see Appendix C.1 for full details and computation of the
mixture weights (10)).

We compare DaC with both non-adaptive and adaptive lightweight mixture
resampling with 2-level NSMC with fully adapted outer level and STPF on data
simulated from the model for $d=2^{5},2^{8},2^{11}$ for $t=100$ time steps. We
use different number of particles $N=100,500,1000$ for Algorithm 2 and the
outer level of NSMC and STPF, while the number of particles for the inner
level of NSMC and the number of particles for each island of STPF is fixed to
$M=100$ as suggested in Næsseth et al. (2015); Beskos et al. (2017).

To evaluate the results, we consider two global measures of accuracy for each
of the $d$ marginals, the Wasserstein-1 distance (see, e.g., Vallender
(1974)), and the Kolmogorov-Smirnov distance

$\displaystyle W_{1,i}:=\int\left\lvert
F_{t,i}(x)-\widehat{F}_{t,i}(x)\right\rvert\mathrm{d}x,\qquad\qquad\textrm{KS}_{i}:=\max_{x}\left\lvert
F_{t,i}(x)-\widehat{F}_{t,i}(x)\right\rvert,$

where $F_{t,i}$ denotes the 1-dimensional cumulative distribution function of
marginal $i$ at time $t$ and $\widehat{F}_{t,i}$ its particle approximation.
Further comparisons, which demonstrate that the mean squared error
($\operatorname{MSE}$) of the filtering mean behaves similarly, are collected
in Appendix D.1.

$W_{1}$$d=32$$d=256$$d=2048$KSRuntime / sRuntime / s Figure 3: Comparison of
DaC, NSMC and STPF for $d=2^{5},2^{8},2^{11}$. Distribution of the average
(over dimension) $W_{1}$ and KS distance at the last time step $t=100$ for 50
runs; the boxes, form left to right, correspond to increasing number of
particles ($N=100,500,1000$). Due to their excessive cost, we do not include
the results for STPF with $d=2048$ and those of the non-adaptive version of
DaC $d=2048,N=1000$.

For lower-dimensional problems (e.g. $d=32$) STPF achieves the best results
both in terms of Wasserstein-1 distance and in terms of KS distance (Figure 3)
and the relative $\operatorname{MSE}$ of the reconstructions is considerably
smaller (almost one order of magnitude smaller; see Appendix D.1). The results
provided by STPF deteriorate quickly as $d$ grows, for $d=256$ the estimates
of $W_{1}$ and KS distance are significantly worse than those provided by NSMC
or DaC without adaptation.

The cost of STPF grows quadratically with $d$, and becomes unmanageable for
large $d$, it is therefore not included in Figure 3 bottom panels. STPF has
the higher cost also for lower dimension (Figure 3 top panel), but in this
case STPF provides the best results. STPF has higher variability than the
other methods, and even when the average results are better than DaC and NSMC
(e.g. $d=32$), $W_{1}$ and KS can take considerably high values.

The results in terms of KS are generally more variable, this is likely due to
the fact that KS is a measure of the worst case mismatch between $F_{t,i}(x)$
and $\widehat{F}_{t,i}(x)$, while for $W_{1}$ the mismatch is averaged over
locations. For large $d$, DaC has the smallest variability among the three
algorithms.

DaC with fixed-cost lightweight mixture resampling generally gives better
results than the adaptive lightweight mixture resampling, however, the cost of
the latter is considerably smaller, making the adaptive version still
manageable for large $N$ whereas the fixed-cost lightweight mixture resampling
becomes too costly for large $N$ and large $d$. As discussed in Section 3.4,
the computational cost of both versions of DaC could be reduced using GPUs. In
particular, both $W_{1}$ and KS decay more quickly with $N$ for DaC without
adaptation than for DaC with adaptive lightweight mixture resampling. The
decay with $N$ is less evident for NSMC.

###  4.2 Spatial Model

We consider a model on a 2D-lattice in which the latent dynamics are simple
but the observation structure is challenging. The components of $X_{t}$ are
indexed by the vertices $v\in V$ of a lattice, where $V=\\{1,\dots,d\\}^{2}$,
and follow a simple linear evolution $X_{t}(v)=X_{t-1}(v)+U_{t}(v)$, where
$U_{t}(v)\overset{\textrm{i.i.d}}{\sim}\mathcal{N}(0,\sigma_{x}^{2})$. The
observations model is $Y_{t}=X_{t}+V_{t}$, where we take $V_{t}$ to be jointly
$t$-distributed with $\nu=10$ degrees of freedom, mean zero and precision
structure encapsulating a spatial component. Let $D$ denote the graph
distance, then the entry in row $v$ and column $j$ of the precision matrix
$\Sigma^{-1}$ is given by $(\Sigma^{-1})_{vj}=\tau^{D(j,v)}$ if $D(j,v)\leq
r_{y}$ and 0 otherwise. We obtain data from the model above with
$\sigma^{2}_{x}=1$, $\tau=-0.25$, $r_{y}=1$ and $t=10$. The observation
density does not factorize, and therefore NSMC and STPF cannot be applied (at
least without approximating $g$ with e.g. a Gaussian or discarding the
covariance information). To validate the correctness of the algorithm, we
compare the results obtained by DaC with those of the standard bootstrap PF in
Appendix D on a small lattice and found the agreement to be excellent.

To decompose the 2D lattice into a binary tree we use the decomposition
described in Lindsten et al. (2017, Section 5.1), which recursively connects
the vertices first horizontally and then vertically. To evaluate the
performances of the algorithm we consider the filtering means obtained with 50
repetitions of DaC-SMC on a $8\times 8$ and a $16\times 16$ grid for
$N=100,500,1000$ and $5000$. To show how the standard bootstrap particle
filter struggles with higher dimensional problems we run a bootstrap PF with
$N=10^{5}$ particles for the $8\times 8$ grid. Figure 4 reports the filtering
means for a corner node and an interior node of the lattice. Both DaC
approaches are in agreement, however the adaptive version of DaC seems to
provide slightly less variable results. The behaviour for different nodes is
similar. As observed for the linear Gaussian model, DaC with adaptive
lightweight mixture resampling has lower cost than the non-adaptive
counterpart and remains feasible for large $N$ (e.g. $N=5000$).

Unsurprisingly, the bootstrap PF struggles to recover the filtering means and
provides high variance estimates for node $(1,1)$ while collapses completely
for node $(8,6)$ failing to recover the filtering mean.

$(1,1)$$(8,6)$$d=8\times 8$$(1,1)$$(8,8)$$d=16\times 16$Runtime / sRuntime / s
Figure 4: Filtering mean estimates for a corner node and a node in the middle
of the grid for a $8\times 8$ and a $16\times 16$ lattice at time $t=10$. The
boxplots from left to right report the distributions over 50 repetitions for
$N=100,500,1000$ and $5000$. The results for the non-adaptive version of DaC
are not included for $N=5000$ due to the excessive cost. The reference lines
for the $8\times 8$ grid show the average value of the filtering mean estimate
and the interquartile range obtained with 50 repetitions of a bootstrap PF
with $N=10^{5}$ particles.

The size of the boxplots in Figure 4 gives an indication of the variance of
the estimator provided by DaC, for small $N$ the decay in variance seems to be
more pronounced (at least for the adaptive version of DaC) than for large $N$,
this is consistent with the decay of the standard deviation which would be
expected from the variance expansions in Kuntz et al. (2021) where for small
$N$ the higher order contributions to the variance are not yet negligible (we
anticipate that a central limit theorem could be obtained by combining the
results of Kuntz et al. (2021) with those for marginal PFs).

##  5 Discussion

We introduced a novel sequential Monte Carlo algorithm, combining ideas from
marginal PFs and divide-and-conquer SMC to extend the latter to the filtering
context. This algorithm is based on a novel space decomposition for high
dimensional SSM which allows to recursively merge low dimensional marginals of
the filtering distribution to obtain the full filtering distribution, taking
into account the mismatch between product of marginals and joint distributions
using importance sampling. In principle, the DaC-SMC approach is amenable to
distributed implementation, although the marginalization technique employed
herein would necessitate significant communication from the node which
computes the overall filtering distribution at time $t-1$ to all nodes
involved in computing at time $t$ and we have not explored that direction
here.

In contrast with Nested SMC and space-time PFs, the DaC-SMC approach to
filtering can be applied when the marginals of the joint density (1) are not
available analytically. The computational cost of this new approach grows
polynomially with the number of particles $N$, however, this cost can be
reduced exploiting GPU routines to reduce the cost of computing the weights as
discussed in Section 3.4.

The experiments in Section 4.1 show that DaC-SMC achieves comparable
performances of NSMC and STPF with a runtime that remains competitive even for
large $d$ (but small $N$), contrary to STPFs. In addition, DaC-SMC can be
applied to filtering problems which do not allow for factorization as shown in
Section 4.2 and Appendix D. The variance decay in Figure 4 and Figure 6 in
Appendix D.1 suggest that this extended DaC-SMC achieves the same convergences
rates as DaC-SMC (Kuntz et al., 2021) for sufficiently large $N$ and we
anticipate that techniques used to analyze the marginal particle filter could
be combined with those in order to provide formal convergence results for the
method developed herein. The adaptive lightweight mixture resampling discussed
in Section 3.3 is a promising route to further reduce the computational cost
of DaC-SMC for filtering, however, as the experiments in Section 4 and
Appendix A suggest, selecting the value at which the target ESS should be set
to obtain the best trade-off between computational cost and accuracy is likely
to be problem dependent and raises interesting theoretical questions that we
leave for future work.

For challenging problems it is likely that tempering and MCMC kernels would be
required to give good performance. As discussed in Johansen (2015) and
Guarniero et al. (2017) the smoothing and filtering distributions (i.e.
$p(x_{1:t}|y_{1:t})$ and $p(x_{t}|y_{1:t})$, respectively) have significantly
different support in the presence of informative observations, especially in
high dimensional settings and so we expect that including the influence of
future observations in the proposals and targets in Section 3.2 (as in
lookahead methods, e.g., Lin et al. (2013); Guarniero et al. (2017); Ruzayqat
et al. (2022)) would lead to considerable improvements in the accuracy of the
estimates and might ultimately be essential in the development of good general
purpose filters for high dimensional problems.

This work focuses on obtaining approximations of the filtering distribution
for high dimensional SSM. In recent years there has been a lot of interest in
obtaining approximations of the smoothing distribution which is a necessary
component of parameter estimation algorithms (e.g., Finke and Singh (2017);
Guarniero et al. (2017)); we anticipate that the DaC-SMC approach to filtering
could be extended to tackle smoothing and parameter estimation dealing with
the marginalization in (6). In principle, algorithms which directly
approximate only marginals of smoothing distributions can be adapted to these
settings (see, e.g. Gerber and Chopin (2017)). We leave this for future work.

## Supplementary Materials

Supplementary materials contains details of the tempering approach, additional
details and results on the experiments. An R package reproducing the
experiments is available at
https://github.com/FrancescaCrucinio/Dac4filtering.

## Acknowledgements

FRC and AMJ acknowledge support from the EPSRC (grant # EP/R034710/1). AMJ
acknowledges further support from the EPSRC (grant # EP/T004134/1) and the
Lloyd’s Register Foundation Programme on Data-Centric Engineering at the Alan
Turing Institute. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

## References

  * Akyildìz et al. (2022) Akyildìz, D., D. Crisan, and J. Miguez (2022).  Space-sequential particle filters for high-dimensional dynamical systems described by stochastic differential equations.  arXiv preprint arXiv:2204.07680.
  * Bengtsson et al. (2008) Bengtsson, T., P. Bickel, and B. Li (2008).  Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems.  In Probability and Statistics: Essays in Honor of David A. Freedman, pp. 316–334. Institute of Mathematical Statistics.
  * Beskos et al. (2017) Beskos, A., D. Crisan, A. Jasra, K. Kamatani, and Y. Zhou (2017).  A stable particle filter for a class of high-dimensional state-space models.  Advances in Applied Probability 49(1), 24–48.
  * Carmi et al. (2012) Carmi, A., F. Septier, and S. J. Godsill (2012).  The Gaussian mixture MCMC particle algorithm for dynamic cluster tracking.  Automatica 48(10), 2454–2467.
  * Charlier et al. (2021) Charlier, B., J. Feydy, J. A. Glaunes, F.-D. Collin, and G. Durif (2021).  Kernel operations on the GPU, with autodiff, without memory overflows.  Journal of Machine Learning Research 22(74), 1–6.
  * Chib et al. (2009) Chib, S., Y. Omori, and M. Asai (2009).  Multivariate stochastic volatility.  In Handbook of Financial Time Series, pp. 365–400. Springer.
  * Chopin and Papaspiliopoulos (2020) Chopin, N. and O. Papaspiliopoulos (2020).  An Introduction to Sequential Monte Carlo.  Cham: Springer.
  * Clarté et al. (2019) Clarté, G., A. Diez, and J. Feydy (2019).  Collective proposal distributions for nonlinear MCMC samplers: Mean-field theory and fast implementation.  arXiv preprint arXiv:1909.08988.
  * Corneflos et al. (2022) Corneflos, A., N. Chopin, and S. Särkkä (2022).  De-Sequentialized Monte Carlo: a parallel-in-time particle smoother.  Journal of Machine Learning Research 23(283), 1–39.
  * Del Moral et al. (2006) Del Moral, P., A. Doucet, and A. Jasra (2006).  Sequential Monte Carlo samplers.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(3), 411–436.
  * Doucet and Johansen (2011) Doucet, A. and A. M. Johansen (2011).  A tutorial on particle filtering and smoothing: Fifteen years later.  In D. Crisan and B. Rozovsky (Eds.), The Oxford Handbook of Nonlinear Filtering, pp. 656–704. Oxford University Press.
  * Evensen (2009) Evensen, G. (2009).  Data Assimilation (2nd ed.).  Berlin, Heidelberg: Springer.
  * Finke and Singh (2017) Finke, A. and S. S. Singh (2017).  Approximate smoothing and parameter estimation in high-dimensional state-space models.  IEEE Transactions on Signal Processing 65(22), 5982–5994.
  * Gerber and Chopin (2017) Gerber, M. and N. Chopin (2017).  Convergence of sequential quasi-Monte Carlo smoothing algorithms.  Bernoulli 23(4B), 2951–2987.
  * Gerber et al. (2019) Gerber, M., N. Chopin, and N. Whiteley (2019).  Negative association, ordering and convergence of resampling methods.  Annals of Statistics 47(4), 2236–2260.
  * Gilks and Berzuini (2001) Gilks, W. R. and C. Berzuini (2001).  Following a moving target—Monte Carlo inference for dynamic Bayesian models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63(1), 127–146.
  * Gordon et al. (1993) Gordon, N. J., D. J. Salmond, and A. F. Smith (1993).  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  In IEE Proceedings F-Radar and Signal Processing, Volume 140, pp. 107–113. IEE.
  * Gray and Moore (2000) Gray, A. and A. Moore (2000).  ‘$N$-body’ problems in statistical learning.  In Advances in Neural Information Irocessing Systems, pp. 521–527.
  * Guarniero et al. (2017) Guarniero, P., A. M. Johansen, and A. Lee (2017).  The iterated auxiliary particle filter.  Journal of the American Statistical Association 112(520), 1636–1647.
  * Han and Nakamura (2021) Han, Y. and K. Nakamura (2021).  The application of Zig-Zag sampler in sequential Markov chain Monte Carlo.  arXiv preprint arXiv:2111.10210.
  * Jasra et al. (2010) Jasra, A., D. A. Stephens, A. Doucet, and T. Tsagaris (2010).  Inference for Lévy-driven stochastic volatility models via adaptive sequential Monte Carlo.  Scandinavian Journal of Statistics 38(1), 1–22.
  * Johansen (2015) Johansen, A. M. (2015).  On blocks, tempering and particle MCMC for systems identification.  IFAC-PapersOnLine 48(28), 969–974.
  * Kitagawa (1996) Kitagawa, G. (1996).  Monte Carlo filter and smoother for non-Gaussian nonlinear state space models.  Journal of Computational and Graphical Statistics 5(1), 1–25.
  * Klaas et al. (2006) Klaas, M., M. Briers, N. De Freitas, A. Doucet, S. Maskell, and D. Lang (2006).  Fast particle smoothing: If I had a million particles.  In Proceedings of the 23rd International Conference on Machine Learning, pp. 481–488.
  * Klaas et al. (2005) Klaas, M., N. De Freitas, and A. Doucet (2005).  Toward practical $N^{2}$ Monte Carlo: The marginal particle filter.  In Proceedings of the 21st Conference in Uncertainty in Artificial Intelligence, pp. 308–315.
  * Kong et al. (1994) Kong, A., J. S. Liu, and W. H. Wong (1994).  Sequential imputations and Bayesian missing data problems.  Journal of the American Statistical Association 89(425), 278–288.
  * Kuntz et al. (2021) Kuntz, J., F. R. Crucinio, and A. M. Johansen (2021).  The divide-and-conquer sequential Monte Carlo algorithm: theoretical properties and limit theorems.  arXiv preprint arXiv:2110.15782.
  * Kuntz et al. (2022) Kuntz, J., F. R. Crucinio, and A. M. Johansen (2022).  Product-form estimators: exploiting independence to scale up Monte Carlo.  Statistics and Computing 32(1), 1–22.
  * Lang et al. (2005) Lang, D., M. Klaas, and N. de Freitas (2005).  Empirical testing of fast kernel density estimation algorithms.  Technical Report TR2005-03, Department of Computer Science, University of British Columbia.
  * Lei et al. (2010) Lei, J., P. Bickel, and C. Snyder (2010).  Comparison of ensemble Kalman filters under non-Gaussianity.  Monthly Weather Review 138(4), 1293–1306.
  * Lin et al. (2013) Lin, M., R. Chen, J. S. Liu, et al. (2013).  Lookahead strategies for sequential Monte Carlo.  Statistical Science 28(1), 69–94.
  * Lin et al. (2005) Lin, M. T., J. L. Zhang, Q. Cheng, and R. Chen (2005).  Independent particle filters.  Journal of the American Statistical Association 100(472), 1412–1421.
  * Lindsten et al. (2017) Lindsten, F., A. M. Johansen, C. A. Næsseth, B. Kirkpatrick, T. B. Schön, J. A. D. Aston, and A. Bouchard-Côté (2017).  Divide-and-Conquer with sequential Monte Carlo.  Journal of Computational and Graphical Statistics 26(2), 445–458.
  * Liu (2001) Liu, J. S. (2001).  Monte Carlo Strategies in Scientific Computing.  New York: Springer.
  * Næsseth et al. (2015) Næsseth, C. A., F. Lindsten, and T. B. Schön (2015).  Nested sequential Monte Carlo methods.  In Proceedings of the 32nd International Conference on Machine Learning, Volume 37, pp. 1292–1301. Proceedings of Machine Learning Research.
  * Næsseth et al. (2019) Næsseth, C. A., F. Lindsten, and T. B. Schön (2019).  High-dimensional filtering using nested sequential Monte Carlo.  IEEE Transactions on Signal Processing 67(16), 4177–4188.
  * Pal and Coates (2018) Pal, S. and M. Coates (2018).  Sequential MCMC with the discrete bouncy particle sampler.  In IEEE Statistical Signal Processing Workshop (SSP), pp. 663–667. IEEE.
  * Pitt and Shephard (1999) Pitt, M. K. and N. Shephard (1999).  Filtering via simulation: Auxiliary particle filters.  Journal of the American Statistical Association 94(446), 590–599.
  * Rebeschini and Van Handel (2015) Rebeschini, P. and R. Van Handel (2015).  Can local particle filters beat the curse of dimensionality?  Annals of Applied Probability 25(5), 2809–2866.
  * Ruzayqat et al. (2022) Ruzayqat, H., A. Er-Raiy, A. Beskos, D. Crisan, A. Jasra, and N. Kantas (2022).  A lagged particle filter for stable filtering of certain high-dimensional state-space models.  SIAM/ASA Journal on Uncertainty Quantification 10(3).
  * Septier et al. (2009) Septier, F., S. K. Pang, A. Carmi, and S. Godsill (2009).  On MCMC-based particle methods for Bayesian filtering: Application to multitarget tracking.  In 3rd IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 360–363. IEEE.
  * Septier and Peters (2016) Septier, F. and G. W. Peters (2016).  Langevin and Hamiltonian based sequential MCMC for efficient Bayesian filtering in high-dimensional spaces.  IEEE Journal of Selected Topics in Signal Processing 10(2), 312–327.
  * Vallender (1974) Vallender, S. S. (1974).  Calculation of the Wasserstein distance between probability distributions on the line.  Theory of Probability & Its Applications 18(4), 784–786.
  * Wang et al. (2020) Wang, L., S. Wang, and A. Bouchard-Côté (2020).  An annealed sequential Monte Carlo method for Bayesian phylogenetics.  Systematic Biology 69(1), 155–183.
  * Xu and Jasra (2019) Xu, Y. and A. Jasra (2019).  Particle filters for inference of high-dimensional multivariate stochastic volatility models with cross-leverage effects.  Foundations of Data Science 1(1), 61–85.
  * Zhou et al. (2016) Zhou, Y., A. M. Johansen, and J. A. Aston (2016).  Toward automatic model comparison: an adaptive sequential Monte Carlo approach.  Journal of Computational and Graphical Statistics 25(3), 701–726.

##  Appendix A Comparison of Resampling Strategies

The divide and conquer SMC algorithm described in Section 2.3 merges the
particle populations of each node’s children using mixture resampling. We
compare the performances of lightweight mixture resampling with its adaptive
version described in Section 3.3, those of full-cost mixture resampling, in
which all possible $N^{2}$ permutations of the particles are created, and
those of a linear cost version of DaC-SMC in which no mixing weights are used
(Lindsten et al., 2017).

To compare the four resampling schemes we consider the mean squared error
($\operatorname{MSE}$) for each component $i$ using the example of Section
4.1,

$\displaystyle\operatorname{MSE}(x_{t}(i)):=\mathbb{E}\left[(\bar{x}_{t}(i)-\mu_{t,i})^{2}\right],$
(12)

where $\bar{x}_{t}(i)$ denotes the estimate of the mean of component $i$ at
time $t$ and $\mu_{t,i}$ denotes the true mean of $x_{t}(i)|y_{1:t}$ obtained
from the Kalman filter. To approximate (12) we consider an empirical average
over 50 repetitions.

Figure 5 shows the average (over component) $\operatorname{MSE}$ of the
estimates obtained for $d=128$ and $t=10$ time steps as a function of runtime.
The linear cost version of DaC-SMC has the smallest runtime, however, the
results in terms of $\operatorname{MSE}$ are the worst among the four
algorithms. In fact, this linear cost version does not use mixture weights at
the point of selection, and therefore does not take into account the mismatch
between $\gamma_{t,\mathcal{C}_{u}}$ and $\gamma_{t,u}$ when resampling (this
is corrected for with a subsequent importance reweighting). The full mixture
resampling has a higher cost, and becomes unmanageable for large $N$ (it is,
in fact, not included for $N=1000,5000$); the non-adaptive version of
lightweight mixture resampling also becomes too expensive for large $N$
although the increase in cost is less steep than that of the full mixture
resampling (the non-adaptive lightweight mixture resampling has manageable
cost for $N=1000$ but it is not included for $N=5000$). The
$\operatorname{MSE}$ achieved by the non-adaptive lightweight mixture
resampling is equivalent to that of the full mixture resampling, at a
considerably lower cost. For small $N$, the adaptive lightweight mixture
resampling also gives comparable results, but, as $N$ increases, the
$\operatorname{MSE}$ stops improving and eventually settles around $0.02$.
This is likely due to the fact that, as it is the case of Figure 2, sometimes
the target ESS is not reached at level 1 and more permutation would be needed
to obtain a good importance sampling proposal. One could then consider to add
tempering steps to obtain better samples or, alternatively, to increase the
target ESS. By how much the target ESS should be increased is a challenging
question which have not investigated this further in this work.

Runtime / s Figure 5: Average (over dimensions) $\operatorname{MSE}$ at $t=10$
over 50 runs for $d=128$ as function of runtime for the linear cost version of
DaC in Lindsten et al. (2017) and Algorithm 2 with lightweight mixture
resampling, adaptive lightweight mixture resampling and full mixture
resampling. The boxes, from left to right, correspond to increasing number of
particles, $N=100,500,1000,5000$. Due to the excessive cost, the results for
mixture resampling are not reported for $N=1000,5000$ and those for
lightweight mixture resampling with no adaptation are not reported for
$N=5000$.

##  Appendix B Tempering Strategy

If the product of the marginals over the two child nodes merged via mixture
resampling provides a poor approximation for $\gamma_{t,u}$ (which in some
circumstances could be detected, for example, by Algorithm 3 failing achieve
the target ESS), one could expect to alleviate this mismatch via a tempering
strategy. Here, we describe an adaptive tempering strategy of the sort
described in Lindsten et al. (2017, Section 4.2) adapted to our context.

We define the following sequence of targets
$\hat{\pi}_{\alpha_{u,j}}=\hat{\gamma}_{\alpha_{u,j}}/\hat{\gamma}_{\alpha_{u,j}}(1)$
where

$\displaystyle\hat{\gamma}_{\alpha_{u,j}}(z_{t,u,j})=\left[\gamma_{t,\mathcal{C}_{u}}(z_{t,u,j})\right]^{1-\alpha_{u,j}}\gamma_{t,u}(z_{t,u,j})^{\alpha_{u,j}}$

for $z_{t,u,j}\in{}^{|\mathcal{V}_{u}|}$ and $\alpha_{u,j}\in[0,1]$, which
interpolates from the proposal $\gamma_{t,\mathcal{C}_{u}}$ to the target at
node $u$, $\gamma_{t,u}$. In a standard approach to tempering we would start
from $\alpha_{u,j}=0$ and select the next value adaptively (see e.g. Wang et
al. (2020); Zhou et al. (2016)). However, Algorithm 3 corresponds to a first
tempering step which moves the proposal $\gamma_{t,\mathcal{C}_{u}}$ closer to
the target $\gamma_{t,u}$; therefore, instead of starting the tempering
schedule at $\alpha_{u,j}=0$, we identify the value of $\alpha^{\star}$ which
corresponds to the intermediate tempered target obtained after mixture
resampling and carry on with tempering from there until $\alpha_{u,j}=1$. To
identify the value of $\alpha$ we consider the ESS

$\displaystyle\textrm{ESS}(\alpha)$
$\displaystyle:=\left(\sum_{n=1}^{N}(\tilde{w}_{t,u}^{n})^{\alpha}\right)^{2}/\sum_{n=1}^{N}(\tilde{w}_{t,u}^{n})^{2\alpha},$

and solve $\textrm{ESS}(\alpha)=\textrm{ESS}^{\star}$, where
$\textrm{ESS}^{\star}$ is the target ESS in Section 3.3. The equally weighted
particles after mixture resampling,
$\\{z_{t,\mathcal{C}_{u}},\omega_{u,j}=1\\}_{n=1}^{N}$, approximate
$\hat{\pi}_{\alpha^{\star}}$. The following values of the tempering sequence
are chosen adaptively as described in, e.g., Zhou et al. (2016, Section 3.3.2
and Algorithm 4) and Wang et al. (2020, Algorithm 4), using the conditional
ESS

$\displaystyle\textrm{CESS}_{j}(\alpha)$
$\displaystyle:=N\left(\sum_{n=1}^{N}\bar{\omega}_{u,j}^{n}(\tilde{w}_{t,u}^{n})^{\alpha-\alpha_{u,j-1}}\right)^{2}/\sum_{n=1}^{N}\bar{\omega}_{u,j}^{n}(\tilde{w}_{t,u}^{n})^{2(\alpha-\alpha_{u,j-1})},$
(13)

where $\bar{\omega}_{u,j}^{n}\propto\omega_{u,j}^{n}$. Given a decay threshold
$\beta>0$, we find the next value $\alpha_{u,j}$ solving
$\textrm{CESS}_{j}(\alpha)=\beta$, then we update the weights, perform a
resampling step if necessary, and rejuvenate the particles using a
$\hat{\pi}_{\alpha_{u,j}}$-invariant kernel $K^{\prime}_{\alpha_{u,j}}$. This
process is repeated until $\alpha_{u,j}=1$. The resulting tempering strategy
is summarized in Algorithm 4.

1: Initialize: set $\alpha_{u,1}$ to be the solution of
$\textrm{ESS}(\alpha)=\textrm{ESS}^{\star}$ and set $\omega_{u,1}^{n}=1$ for
$n\leq N$.

2: Set: $j=1$.

3: while $\alpha_{u,j}<1$ do

4: Adapt: set $j=j+1$ and find $\alpha_{u,j}$ solving
$\textrm{CESS}_{j}(\alpha)=\beta$ for $\alpha\in(\alpha_{u,j-1},1)$ using
bisection.

5: Reweight: compute the weights
$\omega_{u,j}^{n}=\omega_{u,j-1}^{n,N}(\tilde{w}_{t,u}^{n})^{\alpha_{u,j}-\alpha_{u,j-1}}$
for $n\leq N$ and compute the ESS.

6: if $\textrm{ESS}<N/2$ then

7: Resample: draw new $z_{t,u,j-1}^{n,N}$ independently with weights
$\omega_{u,j}^{n}$ and update $\omega_{u,j}^{n}=1$ for $n\leq N$.

8: end if

9: Propose: draw $z_{t,u,j}^{n}\sim
K^{\prime}_{\alpha_{u,j}}(\cdot,z_{t,u,j-1}^{n})$ from a
$\hat{\pi}_{\alpha_{u,j}}$ invariant kernel for $n\leq N$.

10: end while

Algorithm 4 Tempering.

##  Appendix C Further Details on the Experiments

We collect here further details on the models considered in Section 4.

###  C.1 Simple Linear Gaussian Model

We consider the simple linear Gaussian model used in the experiments in
Næsseth et al. (2015, Section 6.1):

$\displaystyle X_{1}$ $\displaystyle\sim\mathcal{N}_{d}(0,\textrm{Id}_{d})$
$\displaystyle X_{t}$ $\displaystyle=0.5AX_{t-1}+U_{t},\qquad
U_{t}\sim\mathcal{N}_{d}(0,\Sigma)$ $\displaystyle Y_{t}$
$\displaystyle=X_{t}+V_{t},\qquad
V_{t}\sim\mathcal{N}_{d}(0,\sigma_{y}^{2}\textrm{Id}_{d}),$

with $A=A_{1}A_{2}^{-1}$ where

$\displaystyle A_{1}=\begin{pmatrix}\tau+\lambda&0&0&\dots&\dots&0&0\\\
0&\tau&0&0&\dots&0&0\\\ 0&\ddots&\ddots&\ddots&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\ 0&0&0&0&0&\tau&0\\\
0&0&0&0&0&0&\tau\end{pmatrix},$ $\displaystyle
A_{2}=\begin{pmatrix}\tau+\lambda&0&0&\dots&\dots&0&0\\\
-\lambda&\tau+\lambda&0&0&\dots&0&0\\\ 0&\ddots&\ddots&\ddots&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
0&0&0&0&-\lambda&\tau+\lambda&0\\\
0&0&0&0&0&-\lambda&\tau+\lambda\end{pmatrix}^{-1},$

and

$\displaystyle\Sigma^{-1}=A_{2}^{T}\begin{pmatrix}\tau&0&0&\dots&\dots&0&0\\\
0&\tau+\lambda&0&0&\dots&0&0\\\ 0&\ddots&\ddots&\ddots&\ddots&0&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
0&0&0&0&0&\tau+\lambda&0\\\ 0&0&0&0&0&0&\tau+\lambda\end{pmatrix}A_{2},$

with $\tau=\lambda=1$ and $\sigma^{2}_{y}=0.5^{2}$.

To derive the mixture weights for our divide-and-conquer approach we write the
joint density of $(X_{1:t},Y_{1:t})$:

$\displaystyle p(x_{1:t},y_{1:t})$
$\displaystyle\propto\prod_{k=1}^{t}\left[\exp\left(-\frac{1}{2}(x_{k}-0.5Ax_{k-1})^{T}\Sigma^{-1}(x_{k}-0.5Ax_{k-1})\right)\right.$
$\displaystyle\qquad\times\left.\prod_{i=1}^{d}\mathcal{N}(y_{k}(i);x_{k}(i),\sigma_{y}^{2})\right].$

The functions $f_{t,u},g_{t,u}$ in the definition of $\gamma_{t,u}$ in (6) are

$\displaystyle f_{t,u}(x_{t-1},x_{t}(\mathcal{V}_{u}))\propto$
$\displaystyle\qquad\exp\left(-\frac{1}{2}(x_{t}(\mathcal{V}_{u})-0.5A_{\mid\mathcal{V}_{u}}x_{t-1}(\mathcal{V}_{u}))^{T}\Sigma_{\mid\mathcal{V}_{u}}^{-1}(x_{t}(\mathcal{V}_{u})-0.5A_{\mid\mathcal{V}_{u}}x_{t-1}(\mathcal{V}_{u}))\right),$
$\displaystyle
g_{t,u}(x_{t}(\mathcal{V}_{u}),(y_{t}(i))_{i\in\mathcal{V}_{u}})\propto\prod_{i\in\mathcal{V}_{u}}\mathcal{N}(y_{t}(i);x_{t}(i),\sigma_{y}^{2}),$

where we denote $x_{t}(\mathcal{V}_{u})=(x_{t}(i))_{i\in\mathcal{V}_{u}}$ and
$A_{\mid\mathcal{V}_{u}},\Sigma_{\mid\mathcal{V}_{u}}^{-1}$ are the
restrictions of $A,\Sigma^{-1}$, respectively, to $\mathcal{V}_{u}$, i.e.
$A_{\mid\mathcal{V}_{u}}$ is the matrix obtained discarding the elements of
$A$ corresponding to components which are not in $\mathcal{V}_{u}$, and
similarly for $\Sigma_{\mid\mathcal{V}_{u}}^{-1}$.

Expanding the quadratic form in the display above, we find that we can
decouple the dependence on the current time step from that on the past and
decompose

$\displaystyle f_{t,u}(x_{t-1},(x_{t}(i))_{i\in\mathcal{V}_{u}})$
$\displaystyle=f_{t,u}^{(1)}(x_{t-1},(x_{t}(i))_{i\in\mathcal{V}_{u}})\prod_{i\in\mathcal{V}_{u},i\neq
j_{u}^{1}}\tilde{f}(x_{t}(i-1),x_{t}(i)),$

with

$\displaystyle f_{t,u}^{(1)}(x_{t-1},(x_{t}(i))_{i\in\mathcal{V}_{u}})$
$\displaystyle\propto\prod_{i\in\mathcal{V}_{u},i\neq
j_{u}^{1}}\left[\exp\left(-\frac{\tau+\lambda}{2}\left(x_{t}(i)-0.5\frac{\tau}{\tau+\lambda}x_{t-1}(i)\right)^{2}\right)\right.$
$\displaystyle\qquad\left.\times\exp\left(-\frac{1}{2}\frac{\lambda\tau}{\tau+\lambda}x_{t}(i-1)x_{t-1}(i)\right)\right]$
$\displaystyle\times\begin{cases}\exp\left(-\frac{\tau+\lambda}{2}\left(x_{t}(j_{u}^{1})-0.5\frac{\tau}{\tau+\lambda}x_{t-1}(j_{u}^{1})\right)^{2}\right)\qquad\textrm{if
}j_{u}^{1}\neq 1\\\
\exp\left(-\frac{\tau}{2}(x_{t}(1)-0.5x_{t-1}(1))^{2}\right)\qquad\textrm{if
}j_{u}^{1}=1\end{cases},$

where $j_{u}^{1}$ denotes the first index associated with node $u$, and

$\displaystyle\tilde{f}(s_{1},s_{2}):=\exp\left(-\frac{1}{2}\left[\frac{\lambda^{2}}{\tau+\lambda}s_{1}^{2}-2\lambda
s_{1}s_{2}\right]\right).$

The mixture weights are then given by

$\displaystyle m_{t,u}(z_{t,\mathcal{C}_{u}})$
$\displaystyle\propto\tilde{f}\left(z_{t,\ell(u)}(i_{\ell(u)}^{n_{\ell(u)}}),z_{t,r(u)}(i_{r(u)}^{1})\right)$
$\displaystyle\times\frac{\sum_{n=1}^{N}f_{t,u}^{(1)}(z_{t-1,\mathfrak{R}}^{n},z_{t,\mathcal{C}_{u}})}{\sum_{n=1}^{N}f_{t,\ell(u)}^{(1)}(z_{t-1,\mathfrak{R}}^{n},z_{t,\ell(u)})\sum_{n=1}^{N}f_{t,r(u)}^{(1)}(z_{t-1,\mathfrak{R}}^{n},z_{t,r(u)})},$

where we recall that $z_{t,u}$ in (6) corresponds to
$z_{t,u}=(x_{t}(i))_{i\in\mathcal{V}_{u}}$, $i_{\ell(u)}^{n_{\ell(u)}}$ is the
last index associated with node $\ell(u)$ and $i_{r(u)}^{1}$ is the first
index associated with $r(u)$.

###  C.2 Spatial Model

The joint density of $(X_{1:t},Y_{1:t})$ is

$\displaystyle p(x_{1:t},y_{1:t})$
$\displaystyle\propto\prod_{k=1}^{t}\left(\prod_{v\in
V}\mathcal{N}\left(x_{k}(v);x_{k-1}(v),\sigma_{x}^{2}\right)\right.$
$\displaystyle\left.\times\left[1+\nu^{-1}\sum_{v\in
V}\left((y_{k}(v)-x_{k}(v))\sum_{j:D(v,j)\leq
r_{y}}\tau^{D(v,j)}(y_{k}(j)-x_{k}(j))\right)\right]^{-(\nu+|V|)/2}\right)$

with initial distribution $X_{1}\sim\prod_{v\in
V}\mathcal{N}\left(x_{k}(v);0,\sigma_{x}^{2}\right)$. The functions
$f_{t,u},g_{t,u}$ in the definition of $\gamma_{t,u}$ in (6) are

$\displaystyle
f_{t,u}(x_{t-1},(x_{t}(i))_{i\in\mathcal{V}_{u}})\propto\prod_{v\in\mathcal{V}_{u}}\mathcal{N}\left(x_{k}(v);x_{k-1}(v),\sigma_{x}^{2}\right)$
$\displaystyle
g_{t,u}((x_{t}(i))_{i\in\mathcal{V}_{u}},(y_{t}(i))_{i\in\mathcal{V}_{u}})\propto$
$\displaystyle\qquad\left[1+\nu^{-1}\sum_{v\in\mathcal{V}_{u}}\left((y_{k}(v)-x_{k}(v))\sum_{j:D(v,j)\leq
r_{y},j\in\mathcal{V}_{u}}\tau^{D(v,j)}(y_{k}(j)-x_{k}(j))\right)\right]^{-(\nu+|\mathcal{V}_{u}|)/2},$

where we recall that $z_{t,u}$ in (6) corresponds to
$z_{t,u}=(x_{t}(i))_{i\in\mathcal{V}_{u}}$. From the above, we obtain that the
mixture weights (10) are given by
$m_{t,u}(z_{t,\mathcal{C}_{u}})=R^{f}_{t,u}(z_{t,\mathcal{C}_{u}})R^{g}_{t,u}(z_{t,\mathcal{C}_{u}})$,
where

$\displaystyle R^{f}_{t,u}(z_{t,\mathcal{C}_{u}})$
$\displaystyle=\frac{N^{-1}\sum_{n=1}^{N}f_{t,u}(z_{t-1,\mathfrak{R}}^{n},z_{t,\mathcal{C}_{u}})}{N^{-1}\sum_{n=1}^{N}f_{t,\ell(u)}(z_{t-1,\mathfrak{R}}^{n},z_{t,\ell(u)})N^{-1}\sum_{n=1}^{N}f_{t,r(u)}(z_{t-1,\mathfrak{R}}^{n},z_{t,r(u)})}$
$\displaystyle\propto\frac{\sum_{n=1}^{N}\prod_{v\in\mathcal{V}_{u}}\mathcal{N}\left(z_{t,\mathcal{C}_{u}}(v);z_{t-1,\mathfrak{R}}^{n}(v),\sigma_{x}^{2}\right)}{\sum_{n_{1}=1}^{N}\prod_{v\in\mathcal{V}_{\ell(u)}}\mathcal{N}\left(z_{t,\ell(u)}(v);z_{t-1,\mathfrak{R}}^{n_{1}}(v),\sigma_{x}^{2}\right)\sum_{n_{2}=1}^{N}\prod_{v\in
r(u)}\mathcal{N}\left(z_{t,r(u)}(v);z_{t-1,\mathfrak{R}}^{n_{2}}(v),\sigma_{x}^{2}\right)}$

and

$\displaystyle R^{g}_{t,u}(z_{t,\mathcal{C}_{u}})$
$\displaystyle=\frac{g_{t,u}(z_{t,\mathcal{C}_{u}},(y_{t}(i))_{i\in\mathcal{V}_{u}})}{g_{t,\ell(u)}(z_{t,\ell(u)},(y_{t}(i))_{i\in\mathcal{V}_{\ell(u)}})g_{t,r(u)}(z_{t,r(u)},(y_{t}(i))_{i\in\mathcal{V}_{r(u)}})}$
$\displaystyle\propto\left[1+\nu^{-1}\sum_{v\in\mathcal{V}_{u}}\left((y_{t}(v)-z_{t,\mathcal{C}_{u}}(v))\sum_{\begin{subarray}{c}j:D(v,j)\leq
r_{y}\\\
j\in\mathcal{V}_{u}\end{subarray}}\tau^{D(v,j)}(y_{t}(j)-z_{t,\mathcal{C}_{u}}(j))\right)\right]^{-\frac{\nu+|\mathcal{V}_{u}|}{2}}$
$\displaystyle\times\left[1+\nu^{-1}\sum_{v\in\mathcal{V}_{\ell(u)}}\left((y_{t}(v)-z_{t,\ell(u)}(v))\sum_{\begin{subarray}{c}j:D(v,j)\leq
r_{y}\\\
j\in\mathcal{V}_{\ell(u)}\end{subarray}}\tau^{D(v,j)}(y_{t}(j)-z_{t,\ell(u)}(j))\right)\right]^{\frac{\nu+|\mathcal{V}_{\ell(u)}|}{2}}$
$\displaystyle\times\left[1+\nu^{-1}\sum_{v\in\mathcal{V}_{r(u)}}\left((y_{t}(v)-z_{t,r(u)}(v))\sum_{\begin{subarray}{c}j:D(v,j)\leq
r_{y}\\\
j\in\mathcal{V}_{r(u)}\end{subarray}}\tau^{D(v,j)}(y_{t}(j)-z_{t,r(u)}(j))\right)\right]^{\frac{\nu+|\mathcal{V}_{r(u)}|}{2}}.$

##  Appendix D Additional Results for the Experiments

###  D.1 Simple Linear Gaussian Model

To further characterize the quality of the estimates obtained with Algorithm
2, we consider the relative mean squared error ($\operatorname{RMSE}$) for
component $i$ at time $t$,

$\displaystyle\operatorname{RMSE}(x_{t}(i)):=\frac{\mathbb{E}\left[(\bar{x}_{t}(i)-\mu_{t,i})^{2}\right]}{\sigma^{2}_{t,i}},$

where where $\bar{x}_{t}(i)$ denotes the estimate of the mean of component $i$
at time $t$ and $\mu_{t,i},\sigma^{2}_{t,i}$ denote the true mean and true
variance of $x_{t}(i)|y_{1:t}$ obtained from the Kalman filter. As for (12),
we approximate the $\operatorname{RMSE}$ using an empirical average over 50
repetitions.

$N=100$$N=500$$N=1000$$d=32$$d=256$$d=2048$Time step Figure 6: Average (over
dimension) $\operatorname{RMSE}$ for 50 runs of DaC, NSMC and STPF. Due to
their excessive cost, we do not include the results for STPF with $d=2048$ and
those of the non-adaptive version of DaC for $d=2048,N=1000$.

As observed in Section 4.1, for small $d$, STPF performs better than DaC and
NSMC (Figure 6 top row), however, as dimension increases, the estimates
provided by STPF become poor and their computational cost is prohibitive for
large $d$ (Figure 6 bottom row). Contrary to the other approaches, NSMC does
not seem to significantly improve as $N$ increases.

###  D.2 Spatial Model

To validate the results for the spatial model in Section 4.2 we compare the
approximations obtained by Algorithm 2 with both adaptive and non-adaptive
lightweight mixture resampling with those of the standard bootstrap particle
filter (see, e.g., Chopin and Papaspiliopoulos (2020, Chapter 10)) on a
$2\times 2$ grid (Figure 7). Given the low dimensionality of the state space,
we expect the bootstrap PF with a large number of particles to provide a good
proxy for the filtering distribution $p(x_{t}|y_{1:t})$. The results for the
three algorithms seem to be in good agreement, and the variance of both DaC
algorithms with large $N=5\cdot 10^{3},10^{4}$ is comparable with that of the
particle filter with $N=10^{5}$.

column $1$column $2$row $1$row $2$$N$$N$ Figure 7: Filtering mean estimates
obtained with DaC on a $2\times 2$ grid. The reference lines show the average
value of the filtering mean estimate and the interquartile range obtained with
50 repetitions of a bootstrap PF with $N=10^{5}$ particles. The boxplots from
left to right report the distributions over 50 repetitions for
$N=100,500,1000,5000$ and $10000$. The results for the non-adaptive version of
DaC and $N=10000$ are not included due to the excessive cost.

To show how the increasing dimension causes instability of the simple
bootstrap PF we report in Figure 8 the filtering mean estimates for two nodes
of a $4\times 4$ grid. Contrary to Figure 7, in which the bootstrap PF with
$N=10^{5}$ particles provides accurate estimates (first quartile, mean and
third quartile coincide), when the dimension slightly increases the variance
of the particle filter with $N=10^{5}$ blows up and the recovered estimates
are poor. On the other hand, both DaC approaches are stable and provide better
and better estimates for values of $N$ which are at least 10 times smaller
than $N=10^{5}$.

$(1,1)$$(2,3)$ Figure 8: Filtering mean estimates obtained with DaC for node
$(1,1)$ and $(2,3)$ of a $4\times 4$ grid. The reference lines show the
average value of the filtering mean estimate and the interquartile range
obtained with 50 repetitions of a bootstrap PF with $N=10^{5}$ particles. The
boxplots from left to right report the distributions over 50 repetitions for
$N=100,500,1000,5000$ and $10000$. The results for the non-adaptive version of
DaC and $N=10000$ are not included due to the excessive cost.