Medical Imaging with Deep LearningMedical Imaging with Deep Learning 2019
Extended Abstract – MIDL 2019 Submission Raphael Meier1<EMAIL_ADDRESS>
Michael Rebsamen1<EMAIL_ADDRESS>
Urspeter Knecht2<EMAIL_ADDRESS>
Mauricio Reyes2,3<EMAIL_ADDRESS>
Roland Wiest1<EMAIL_ADDRESS>
Richard McKinley1<EMAIL_ADDRESS>
1 Support Center for Advanced Neuroimaging (SCAN), University Institute of
Diagnostic and Interventional Neuroradiology, University of Bern, Inselspital,
Bern University Hospital, Bern, Switzerland
2 Institute for Surgical Technology and Biomechanics, University of Bern,
Bern, Switzerland and 3 Healthcare Imaging A.I. Lab, Insel Data Science
Center, Inselspital, Bern University Hospital, Bern, Switzerland and

# Stratify or Inject: Two Simple Training Strategies to Improve Brain Tumor
Segmentation

###### Abstract

Deep learning methods for brain tumor segmentation are typically trained in an
ad hoc fashion on all available data. Brain tumors are tremendously
heterogeneous in image appearance and labeled training data is limited. We
argue that incorporation of additional prior information, specifically tumor
grade, associated with tumor imaging phenotypes during model training can
significantly improve segmentation performance. Two strategies for
incorporation of tumor grade during model training are proposed and their
impact on segmentation performance is demonstrated on the BRATS 2018 dataset.

††editors: Accepted for MIDL 2019

##  1 Introduction

The segmentation of brain tumors has been a long standing problem in medical
image analysis. Research on this topic has been accelerated greatly through
the availability of public datasets such as the Brain Tumor Segmentation
(BRATS) Challenge dataset [Menze et al.(2015)Menze, Jakab, and et al., Bakas
et al.(2017a)Bakas, Akbari, and et al., Bakas et al.(2017b)Bakas, Akbari, and
et al., Bakas et al.(2017c)Bakas, Akbari, and et al.]. Currently, the best
performing methods for brain tumor segmentation are based on deep learning
[Bakas et al.(2018)Bakas, Reyes, and et al.], with first approaches being
applied in clinically critical areas such as tumor response assessment
[Kickingereder et al.(2019)Kickingereder, Isensee, and et al.] or radiation
therapy planning [Jungo et al.(2018)Jungo, Meier, and et al.]. A general view
in deep learning is that more training data yields better generalization
performance. For tasks in computer vision it was shown that model performance
increases logarithmically based on volume of training data [Sun et
al.(2017)Sun, Shrivastava, and et al.]. Consequently, deep learning
segmentation models are trained on all available data, often neglecting
peculiarities of the data at hand.

The BRATS Challenge has been concerned so far with the segmentation of glioma,
which are primary tumors of the central nervous system. Glioma can be
classified into different tumor grades based on the underlying molecular
characteristics and histology [Louis et al.(2016)Louis, Perry, and et al.]. A
higher grade reflects increasing malignancy and glioma are commonly grouped
into high-grade (grade III/IV) and low-grade glioma (grade I/II). Furthermore,
they exhibit a tremendous genetic and molecular heterogeneity which spans
across tumor grades but also manifests itself within a particular type such as
glioblastoma (grade IV) [Verhaak et al.(2010)Verhaak, Hoadley, and et al.,
Sottoriva et al.(2013)Sottoriva, Spiteri, and et al.]. The underlying
biological configuration of a tumor has been associated with distinctively
different tumor imaging phenotypes [Grossmann et al.(2016)Grossmann, Gutman,
and et al.]. In general, low-grade glioma present much less or no contrast-
enhancement compared to high-grade glioma [Forst et al.(2014)Forst, Nahed, and
et al]. Deep learning methods are confronted with the challenge to
successfully generalize across these different imaging phenotypes.

We hypothesize that brain tumor segmentation performance of deep learning
methods can be improved by utilizing prior information associated with tumor
imaging phenotypes during model training. Thus, we propose two simple training
strategies targeted at tumor grade and evaluate their effectiveness on the
BRATS 2018 dataset using a recently proposed, top-ranked method [McKinley et
al.(2019)McKinley, Meier, and Wiest]. This work is part of a more extensive
study currently being submitted to a journal.

##  2 Methods

Model architecture. The deep learning method corresponds to a shallow U-Net
style model of down and upsampling connections featuring densely connected
blocks of dilated convolutions. For more details on the model architecture we
refer to [McKinley et al.(2019)McKinley, Meier, and Wiest].

Incorporation of tumor grade. We propose two strategies to utilize information
on tumor grade at the stage of model training. The first strategy consists of
stratifying the training data into high-grade glioma (HGG) and low-grade
glioma (LGG) cases and training two separate models. During testing the
respective model is applied to testing data with corresponding tumor grade. As
a second strategy, we propose an injection of the tumor grade. In addition to
feeding the model with imaging data consisting of the co-registered Magnetic
Resonance (MR) sequences, we provide it with a binary input indicating if the
case at hand is either a LGG or HGG. The tumor type is injected as an image
volume, with dimensions identical to the MR sequences and all voxel values set
either to zeros or ones. The model is then trained on this enlarged dataset
(type-aware network). For both strategies, the model architecture and
hyperparameter setting remains unchanged.

In the following, we utilized training data of the BRATS 2018 Challenge, which
includes 75 patients with LGGs and 210 patients with HGGs. We did not include
the testing data in the analysis since tumor grade is blinded for those cases.
A detailed description of the imaging data can be found in [Bakas et
al.(2018)Bakas, Reyes, and et al.]. Four different models were trained using
five-fold cross-validation: i) a baseline model using all available cases
(N=285), ii) a model trained only on HGG data (N=210), iii) a model trained
only on LGG data (N=75), and iv) the type-aware network trained on all cases
(N=285).

##  3 Results & Conclusion

The segmentation performance of the four different deep learning models in
terms of fraction of cases with improved Dice coefficient is shown in Table 1.
If we look at the HGG data alone, the HGG model yields a significantly
improved performance over the baseline model in 58.4% of the cases for
contrast-enhancing tumor and 70.3% of the cases for the tumor core
segmentation. Looking at all the data including also LGG, we see a significant
improvement for the tumor core segmentation in 64.9% of the cases when using
the two models trained on stratified data compared to the baseline model
trained on all available data. Finally, the type-aware network yields a minor
but significant improvement for the segmentation of the tumor core in 53.9 %
of the cases. Figure LABEL:fig:epochs shows that observed improvements are
manifested relatively consistent across all epochs of model training.

Table 1: Ratio in % of better performing subjects compared to baseline. p-values from a one-sided Wilcoxon signed rank test. Bold numbers indicate statistically significant ($p<0.05$) results. (CE = contrast-enhancing tumor) | CE | Core | Tumor
---|---|---|---
LGG vs. Baseline | 41.7 (p=0.877) | 49.3 (p=0.454) | 54.7 (p=0.208)
HGG vs. Baseline | 58.4 (p=0.005) | 70.3 (p=5.659e-09) | 46.7 (p=0.877)
HGG/LGG vs. Baseline | 54.6 (p=0.127) | 64.9 (p=1.441e-05) | 48.8 (p=0.725)
Type-aware vs. Baseline | 53.4 (p=0.321) | 53.9 (p=0.028) | 52.6 (p=0.231)

fig:epochs

Figure 1: Impact of proposed strategies on model performance over all epochs.
From top left to bottom right: LGG vs. Baseline, HGG vs. Baseline, HGG/LGG vs.
Baseline, Type-aware network vs. Baseline.

We have proposed two different strategies on incorporating information of
tumor grade during model training, which are straightforward to apply, and
demonstrated their effectiveness on the BRATS 2018 dataset. While data
stratification yields clear improvements, more advanced network architectures
incorporating prior information about tumor grade beyond injecting it as an
additional input should be investigated further. In addition, our strategies
could also be used in conjunction with networks for tumor typing.

This work was supported by the Swiss National Foundation, grant number 169607.

## References

  * [Bakas et al.(2017a)Bakas, Akbari, and et al.] Spyridon Bakas, Hamed Akbari, and et al.  Advancing the cancer genome atlas glioma MRI collections with expert segmentation labels and radiomic features.  _Scientific Data_ , 4:170117, September 2017a.  10.1038/sdata.2017.117.  URL https://doi.org/10.1038/sdata.2017.117.
  * [Bakas et al.(2017b)Bakas, Akbari, and et al.] Spyridon Bakas, Hamed Akbari, and et al.  Segmentation labels for the pre-operative scans of the tcga-gbm collection, 2017b.  URL https://wiki.cancerimagingarchive.net/x/KoZyAQ.
  * [Bakas et al.(2017c)Bakas, Akbari, and et al.] Spyridon Bakas, Hamed Akbari, and et al.  Segmentation labels for the pre-operative scans of the tcga-lgg collection, 2017c.  URL https://wiki.cancerimagingarchive.net/x/LIZyAQ.
  * [Bakas et al.(2018)Bakas, Reyes, and et al.] Spyridon Bakas, Mauricio Reyes, and et al.  Identifying the best machine learning algorithms for brain tumor segmentation, progression assessment, and overall survival prediction in the BRATS challenge.  _CoRR_ , abs/1811.02629, 2018.  URL http://arxiv.org/abs/1811.02629.
  * [Forst et al.(2014)Forst, Nahed, and et al] D. A. Forst, B. V. Nahed, and et al.  Low-grade gliomas.  _The Oncologist_ , 19(4):403–413, mar 2014.  10.1634/theoncologist.2013-0345.  URL https://doi.org/10.1634/theoncologist.2013-0345.
  * [Grossmann et al.(2016)Grossmann, Gutman, and et al.] Patrick Grossmann, David A. Gutman, and et al.  Imaging-genomics reveals driving pathways of mri derived volumetric tumor phenotype features in glioblastoma.  _BMC Cancer_ , 16(1):611, Aug 2016.  ISSN 1471-2407.  10.1186/s12885-016-2659-5.  URL https://doi.org/10.1186/s12885-016-2659-5.
  * [Jungo et al.(2018)Jungo, Meier, and et al.] Alain Jungo, Raphael Meier, and et al.  Uncertainty-driven sanity check: Application to postoperative brain tumor cavity segmentation.  _CoRR_ , abs/1806.03106, 2018.  URL http://arxiv.org/abs/1806.03106.
  * [Kickingereder et al.(2019)Kickingereder, Isensee, and et al.] Philipp Kickingereder, Fabian Isensee, and et al.  Automated quantitative tumour response assessment of MRI in neuro-oncology with artificial neural networks: a multicentre, retrospective study.  _The Lancet Oncology_ , apr 2019.  10.1016/s1470-2045(19)30098-1.  URL https://doi.org/10.1016/s1470-2045(19)30098-1.
  * [Louis et al.(2016)Louis, Perry, and et al.] David N. Louis, Arie Perry, and et al.  The 2016 world health organization classification of tumors of the central nervous system: a summary.  _Acta Neuropathologica_ , 131(6):803–820, Jun 2016.  10.1007/s00401-016-1545-1.  URL https://doi.org/10.1007/s00401-016-1545-1.
  * [McKinley et al.(2019)McKinley, Meier, and Wiest] Richard McKinley, Raphael Meier, and Roland Wiest.  Ensembles of densely-connected cnns with label-uncertainty for brain tumor segmentation.  In Alessandro Crimi, Spyridon Bakas, and et al., editors, _Brainlesion: Glioma, Multiple Sclerosis, Stroke and Traumatic Brain Injuries_ , pages 456–465, Cham, 2019. Springer International Publishing.  ISBN 978-3-030-11726-9.
  * [Menze et al.(2015)Menze, Jakab, and et al.] Bjoern H. Menze, Andras Jakab, and et al.  The multimodal brain tumor image segmentation benchmark (BRATS).  _IEEE Transactions on Medical Imaging_ , 34(10):1993–2024, oct 2015.  10.1109/tmi.2014.2377694.
  * [Sottoriva et al.(2013)Sottoriva, Spiteri, and et al.] Andrea Sottoriva, Inmaculada Spiteri, and et al.  Intratumor heterogeneity in human glioblastoma reflects cancer evolutionary dynamics.  _Proceedings of the National Academy of Sciences_ , 110(10):4009–4014, 2013.  ISSN 0027-8424.  10.1073/pnas.1219747110.  URL https://www.pnas.org/content/110/10/4009.
  * [Sun et al.(2017)Sun, Shrivastava, and et al.] Chen Sun, Abhinav Shrivastava, and et al.  Revisiting unreasonable effectiveness of data in deep learning era.  In _IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017_ , pages 843–852, 2017.  10.1109/ICCV.2017.97.  URL https://doi.org/10.1109/ICCV.2017.97.
  * [Verhaak et al.(2010)Verhaak, Hoadley, and et al.] Roel G.W. Verhaak, Katherine A. Hoadley, and et al.  Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in pdgfra, idh1, egfr, and nf1.  _Cancer Cell_ , 17(1):98 – 110, 2010.  ISSN 1535-6108.  https://doi.org/10.1016/j.ccr.2009.12.020.  URL http://www.sciencedirect.com/science/article/pii/S1535610809004322.


# A High-Dimensional Particle Filter Algorithm

Jameson<EMAIL_ADDRESS>[ Harvard Universitym1

###### Abstract

Online data assimilation in time series models over a large spatial extent is
an important problem in both geosciences and robotics. Such models are
intrinsically high-dimensional, rendering traditional particle filter
algorithms ineffective. Though methods that begin to address this problem
exist, they either rely on additional assumptions or lead to error that is
spatially inhomogeneous. I present a novel particle-based algorithm for online
approximation of the filtering problem on such models, using the fact that
each locus affects only nearby loci at the next time step. The algorithm is
based on a Metropolis-Hastings-like MCMC for creating hybrid particles at each
step. I show simulation results that suggest the error of this algorithm is
uniform in both space and time, with a lower bias, though higher variance, as
compared to a previously-proposed algorithm.

,

t1Thanks to Mira Bernstein, Pierre Jacob, and Luke Miratrix for constructive
criticism of this manuscript.

##  1 Background

Filtering problems arise in many applied contexts, whenever noisy observations
over time must be combined, using an explicit dynamical model, into a best-
guess distribution of a current state. In cases where the system being modeled
involves processes over a large spatial extent, such as models of weather or
other large-scale fluid dynamics, filtering is also called data
assimilation.[3] This is an active area of research, with broad applications
in predictive geoscience[4][20] and robotics[19]. In fact, it is considered to
be among the central problems in both of these disciplines.111 According to
the papers cited above, there is “much focus in the [geoscience] literature on
the assimilation of data and numerical models pertain[ing] to the sampling of
high-dimensional probability density functions” [4], and “The SLAM
[simultaneous location and mapping] problem is generally regarded as one of
the most important problems in the pursuit of building truly autonomous mobile
robots.”[19]

The basic filtering problem is as follows. We model the state of the system at
time $t$ as a random variable $X_{t}$. In our context, $X_{t}$ will have
values
$x_{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr
l\\\ t\crcr}}}$ at each spacial locus $l$, for a large number of loci. We
assume $X_{0},...,X_{T}$ form a Markov chain with known and sampleable
densities for both the initial state ($\pi_{0}$) and transition function
($\mathsf{P}$, which maps states or densities at time $t-1$ to densities at
time $t$). We also have a series of observations $Y_{1},...,Y_{T}$, and we
assume that each $Y_{t}$ depends only on the corresponding $X_{t}$ according
to a known and sampleable observation density $f(Y_{t}|X_{t})$. This is shown
graphically in Figure 1.

$X_{0}$$X_{1}$$\dots$$X_{T-1}$$X_{T}$$Y_{0}$$Y_{1}$$\dots$$Y_{T-1}$$Y_{T}$
Figure 1: Graphical model of a low-dimensional filtering problem.

At each time step $t$, we wish to estimate the filtering distribution: that
is, the probability density $\pi_{t}$ of $X_{t}|\\{Y_{1},...,Y_{t}\\}$.
Because the $X_{t}$ are Markovian, and $Y_{t}$ depends only $X_{t}$, we can
write $\pi_{t}$ recursively as

$\pi_{t}(\cdot)=\bm{E}_{X_{t-1}\sim\pi_{t-1}}[P(X_{t}\in\cdot|Y_{t},X_{t-1})]$
(1.1)

To minimize subscripts, we adopt the following notation:

  * •

We abbreviate the filtering distributions $\pi_{t-1}$ and $\pi_{t}$ by $\tau$
and $\pi$ respectively.

  * •

We abbreviate samples from $X_{t-1}$ and $X_{t}$ by $\bm{x}$ and $\bm{z}$
respectively.

  * •

When $\bm{y}_{t-1}$ is not relevant, we write $\bm{y}$ for $\bm{y}_{t}$.

  * •

Superscripts should not be read as exponentiation for these and similar
entities.

The recursive formula for $\pi_{t}$ suggests the possibility of online
calculation, with only a constant computing time required to update from
$\tau=\pi_{t-1}$ to $\pi=\pi_{t}$. However, unless we assume a particular
parametric form for $\pi$, there is no finite set of sufficient statistics
that could stand in for the full distribution. Thus, aside from very simple
special cases, exact calculation is impossible; we look for an approximation
instead.

A widely-used recursive algorithm for approximating the filtering distribution
is the bootstrap particle filter. Assuming we have a sampleable distribution
$\hat{\tau}$ at time $t-1$ that approximates the true filtering distribution
$\tau$, we proceed as follows:

  1. 1.

Sample $M$ iid particles $\bm{x}^{1..M}$ from $\hat{\tau}$. (Note that if we
have been following the algorithm up to step $t-1$, then $\hat{\tau}$ takes
the form of step (3) below.)

  2. 2.

For each $\bm{x}^{i}$, progress it to get candidate particle
$\bm{z}^{i}\sim\mathsf{P}\bm{x}^{i}$.

  3. 3.

Find weights $w^{i}\equiv f(\bm{y}|\bm{z}^{i})$. The set of weighted particles
forms

$\hat{\pi}\equiv\frac{\sum_{i=1}^{M}w^{i}\delta(\bm{z}^{i})}{\sum_{i=1}^{M}w^{i}}.$

Here $\delta(a)$ is the Dirac delta density; for example,
$\frac{1}{2}(\delta(0)+\delta(1))$ is the Bernoulli distribution with
$p=\frac{1}{2}$.

A key property of the particle filter algorithm is that, for large enough $M$,
the Monte Carlo error remains under control, even as the time steps
accumulate.[6] Specifically, let $\hat{\pi}_{t}^{M}$ be the approximation to
$\pi_{t}$ obtained using $M$ particles, as above. Let $\mathcal{F}$ be the set
of functions from the domain of $\pi_{t}$ to $(-1,1)$. For each
$\mathrm{f}\in\mathcal{F}$, denote $\bm{E}_{X\sim{\pi_{t}}}(\mathrm{f}(X))$
and $\bm{E}_{X\sim{\hat{\pi}_{t}^{M}}}(\mathrm{f}(X))$ by
$\pi_{t}(\mathrm{f})$ and $\hat{\pi}_{t}^{M}(\mathrm{f})$ respectively. Then

$\sup_{\mathcal{F}}\bm{E}|\pi_{t}(\mathrm{f})-\hat{\pi}^{M}_{t}(\mathrm{f})|\leq\frac{C}{\sqrt{M}},$
(1.2)

for some constant $C$ that does not depend on $t$. The outer expectation here
is taken over the randomness of the algorithm itself; that is, considering the
distribution $\hat{\pi}$ as itself a random variable, while $\pi$ is
fixed.[14, p. 2814]

Now suppose that we are interested in modeling processes with a large spatial
extent. For instance, in a weather model, one might use a lattice of points
that cover the region of interest, with various continuous values
(temperature, humidity, pressure, wind, etc.) recorded at each locus. If
$X_{t}$ contains information about $L$ separate spatial loci, and the state
space at each locus has dimension $K$, then the full state space of $X_{t}$
has dimension $KL$. In practical applications, this can easily be $10^{7}$ or
more.[21]

To see why, consider a schematic diagram of the model (Figure 2), where the
state of the system at time $t$ and locus $l$ is denoted
$x_{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\cr
l\\\ t\crcr}}}$.

$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{0}\end{subarray}}$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{1}\end{subarray}}$
$\dots$$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{T}\end{subarray}}$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{0}\end{subarray}}$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{1}\end{subarray}}$ $\dots$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{T}\end{subarray}}$ $\dots$ $\dots$ $\dots$ $\dots$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{L}\\\
\scriptstyle{0}\end{subarray}}$$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{L}\\\
\scriptstyle{1}\end{subarray}}$ $\dots$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{L}\\\
\scriptstyle{T}\end{subarray}}$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{0}\end{subarray}}$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{1}\end{subarray}}$ $\dots$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{0}\\\
\scriptstyle{T}\end{subarray}}$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{0}\end{subarray}}$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{1}\end{subarray}}$ $\dots$
$\scriptstyle{x}_{\begin{subarray}{c}\scriptstyle{1}\\\
\scriptstyle{T}\end{subarray}}$ $\dots$ $\dots$ $\dots$ $\dots$
$\scriptstyle{y}_{\begin{subarray}{c}\scriptstyle{L}\\\
\scriptstyle{t}\end{subarray}}$ Figure 2: Graphical model of a (simple) high-
dimensional filtering problem.

Note the following assumptions implicit in the diagram:

  * •

$(y_{\begin{subarray}{c}l\\\ t\end{subarray}}|x_{\begin{subarray}{c}l\\\
t\end{subarray}})\perp\\!\\!\\!\perp x_{\begin{subarray}{c}k\neq l\\\
t\end{subarray}}$. This assumption will be used in the following for
simplicity, although I believe it can be relaxed in practice with only minor
additional complications.

  * •

More crucially, $x_{\begin{subarray}{c}l\\\ t\end{subarray}}$ depends on
$x_{\begin{subarray}{c}k\\\ t-1\end{subarray}}$ only for $k$ in some small
spatial “neighborhood” $\mathcal{N}(l)$ of $l$. The precise composition of
$\mathcal{N}(l)$ depends on model assumptions as well as the way that the $L$
loci are positioned in space. The diagram depicts the case where the loci are
all laid out along a single spatial dimension, so that the immediate neighbors
of locus $l$ are loci $l-1$ and $l+1$. In practical applications, the loci
would more likely be connected in a 2D or 3D grid.

I will discuss this locality of dynamics assumption further below, as it is
key to the performance of the algorithm I propose in this paper.

Technically speaking, the error bounds in 1.2 still apply: errors are stable
over time and inversely proportional to the square root of the number of
particles ($\propto 1/\sqrt{M}$). But it is widely recognized that the
bootstrap particle filter is no longer a practical solution in this context,
due to weight degeneracy.[18][5][1] The problem is that, in the resampling
step, the majority of the weight will tend to be carried by only a small
fraction of the proposed particles. To see why, note that the log likelihood
of each particle is the sum of its log likelihood at each locus. Since these
terms are roughly independent, as $L$ increases, the empirical distribution of
log likelihoods of the particles comes to resemble a Gaussian distribution
with variance that scales linearly with $L$. Thus the weights come to be
distributed approximately according to a log-normal distribution, whose
skewness increases exponentially with $L$. Thus, the fraction of particles
with above-average weight will shrink exponentially with $L$. 222Bickel et
al.[5] formalize this line of argument, showing that the variance of the log
likelihood may be seen as an estimate of the effective state dimension
depicted by the measurements $\bm{y}$.

In fact, without an exponentially large number of particles, not only will one
of them tend to have more weight than all the others, but it is likely that
there is some missing value whose weight would dominate even our best
particle. At this point the particle filter ceases to be a useful
approximation of the filtering distribution. That is, although the constant
$C$ in Equation 1.2 does not depend on time, it does grow exponentially with
$L$. This is what is known as curse of dimensionality for particle filters.

###  1.1 Existing state of the art (Rebeschini and van Handel, 2015)

There have been various proposals for dealing with this curse of
dimensionality in general. In fact, there are three different recent survey
articles reviewing and comparing these, by Septier and Peters[16], Morzfield
et al.[11], and Farchi and Bocquet.[7] Some of these prior methods do not use
the locality of dynamics assumption, which I believe limits their
effectiveness. This includes Gilks and Berzuini,[8] who suggest rejuvenating
particles with MCMC steps, targeted to the filtering distribution, to avoid
the duplication problem from resampling; and Goodsill and Clapp,[9] who
propose using bridging densities such as annealed quasi-filtering densities to
solve the problem of lack of overlap between the progressed density
$\mathsf{P}\tau\equiv(\bm{z}|\bm{y}_{1},...,\bm{y}_{t-1})$ and the likelihood
$f(\bm{y}_{t}|\bm{z})$.

The two previous proposals that do use locality of dynamics come from
Poterjoy[12][13] and from Rebeschini and van Handel.[14] Of these two,
Rebeschini and van Handel’s proposal is more generally applicable, so I will
explain it further below. Poterjoy, on the other hand, suggests a scheme that,
as given, is limited to situations of sparse observations; it uses estimated
covariance matrices to blend resampled particle filter values at a local scale
with unresampled values at a meso-scale.

Rebeschini and van Handel’s proposal is called the block particle filter; also
sometimes termed the localized particle filter. In simple terms, they replace
the global resampling step of the bootstrap particle filter with a local
resampling step which constructs new particles by resampling neighborhoods
independently.

They begin their theoretical discussion by offering an overall point of view
of the problem which very much inspired the current work. They focus on the
decay of correlations between local values as spatial distance increases,
which they say is “in essence a spatial counterpart of the much better-
understood [temporal] stability property of nonlinear filters”. This decay of
correlations, discussed further below, is a product of the locality-of-
dynamics assumption shown in the diagrams above.

Rebeschini and van Handel partition the loci ${1,...,L}$ of the progressed
particles into $J$ zones $\\{\mathcal{Z}_{j}\\}$, where each zone consists of
a small number of (contiguous) loci. They then weight and resample values from
each zone independently. This produces what I would call “Frankenstein”
particles, sewn together from pieces which tend to fit well with observations
locally, but which may often come from progressing different particles from
time $t-1$.

The precise steps of the block particle filter algorithm are as follows:

  1. 1.

Given $\hat{\tau}$, a sampleable distribution that approximates the ideal
filtering distribution $\tau$, sample $M$ iid particles $\bm{x}^{1..M}$ from
$\hat{\tau}$. (Note that if $\hat{\tau}$ takes the form of the output of step
(3) below, then this amounts to sampling, independently with replacement, a
$k^{i}_{j}$ for each zone $Z_{j};1\leq j\leq J$ and particle $i$ with
probability $w^{k^{i}_{j}}_{Z_{j}}/\Sigma_{n}w^{n}_{Z_{j}}$; then putting
those together to make the final particles, so that $l\in Z_{j}\Rightarrow
z^{i}_{l}=z^{k^{i}_{j}}_{l}$.)

  2. 2.

For each $\bm{x}^{i}$, progress it to get $\bm{z}^{i}\sim\mathsf{P}\bm{x}$.

  3. 3.

Find weights for each particle for each zone: $w^{i}_{Z_{j}}={\prod_{l\in
Z_{j}}}f(y_{l}|z^{i}_{l})$. Then

$\hat{\pi}\equiv\operatorname*{{{{{\mathbf{\otimes}}}}}}_{j=1}^{J}\frac{\sum_{i=1}^{M}w^{i}_{Z_{j}}\delta(\bm{x}^{i}_{Z_{j}})}{\sum_{i=1}^{M}w^{i}_{Z_{j}}}$

Intuitively, cutting the progressed particles into zones and reconstituting
them is a way to solve the exponential curse of dimensions. However, it gives
rise to a different problem: it breaks any covariances across zone boundaries,
whether those come from covariances at the $t-1$ time step or are induced by
the transition kernel using a history that overlaps the boundary. These broken
covariances lead to error at the boundaries, which does not disappear even as
number of particles goes to infinity.

It could also then lead to unrealistic dynamics near the boundaries at later
time steps, especially if the forward density operator $\mathsf{P}$ is
nonlinear. [14, p. 2829] For instance, imagine a weather model in which the
hypothetical air pressure in a particle varied reasonably within each zone,
but a discontinuity at the zone boundary led to a prediction of a tornado
forming in the next time step. Note that the algorithm proposed in this paper
avoids such discontinuities, but for unrelated reasons can be ill-suited to
modeling models with nonlinear dynamics; I will address this issue in later
work with a proposed extension to my algorithm.

For the block particle filter method to be useful, covariances between local
values must tend to decay with distance. If such decay of correlations holds,
then, far enough from a zone boundary, the dynamics return to normal.

This intuition helps explain the error bounds that these researchers prove
their algorithm obeys. They show that the error for the value at a given locus
$l$ satisfies

$|||\pi-\hat{\pi}^{M}|||_{\\{l\\}}\leq\alpha(\frac{e^{\beta|Z_{l}|}}{\sqrt{M}}+e^{-\gamma\inf|l-b\in
Z_{l}^{C}|})$ (1.3)

where the constants $\alpha,\beta,\gamma>0$ do not depend on $t$. (They define
a norm for this purpose which measures the distance between random
distributions; I will not reproduce this definition here, as in this this
result is merely a guide for intuition.)

One can see that there is a tradeoff: using smaller zones and/or more
particles allows better guesses for a given zone to control the term
$\frac{e^{\beta|Z_{l}|}}{\sqrt{M}}$, while using larger zones makes it
possible to move away from boundary effects and control the term
$e^{-\gamma\inf|m-b\in Z_{l}^{C}|}$. In practice, Rebeschini and van Handel
give a simple example where it would still be possible to control average
error per locus, based on finding an optimal balance between adding particles
and increasing neighborhood size, but in that example the inverse of their
error bound grows only logarithmically with the computing resources/number of
particles — in other words, as tolerances tighten, the required number of
particles can still grow exponentially. Thus, although their algorithm in
practice gives error far lower than that of the bootstrap particle filter, it
has not fully overcome the problem of needing exponential computing cost as
dimension grows, especially if one wishes to achieve fixed error bounds that
are lower than what comes easily from a moderate neighborhood size.

##  2 The Finkelstein solution

In this section, I will sketch out the basic outlines of a recursive algorithm
in which each particle at time $t$ is composed of values at different loci
which are drawn from state vectors progressed from different particles at time
$t-1$. The choice of which values for a given locus combine with which values
at other loci is made by running a separate Metropolis-Hastings MCMC to create
each composite particle, proposing to replace one locus value at each MCMC
step. What acceptance ratio $\rho$ to use for those proposals will be
discussed in later sections.

I name this the Finkelstein algorithm, after the character Sally Finkelstein
from the movie “The Nightmare Before Christmas”. I have already compared the
block particle filter to a Frankenstein solution, in which the progressed
particles are chopped up and then randomly sewn back together. In that
algorithm, the suitability of the values in each zone of each particle is
measured against observation, but not against the other zones on which it
borders. Finkelstein can improve on this. Though herself originally a
Frankenstein’s-monster-like creation of a stereotypical mad doctor, now that
she has been animated, Sally is able to lose her body parts and sew them back
on, and thus presumably to choose for herself only those body parts that best
fit together. In my terms, Finkelstein would be able to run an MCMC process on
her own body, targeting whatever distribution she pleases.

A key aspect of this algorithm is the Metropolis-Hastings acceptance ratio
$\rho$. When choosing a formula for such a ratio, the key question is, what
distribution do we wish to target? I’ll begin by showing an algorithm that
targets the natural unnormalized density:

$\displaystyle\begin{split}f_{\bm{x}\sim\tau}(\bm{z}|\bm{y})&=\int
f(\bm{z}|\bm{y},\bm{x})\tau(\bm{x})d\bm{x}\\\ &\propto\int
f(\bm{y},\bm{z})f(\bm{z},\bm{x})\tau(\bm{x})d\bm{x}\\\
&=[\Pi_{l}f(y_{l}|z_{l})]\int[\Pi_{l}f(z_{l}|\bm{x})]\tau(\bm{x})d\bm{x}\end{split}$
(2.1)

It will turn out that targeting this density still suffers a similar curse of
dimensionality as the bootstrap particle filter, so I will modify the
algorithm, such that its stationary distribution is not precisely the above
expression. Still, this expression is still the starting point; by
approximately targeting it, I approximately target $\pi$.

Here are the steps of the basic Finkelstein algorithm. I do not include a
formula for the acceptance probability $\rho$ here; I will develop and discuss
several alternatives for such a formula, based on modifications of the density
above, in the following sections.

  1. 1.

Assume we have $\hat{\tau}^{M}$, a sampleable distribution which in some sense
approximates the ideal filtering distribution $\tau$, and which must be of the
form $\frac{1}{M}\sum_{i=1}^{M}\delta(\bm{x}^{i})$. Note that unlike the cases
of the bootstrap and block particle filters, the final $\hat{\pi}$ produced by
this algorithm already has equally-weighted particles; resampling is not
required.

  2. 2.

For each particle $\bm{x}^{i}$, progress it to get a full particle
$\bm{\tilde{z}}^{i}\sim\mathsf{P}\bm{x}^{i}$ whose local values are known as
${\tilde{z}}^{i}_{l}$. (In later steps, I will assume for simplicity that
there are no duplicate values at any locus, so $i\neq
j\Rightarrow{\tilde{z}}^{i}_{l}\neq{\tilde{z}}^{j}_{l}$, but relaxing this
assumption should be straightforward if necessary.)

  3. 3.

Find likelihood weights for each such local value, denoted $w^{i}_{l}\equiv
f(y_{l}|\tilde{z}^{i}_{l})$; and forward densities conditional on $\bm{x}^{j}$
for all $j$ (including $j=i$), denoted $f^{j\rightarrow i}_{l}\equiv
f_{\mathsf{P}}(\tilde{z}^{i}_{l}|\bm{x}^{j})$.

  4. 4.

In parallel, for $k=1,\dots,M$, do the following:

    1. (a)

Create a new proposal particle by independently sampling each locus of a
vector $\bm{\iota}^{0}\in\\{1..M\\}^{L}$. This vector specifies the source for
the value of $\tilde{z}$ that is being considered at each locus; that is,
after running the MCMC for $S$ steps, the final value $\bm{\iota}^{S}$ will be
used to define $\bm{z}^{k}$ by setting
$z_{l}^{k}=\tilde{z}_{l}^{\iota^{S}_{l}}$. The initial sampling at each locus
uses the probabilities

$P(\iota^{0}_{l}=i)=w^{i}_{l}/\sum_{j}w^{j}_{l}.$

(Note that these initial sampling probabilities are arbitrary and, if the MCMC
successfully runs until convergence, irrelevant. The probabilities above
represent a reasonable starting point that should converge reasonably well,
but it may be possible to get even faster convergence through some sampling
scheme that is not independent across loci.)

    2. (b)

Run a Metropolis-Hastings MCMC chain targeting an approximation of the
filtering distribution, for $s=1,...,S$ steps to (assumed) convergence:

      1. i.

Choose a spatial locus $\lambda(s)\in\\{1,\dots,L\\}$ uniformly at random. For
brevity, I will refer to this as $\lambda$, suppressing the dependency on $s$,
in the steps that follow.

      2. ii.

Sample a proposed particle $\iota^{*}$ from which to draw the replacement
value $\tilde{z}^{\iota^{*}}_{\lambda}$ for locus $\lambda$, with probability

$P(\iota^{*}=i)=w^{i}_{\lambda}/\sum_{j}w^{j}_{\lambda}.$

As with $\lambda$ itself, I am omitting here the subscript $s$, even though
this will be resampled at each step. Note that unlike $\bm{\iota}^{s}$, which
is a vector of one integer per locus, this $\iota^{*}$ is only one integer,
which determines the source for the proposed value only at locus $\lambda$. So
for convenience in the formulas below, I will also define the vector
$\bm{\iota}^{**}$ such that $\iota^{**}_{\lambda}=\iota^{*}$ and $\forall
k\in\\{1,...,\lambda-1,\lambda+1,...,L\\}:\iota^{**}_{k}=\iota^{s-1}_{k}$

      3. iii.

Accept this proposed change with M-H probability:
$1\wedge\rho(\bm{\iota}^{s-1},\lambda,\iota^{*})$, where $\rho$ is defined
below. In case of acceptance, $\bm{\iota}^{s}:=\bm{\iota}^{**}$; Otherwise, in
case of rejection, make no change, so that $\bm{\iota}^{s}:=\bm{\iota}^{s-1}$.

    3. (c)

Let $z_{l}^{k}=\tilde{z}_{l}^{\iota^{S}_{l}}$.

We would like to tune the proposal and acceptance probabilities so that this
algorithm targets the distribution $(\bm{z}|\bm{y},\\{\bm{x}^{i}\\})$. Insofar
as we succeed, the full algorithm will be very similar to a bootstrap particle
filter, which similarly targets that distribution. Thus, on a purely intuitive
level, it is unsurprising that this should work if the MCMC does. But the
whole process depends on the validity of the MCMC, which in turn depends on
the M-H acceptance probability $\rho(\bm{\iota}^{s-1},\lambda,\iota^{*})$,
left unspecified above.

##  3 Developing a working formula for $\rho$

In this section, I’ll first develop some motivating understanding of what
$\rho$ should be like, then develop two specific formulas for $\rho$:

  1. 1.

The first formula, $\rho_{\mathrm{full}}$, is for illustrative purposes only.
Though it is, by construction, asymptotically correct–that is, the algorithm
using $\rho_{\mathrm{full}}$ approaches the correct filtering distribution as
the number of particles $M$ approaches infinity–it is unsuitable for use in
practice. Not only does it lead to impractically high computational costs for
a specific $M$, it also suffers from similar dimensionality problems as the
bootstrap particle filter.

  2. 2.

The second formula $\rho_{\mathrm{local}}$ only considers values of the MCMC
in some local neighborhood of $\lambda(s)$. It thus does not suffer the weight
degeneracy problem of the bootstrap particle filter or $\rho_{\mathrm{full}}$,
while still being possible to calculate in computing time that’s polynomial in
number of particles, and linear in dimension and number of time steps.

In the next section, I’ll develop a further formula for $\rho$ that improves
the computational characteristics, as well as proposing some other
computational optimizations.

###  3.1 $\rho_{\mathrm{full}}$: full acceptance probability

I will begin by using a standard Metropolis-Hastings ratio to derive a
$\rho_{\mathrm{full}}$ that targets the (unnormalized) density 2.1. That
expression is promising in one sense: it suggests that one can judge the fit
of local proposals drawn from two different particles using an expression
involving the transition kernel forward density. However, the fact that it
involves a sum over all particles, of a product over all loci, makes
$\rho_{\mathrm{full}}$ computationally unworkable in practice. And even if
there were enough available computational power to calculate this at every
step of an MCMC, taking a product over loci would lead to similar problems as
the naive high-dimensional bootstrap particle filter has.

If we propose replacement $z^{i}_{l}$ for the value at a given locus $l$ with
probability proportional to some value ${v}^{i}_{l}$, we can construct a
Metropolis-Hastings acceptance probability in the usual way, by multiplying a
ratio of target densities (proposed over current) by a ratio of proposal
densities (current over proposed).

(To avoid nested subscripts/superscripts in the following, I leave out
redundant indices for locus; thus using the notation
$f^{j\rightarrow\iota^{**}}_{k}$ rather than
$f^{j\rightarrow\iota^{**}_{k}}_{k}$, using $f^{j\rightarrow\iota(s-1)}_{k}$
rather than $f^{j\rightarrow\iota^{s-1}_{k}}_{k}$, and using
$w^{\iota(s-1)}_{k}$ rather than $w^{\iota^{s-1}_{k}}_{k}$.)

$\rho_{\mathrm{full}}\equiv\frac{[{w^{\iota^{*}}_{\lambda}}{\Pi_{k\neq\lambda}w^{\iota(s-1)}_{k}}]~{}\Sigma_{j=1}^{M}[\Pi_{k}f^{j\rightarrow\iota^{**}}_{k}]}{[{w^{\iota(s-1)}_{\lambda}}{\Pi_{k\neq\lambda}w^{\iota(s-1)}_{k}}]~{}\Sigma_{j=1}^{M}[\Pi_{k}f^{j\rightarrow\iota(s-1)}_{k}]}\cdot\frac{{v}^{\iota(s-1)}_{\lambda}~{}\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s-1)}_{\lambda}}{{v}^{\iota^{**}}_{\lambda}~{}\Sigma_{j=1}^{M}f^{j\rightarrow\iota{**}}_{\lambda}}$
(3.1)

Let’s consider these terms from left to right, taking the version of each term
as it appears in the numerator:

  1. 1.

${w^{\iota^{*}}_{\lambda}}$: The likelihood at the locus in question; the term
which depends on $\bm{y}$. I will set ${v}^{\iota^{*}}_{\lambda}$ so as to
cancel this out.

  2. 2.

${\Pi_{k\neq\lambda}w^{\iota(s-1)}_{k}}$: The likelihoods at other loci. These
cancel out naturally.

  3. 3.

$\Sigma_{j}[\Pi_{k}f^{j\rightarrow\iota^{**}}_{k}]$: This sum of products term
is the heart of the calculation at each MCMC step. Each product is the forward
likelihood of the current/proposed hybrid particle conditional on a given
history; the sum of products is proportional to the forward likelihood of the
current/proposed hybrid particle conditional on $\hat{\tau}^{M}$.

  4. 4.

${v}^{\iota(s-1)}_{\lambda}$: Weights that define the proposal distribution
and can be chosen arbitrarily (up to normalization).

  5. 5.

$\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s-1)}_{\lambda}$: The part of the
proposal density that’s due to $\hat{\tau}^{M}$; the probability density that
a given value would have been in $\\{\bm{\tilde{z}}^{i}_{\lambda}:0<l<L\\}$ to
be available to be sampled. To someone used to other variations of particle
filtering algorithms, it may seem counterintuitive to include this term;
usually, taking advantage of the forward density is a key aspect of how the
algorithm works, not something that needs canceling out. However, from the
perspective of a Metropolis-Hastings construction, it is necessary to include
this in order to target the intended (unnormalized) density 2.1. On a more
intuitive level, one might note that the forward density values
$f^{j\rightarrow\iota(s-1)}_{\lambda}$ for each $j$ appear both here and in
the sum of products term in the denominator; so if one did not include this
term, that would in a sense be double counting these forward density values by
allowing them to cause an increased proposal density and then also increase
the acceptance ratio.

The unnormalized proposal weights ${v}^{i}_{l}$ can be set at will; I will let
${v}^{i}_{l}\equiv w^{i}_{l}$, so that these terms cancel out. Now $\rho$ is
just a ratio of sums of products, multiplied by a ratio of the forward mean
proposal densities $\bar{f}^{i}_{l}$. I could have set ${v}^{i}_{l}$ to also
cancel out $\bar{f}^{i}_{l}$, but as seen later, this would slow convergence.

The quantities $\bar{f}^{i}_{l}$ do not change for different MCMC chains or
for different steps in each chain and can be precalculated in $O(LM)$ time. We
thus have

$\displaystyle\rho_{\mathrm{full}}$
$\displaystyle\equiv\frac{\Sigma_{j=1}^{M}[\Pi_{k=1}^{L}f^{j\rightarrow\iota^{**}}_{k}]}{\Sigma_{j=1}^{M}[\Pi_{k=1}^{L}f^{j\rightarrow\iota(s-1)}_{k}]}\frac{\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s-1)}_{l}}{\Sigma_{j=1}^{M}f^{j\rightarrow\iota^{**}}_{l}}$
(3.2)

By construction, this acceptance probability satisfies the conditions for
detailed balance of a Metropolis-Hastings MCMC, and thus converges to the
desired target stationary distribution 2.1 under standard regularity
assumptions. At convergence, the algorithm will give $M$ samples from the
correct target distribution, but with the constraint that the value for each
sample at each locus must be available in $\tilde{\bm{z}}^{1..M}$. As
$M\rightarrow\infty$, the set of values available at each locus will become
dense, so the conditionality will not be restrictive. Thus, asymptotically,
one would expect that each $\bm{z}^{i}$ should be a sample from the correct
filtering distribution, conditional on $\bm{x}\sim\hat{\tau}^{M_{t-1}}$.

Yet $\rho_{\mathrm{full}}$ is not useful in practice, for two reasons. First,
on a relatively trivial level, actually calculating the sum of products terms
once for each step of the MCMC would be computationally prohibitive; though
not exponential, the resources required would be extreme.

More importantly, unless $M_{t-1}$ is exponentially high,
$\mathsf{C}_{t}\mathsf{P}\hat{\tau}^{M_{t-1}}$ is not a good approximation of
$\mathsf{C}_{t}\mathsf{P}\tau$, because of a curse of dimensionality very
similar to that which causes the bootstrap particle filter to fail in high
dimensions. In $\rho_{\mathrm{full}}$, the acceptance probability is based on
a sum over histories of products over loci of likelihoods. Following a similar
logic as Bickel et al.[5], discussed above, for showing weight degeneracy in
the bootstrap particle filter, we see that, assuming that the likelihoods
associated with distantly-separated loci are roughly independent, then as
dimension increases the distribution of these products over loci will approach
a log-normal. Since the variance of that distribution will grow with
dimension, the sum is likely to be degenerate unless number of particles grows
exponentially with dimension; just one history particle will contribute more
to the sum than all others put together.

Consider the example of a weather model of the continental United States,
where imperfect measurements of atmospheric conditions are taken daily over a
set of cities. In this case, the sum would be degenerate because, although any
given proposal particle (weather map for today) might accord well with a given
history (possible weather map for yesterday) for some cities, you’d
nevertheless need to consider an exponentially large number of possible
histories before finding one which accords well across all cities with a given
realistic present.

In essence, rather than resolving the high-dimensionality problem at time $t$,
we’ve merely pushed it off to time $t-1$; because of the high variance of the
product terms $\Pi_{k=1}^{L}f^{j\rightarrow\iota^{**}}_{k}$, the sums will
tend to be dominated by the product term for a single history $j$, losing most
of the benefits of a high number of particles.

###  3.2 $\rho_{\mathrm{local}}$: Simplifying the product terms by focusing on
local neighborhoods

The main computational burden of calculating $\rho_{\mathrm{full}}$ comes from
the sum over histories of a product over loci. To deal with these
computational issues, as well as with the degeneracy of the sum, it would be
good to take this product over fewer terms. To do so, I restrict the product
over loci to only consider loci in some neighborhood of the locus $l$ which
the proposal would change.

This idea gains some support from the decay of correlations property of that
Rebeschini and van Handel (2015) demonstrate. This is a complex issue which
occupies a significant portion of their paper, but to summarize briefly: they
assume particle filters with local dynamics, and both forward densities and
observation likelihoods that are strongly bounded away from zero and infinity,
Given those assumptions (which they argue are probably stronger than necessary
in most practical cases), they show that changing the value at locus $k$
cannot change the conditional distribution of the value at $l$ by more than a
quantity that falls exponentially as the distance between $k$ and $l$
increases. Thus, it would seem logical that, in calculating an acceptance
probability to target the distribution at $l$, one may safely ignore faraway
loci $k$.

To use this idea for the Finkelstein algorithm, assume there is a natural
distance metric $d(l,k)$ over loci; for instance, if loci were arranged in a
square lattice, $d(l,k)$ could be the $\ell_{1}$-distance. Use this distance
to define neighborhood balls
$\mathcal{B}_{r}(\lambda)\equiv\\{l:d(\lambda,l)\leq r\\}$, and use the
natural notation that
$\bm{x}_{\mathcal{B}_{r}(\lambda)}\equiv\\{x_{l}:l\in\mathcal{B}_{r}(\lambda)\\}$.
Thus, the new $\rho$ would be:

$\rho_{\mathrm{local}}\equiv\frac{\Sigma_{j=1}^{M}[\Pi_{k\in\mathcal{B}_{r}(\lambda)}f^{j\rightarrow\iota^{**}}_{k}]}{\Sigma_{j=1^{M}}[\Pi_{k\in\mathcal{B}_{r}(\lambda)}f^{j\rightarrow\iota(s-1)}_{k}]}\frac{\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s-1)}_{\lambda}}{\Sigma_{j=1}^{M}f^{j\rightarrow\iota^{*}}_{\lambda}}$
(3.3)

If this works, it will have finally conquered the curse of dimensions. For any
locus $l$, there are only $|\mathcal{B}_{r}(l)|$ terms in each product; a
quantity which does not depend on the overall dimension of the problem, only
on the local connectivity. We can therefore choose a fixed $M$ large enough
such that this sum of products is not degenerate (not dominated by just one of
the products) for $N\equiv\max_{l}|\mathcal{B}_{r}(l)|$ loci.

However, it should be noted that with this acceptance probability, the
algorithm is no longer strictly speaking Metropolis-Hastings. In particular,
the overall MCMC is no longer guaranteed to obey detailed balance. If one
repeatedly replaced the values of a single locus $l$, the MCMC would, by the
standard Metropolis-Hastings construction, show detailed balance at a unique
stationary distribution with a density proportional to:

$f_{l}(z_{l}|z_{1},...,z_{l-1},z_{l+1},...,z_{L})\equiv(\Pi_{\lambda=1}^{L}\mathds{1}_{z_{\lambda}\in\\{\tilde{z}_{\lambda}\\}})\Sigma_{j=1}^{M}[\Pi_{k\in\mathcal{B}_{r}(l)}f_{\mathsf{P}}(z_{k}|\bm{x}^{j})]$
(3.4)

This function is not only of $z_{l}$, but of all $z_{k}$ such that
$k\in\mathcal{B}_{r}(l)$. However, since this density is not the same for two
different values of $l$, this detailed balance can and will break down. The
MCMC is still uniformly ergodic, so a unique stationary distribution still
exists; but without detailed balance, we lack the nice guarantees that
Metropolis-Hastings would offer as to what that target distribution is. At
present, then, my use of $\rho_{\mathrm{local}}$, and all later versions of
$\rho$ that build on it, is based on empirical validation, as seen below in
the simulation section, not rigorous theory.

##  4 Computational optimizations

###  4.1 $\rho_{\mathrm{sampled}}$: a version of $\rho_{\mathrm{local}}$ which
replaces numerator and denominator by unbiased estimators

Running the Finkelstein algorithm with acceptance probability
$\rho_{\mathrm{local}}$ does not require exponential computation, but even
polynomial amounts of computation can be daunting in practice. Recall that the
of sums over all histories in $\rho_{\mathrm{local}}$ are proportional to the
forward likelihood conditional on $\hat{\tau}^{M}$, the $M$-particle
approximation of the filtering distribution $\tau$. At each step, we can save
computation by, effectively, estimating $\tau$ with only an arbitrary fixed
number $H<<M$ particles; that is, by using only $H$ history terms for these
sums and using those to get unbiased Horvitz-Thompson estimators of the
totals. This leads to $\rho_{\mathrm{sampled}}$. Note that the specific $H$
history particles used will change from step to step and locus to locus in the
MCMC, thus taking advantage of the full $M$ particles in equilibrium.

The idea of using unbiased estimators to calculate the Metropolis-Hastings
acceptance ratio is not new, and, as Andrieu and Roberts 2009[2] show, this
can be made to conserve the stationary distribution, provided that any
randomness used in finding the estimator is maintained as part of an expanded
Metropolis-Hastings parameter space. This could work for
$\rho_{\mathrm{full}}$. But now that we are working from
$\rho_{\mathrm{local}}$, this is impossible because the MCMC is already not
true Metropolis-Hastings with a single common parameter space. As discussed
above, the target distribution of the particle as a whole is different when
changing values at different loci, although there’s reason to hope that the
difference for nearby loci is small. Thus, $\rho_{\mathrm{sampled}}$ will
inevitably have a different stationary distribution from
$\rho_{\mathrm{local}}$. Nevertheless, in the algorithm below, in an attempt
to ensure that the stationary distribution changes as little as possible, I
expand the parameter space of $\rho_{\mathrm{sampled}}$ with a matrix
$\bm{\eta}^{s}$. This ensures that the same $H$ histories used when accepting
a value at a locus are also used when deciding whether to change that value
later.

We revise the algorithm from section 2 as follows:

  1. 1.

As before, assume we have $\hat{\tau}^{M}$.

  2. 2.

As before, for each particle $\bm{x}^{i}$, progress it to get a full particle
$\bm{\tilde{z}}^{i}\sim\mathsf{P}\bm{x}^{i}$.

  3. 3.

As before, find likelihood weights $w^{i}_{l}\equiv f(y_{l}|z^{i}_{l})$ and
forward densities $f^{j\rightarrow i}_{l}\equiv
f_{\mathsf{P}}(z^{i}_{l}|\bm{x}^{j})$ for all $i,j\in\\{1,\dots,M\\}$ and
$l\in\\{1,\dots,L\\}$.

  4. 4.

In parallel, for $k=1,\dots,M$, do the following:

    1. (a)

As before, sample $\bm{\iota}^{0}\in\\{1..M\\}^{L}$ to initialize the state of
the new proposal particle.

    2. (b)

In addition, initialize an $L\times H$ matrix $\bm{\eta}^{0}$ with entries in
$\\{1...M\\}$, sampled iid with probabilities

$P(\eta_{l,h}^{0}=i)=\frac{g(f^{i\rightarrow\iota^{0}_{l}}_{l})}{\Sigma_{j}g(f^{j\rightarrow\iota^{0}_{l}}_{l})},$

where $g$ is an arbitrary monotonically increasing function.

Each entry $\eta_{l,h}^{0}$ gives the index $i$ of one history $\bm{x}^{i}$
which we will later use to estimate the denominator of $\rho$. The
significance of $g$ will be explained in more detail below.

    3. (c)

As before, run a Metropolis-Hastings MCMC chain, for steps $s=1,...,S$,
updating $\bm{\eta}$ and $\bm{\iota}$ at each step:

      1. i.

As before, choose a spatial locus $\lambda(s)$ (aka $\lambda$) uniformly at
random.

      2. ii.

As before, sample a proposed replacement $\iota^{*}(s)$ (aka $\iota^{*}$) for
locus $\lambda$, with probability

$P(\iota^{*}=i)=w^{i}_{\lambda}/\sum_{j}w^{j}_{\lambda}.$

Also, define $\bm{\iota}^{**}$ as before.

      3. iii.

In addition, sample (iid) a set $(\eta^{*}_{1},...,\eta^{*}_{H})$ of histories
by which to judge this proposed replacement, where

$P(\eta^{*}_{h}=i)=\frac{g(f^{i\rightarrow\iota^{*}}_{\lambda})}{\Sigma_{j}g(f^{j\rightarrow\iota^{*}}_{\lambda})}.$

Define the matrix $\bm{\eta}^{**}$ by

$\eta^{**}_{l,h}:=\begin{dcases}\eta^{*}_{h}&\text{if }l=\lambda\\\
\eta^{s-1}_{l,h}&\text{if }l\neq\lambda\end{dcases}$

      4. iv.

Finally, define

$\rho_{\mathrm{sampled}}\equiv\frac{\Sigma_{h\in\\{1..H\\}}\frac{1}{g_{\lambda}(\eta^{*}_{h},\iota^{**})}[\Pi_{l\in\mathcal{B}_{r}(\lambda)}f^{\eta^{*}_{h}\rightarrow\iota^{**}}_{l}]}{\Sigma_{h\in\\{1..H\\}}\frac{1}{g_{\lambda}(\eta^{*}_{h},\iota^{s-1})}[\Pi_{l\in\mathcal{B}_{r}(\lambda)}f^{\eta^{*}_{h}\rightarrow\iota^{s-1}}_{l}]}\frac{\Sigma_{j=1}^{M}f^{j\rightarrow\iota^{s-1}}_{\lambda}}{\Sigma_{j=1}^{M}f^{j\rightarrow\iota^{**}}_{\lambda}},$
(4.1)

where

$g_{l}(i,j)\equiv\frac{g(f^{i\rightarrow j}_{l})}{\Sigma_{k}g(f^{k\rightarrow
j}_{l})}.$

As usual, we accept the proposed replacement with M-H probability
$1\wedge\rho_{\mathrm{sampled}}$. In case of acceptance, let
$\bm{\iota}^{s}:=\bm{\iota}^{**}$ and $\bm{\eta}^{s}:=\bm{\eta}^{**}$.
Otherwise, make no change, so that $\bm{\iota}^{s}:=\bm{\iota}^{s-1}$ and
$\bm{\eta}^{s}:=\bm{\eta}^{s-1}$.

    4. (d)

As before, set $z_{l}^{k}=\tilde{z}_{l}^{\iota^{S}_{l}}$.

In addition to showing that using unbiased estimators can conserve the target
distribution when the parameter space is expanded, Andrieu and Roberts also
discuss the case where the parameter space is not expanded. They show that
this case still has a stationary distribution, which converges to the original
target distribution as more samples are taken. With minor modifications, the
same proof applies to the MCMC using $\rho_{\mathrm{sampled}}$; the stationary
distribution of $\rho_{\mathrm{sampled}}$ is not the same as that of
$\rho_{\mathrm{local}}$, but converges to it as $H\rightarrow\infty$.

A few words on the function $g$. First, note that the sampling probabilities
for $\bm{\eta}^{s}$ are in principle arbitrary, so could be a function of the
full current vector $\bm{\iota}^{s}$. In the above discussion, however, they
are a function of only $f^{i\rightarrow\iota^{*}}_{l}$; this simplifies both
notation and computation. Second, $g$ should be chosen to be some non-
decreasing function of $f^{i\rightarrow\iota^{*}}_{\lambda}$, so that the
variance in
$\Pi_{l\in\mathcal{B}_{r}(\lambda)}f^{\eta^{*}_{h}\rightarrow\iota^{**}}_{l}$
is at least partially offset by that in $g$, increasing the efficiency of the
estimation process. Another way of saying this is that we should make it more
likely to sample plausible histories than implausible ones. Possible choices
for $g$ are discussed in Section 4.3 below. Whatever $g$ is chosen, the
denominator ${\Sigma_{j}g(f^{j\rightarrow\iota^{*}}_{\lambda})}$ can be
precalculated for each possible choice of $\iota^{*}$, meaning that this does
not meaningfully increase computing requirements per MCMC step.

How much computation does $\rho_{\mathrm{sampled}}$ save? The
$\rho_{\mathrm{local}}$ algorithm requires running $M$ different MCMC chains,
with each of $L$ loci going through $S$ steps, and at each step calculating a
$\rho$ using a sum over $M$ histories of a product over the up to $N$
locations in the relevant $\mathcal{B}_{r}(l)$. The total computation cost is
at least $O(M^{2}LNS)$. This is better than $O(M^{2}L^{2}S)$ that
$\rho_{\mathrm{full}}$ would have taken, but still somewhat burdensome. To get
$\rho_{\mathrm{sampled}}$, on the other hand, we only calculate the product of
locus likelihoods for an arbitrary number $H$ of histories rather than all $M$
of them. Thus, the total computation cost falls to $O(MHLNS)$; since $H$ and
$N$ are arbitrary constants that can be set independently from the full size
$M$ and $L$ respectively, this is a substantial improvement.

###  4.2 Discussion of proposal weights

In the above discussion, the proposal weights ${v}_{l}^{i}$, used by all
versions of $\rho$, are arbitrary. That is, any proposals with nonzero weights
could be used; the ${v}_{l}^{i}$ are accounted for out in the $\bar{F}$ term.

Ideally, these proposal weights should both be tuned for maximum efficiency of
the MCMC; that is, insofar as it does not substantially increase the
computational costs per step, to try to ensure that the variance of the
acceptance probability is as low as possible (approaching the Gibbs sampling
case where it’s uniformly 1) while maintaining disperse (high-variance)
proposal values for good ergodicity/mixing.

For ${v}_{l}^{i}$, that is similar to the idea of an optimal proposal
distribution, which is common in the particle filter literature.[17] The low-
dimensional bootstrap particle filter uses $(z|x^{1..M})$ as a proposal, then
reweights using $(y|z)$. In such a case, the idea of an ideal proposal
distribution is that if you could propose from $(z|x^{1..M},y)$, the
reweighting step would not be necessary.

Applying a similar idea to ${v}_{l}^{i}$, it becomes clear why I have set it
equal to $w_{l}^{i}$. Of course, including a factor of $w_{l}^{i}$ in
${v}_{l}^{i}$ helps these terms cancel and thus simplifies the calculation of
$\rho$. But if simplicity of calculating $\rho$ were the only consideration, I
could have set ${v}_{l}^{i}$ to $w_{l}^{i}{\Sigma_{j=1}^{M}f^{j\rightarrow
i}_{l}}$, so that the ratio
$\frac{\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s-1)}_{l}}{\Sigma_{j=1}^{M}f^{j\rightarrow\iota(s)*}_{l}}$
would cancel out too.

But setting ${v}_{l}^{i}\equiv w_{l}^{i}$ ensures that the proposal density
for $z_{l}$, conditional on $\bm{x}^{1..M}$ and $\bm{y}$, is
$(z_{l}|\bm{x}^{1..M},y_{l})$ — not too far from the ideal
$(z_{l}|\bm{x}^{1..M},\bm{y},\bm{\zeta}^{s-1}_{\mathcal{B}_{r}(l)\setminus
l})$ which would allow an acceptance probability of uniformly 1. That’s
because the density of $(z_{l}|\bm{x}^{1..M})$ is included implicitly through
the progression procedure, while that of $(z_{l}|y_{l})$ is handled explicitly
through $w_{l}^{i}$.

###  4.3 Refining the history sampling weights

What about the history sampling weights $g_{l}^{i\rightarrow j}$? As above,
these are arbitrary. Ideally, to minimize the variance of the acceptance
probability, they would approximate:

$\displaystyle g_{l}^{i\rightarrow
j}(\bm{\zeta}^{s-1}_{\mathcal{B}_{r}(l)\setminus l})$ $\displaystyle\propto
f_{\mathsf{P}}(\tilde{z}^{j}_{l}|\bm{x}^{i},\bm{\zeta}^{s-1}_{\mathcal{B}_{r}(l)})$
(4.2) $\displaystyle=\Pi_{\kappa\in\mathcal{B}_{r}(l)}f^{i\rightarrow
j}_{\kappa}$ (4.3)

… because this expression in the Horvitz-Thompson inverse sampling weight term
would cancel exactly with its relative contribution to the estimated sum of
procucts, so that the overall estimator would be governed solely by the sum of
weights term in the denominator.

It is computationally infeasible to calculate these quantities exactly for
each step of the MCMC, so I use $g_{l}^{i\rightarrow j}$. In the simulation
below, I’ve tested two formulas for these weights:

  1. 1.

$g_{\mathrm{uniform}}(x)=1$

  2. 2.

$g_{\mathrm{bentlog}}(x)=\frac{\log(f^{i\rightarrow
j}_{\lambda})-\log(\min(f^{i\rightarrow
j}_{\lambda})}{{\alpha}}+max(0,\log(x)-\log(\max(x)+{\beta})$, where $\min(x)$
and $\max(x)$ are the precalculated minimum and maximum values of
$f^{i\rightarrow j}_{\lambda}$ and ${\alpha},{\beta}$ are positive constants.

Both of these options are simply computationally-convenient first attempts;
though simulations show $g_{\mathrm{bentlog}}$ is an improvement over
$g_{\mathrm{uniform}}$, it is surely not optimal in this regard. In further
work, I will look into using proposal distributions that are conditional on
the current values at other loci, not just on the observations at the current
locus.

###  4.4 Theoretical limitations of the algorithms in this paper

Both the block particle filter and the Finkelstein particle filter proposed
here are intended to deal with the curse of dimensionality. However, both may
fail in cases where forward densities — that is, the relative probabilities of
given states at time $t$ conditional on the state at time $t-1$ — are
concentrated around particular values, and thus insufficiently ergodic.
Rebeschini and Van Handel’s error bounds rely on a strong ergodicity
assumption, bounding the forward density away from $0$ in a way that they
themselves acknowledge is unrealistic in real-world cases. In a separate
paper, they explore further the kind of problems that can arise when this
assumption does not apply, and the regime where that failure occurs in
practice.[15]

For the Finkelstein particle filter, I am not giving any formal proofs of
performance, but it is clear that if ergodicity is poor enough, my algorithm
will also break down. For example, suppose that the forward density from
history $\bm{x}^{j}$ to raw locus value $z^{i}_{l}$ is less than $\epsilon$ if
$i\neq j$, and otherwise greater than $100\epsilon$. In that case, the MCMC
will strongly tend to get stuck in states where all locus values come from the
same history. Metaphorically, Sally Finkelstein would be too picky, never
selecting nearby body parts that didn’t match perfectly. The result would then
reduce to the bootstrap particle filter, with more useless computational cost.

Do the real-world problems to which these algorithms are applied have enough
ergodicity for them to function? For these algorithms to be appropriate, we’d
need a situation with enough nonlinear effects that simple Kalman filters
don’t suffice; yet also one which still has plenty of new randomness at each
time step, such that even if the forward density is not actually strongly
ergodic, it is at least diffuse enough for these algorithms to work. In SLAM
(simultaneous location and mapping) models for robotics applications, such
situations arise. But in fluid dynamics models, chaotic dynamics are the rule.
Such models can be deterministic or nearly so, with highly concentrated
forward densities, yet still have interesting dynamics. Uncertainty in the
initial conditions is amplified at each time step, so even in a deterministic
model with new measurement tending to reduce the uncertainty at each time
step, overall uncertainty can remain in equilibrium.

Due to the “picky Sally” problem explained just above, the Finkelstein
algorithm as explained in current paper does not deal well with such
deterministic or nearly-deterministic models. However, in a follow-up paper, I
will offer a modification of this algorithm to address such models.

##  5 Numerical simulations

###  5.1 Setup

Filtering algorithms cannot be expected to precisely infer the underlying true
state of the hidden Markov model. Instead, the goal is merely to infer its
conditional distribution, and most of the interest of the problem lies in the
fact that this distribution remains non-degenerate; as we acquire more
information in order to narrow down the possible states, the state itself
evolves, so we never catch up.

Thus, we cannot simply follow the recipe of a more traditional simulation
study of a technique for parameter inference. In traditional parameter
inference, the object of interest is the true parameter value(s), which can be
arbitrarily chosen when running a simulation. Although inference algorithms
may yield an inferred distribution for the parameter(s), the interpretation of
this distribution as a confidence distribution (for frequentist methods) or a
credible distribution (for Bayesian methods) is in some sense not inherent to
the problem; for example, in the case of Bayesian methods, a credible interval
is only as valid as the priors that produce it. However, in this case, we are
not making arbitrary assumptions in order to get the best or the most robust
performance; the assumptions are given by the problem, and the aim is to
calculate a true probability distribution. In order to efficiently measure an
algorithm’s performance, we’d like a setting where the correct value of the
object of interest – not the point value, but the conditional distribution –
is known.

So, in order to run a simulation study, I fall back on a linear Gaussian
model, where the Kalman filter algorithm gives an analytically-correct
filtering distribution. Of course, given that such an analytic solution does
exist, one would never in practice use an inexact filtering algorithm such as
those discussed by this paper. However, the ability of our more-general
algorithm to roughly reproduce the results of a Kalman filter is, at the
least, encouraging.

The particular linear Gaussian model I use is a model based on a progression
matrix $P$, a novelty matrix $N$, and a measurement error covariance matrix
$E$:

$\displaystyle z$ $\displaystyle=Px+\delta$ (5.1) $\displaystyle y$
$\displaystyle=z+\epsilon$ (5.2) $\displaystyle\delta$
$\displaystyle\sim\mathcal{N}(0,N)$ (5.3) $\displaystyle\epsilon$
$\displaystyle\sim\mathcal{N}(0,E)$ (5.4)

$P$ is a tridiagonal matrix. $N$ and $E$ are diagonal matrices with periodic
structure, so that some loci are best learned about through their neighbors.
Specifically, $N$’s elements alternate between a higher and a lower variance,
each value occurring at every 2nd locus, while $E$ has a lower variance at
every 5th locus.:

$\displaystyle P$ $\displaystyle\equiv\begin{bmatrix}{{b}}&{{c}}&{}&{}&{0}\\\
{{a}}&{{b}}&{{c}}&{}&{}\\\ {}&{{a}}&{{b}}&\ddots&{}\\\
{}&{}&\ddots&\ddots&{{c}}\\\ {0}&{}&{}&{{a}}&{{b}}\\\ \end{bmatrix}$ (5.5)
$\displaystyle{a}$
$\displaystyle\equiv.4;{b}\equiv.35;{c}\equiv.05;{a}+{b}+{c}=.8<1$ (5.6)
$\displaystyle N$ $\displaystyle\equiv\begin{bmatrix}{1}&{}&{}&{}&{}&{0}\\\
{}&{{q}}&{}&{}&{}&{}\\\ {}&{}&{1}&{}&{}&{}\\\ {}&{}&{}&{{q}}&&{}\\\
{}&{}&{}&&\ddots&{}\\\ {0}&{}&{}&{}&{}&{{q}}\\\ \end{bmatrix}$ (5.7)
$\displaystyle{q}$ $\displaystyle=.25$ (5.8) $\displaystyle E$
$\displaystyle\equiv\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{2}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\dots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{4}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{5}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{6}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\dots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{d}\\\\{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{e}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{0}\\\\{2}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\{4}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\{5}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{e}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\{6}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}\\\\{d}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{0}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle{1}\\\$\hfil\kern
5.0pt\crcr}}}}\right]$}}\\\\{e}$ $\displaystyle=.16\\\$ (5.10)

The state was initialized at mean 0 and variance 5 independently at each
locus. The model was run for 10 time steps, and for the particle filtering
algorithms the outcome variables of 5 separate runs were averaged.

A number of parameters were tried for the algorithms, but a good set of
numbers for comparing different models was: $400$ particles for the
Finkelstein variants, $400^{2}=160000$ particles for the bootstrap particle
filter, and $400^{2}/5=32000$ particles for the block particle filter
algorithm. These numbers were chosen so that each algorithm would take roughly
comparable computing time; the only step that requires computing power that is
quadratic in the number of particles is precalculating the forward densities
$f^{j\rightarrow i}_{l}$ in the Finkelstein algorithm.

All results for Figures 4 and 5 are for a 30-dimensional model. Results for
both Finkelstein and block particle filter algorithms remained materially
similar as model dimension was varied from 30 to 90, demonstrating that the
Finkelstein and Frankenstein algorithms’ errors are roughly independent of
dimension, as expected.

The Finkelstein algorithm was used with $\rho_{\mathrm{sampled}}$, with 45
histories per location and two formulas for the history sampling probabilities
$g$ ($g_{\mathrm{uniform}}$, and $g_{\mathrm{bentlog}}$ with
${\alpha}={\beta}=5$). The neighborhood width was $r\equiv 1$, which is to say
that $\max|\mathcal{B}_{r}(l)|=3$. Similarly, the zone size for the block
particle filter algorithm was 3.

Note that I am not the first to simulate outcomes for the block particle
filter. Although Rebeschini and van Handel’s paper originally proposing it
relied on proofs rather than simulations, more recent papers have implemented
it and given results.[10] [7] [22] The results there are not directly
comparable with those given here due to different models used.

###  5.2 Results

To get an intuition for this situation, I will begin by showing the evolution
of a single run of the model. The top panel of figure 3 shows the evolution
over time of the true value ($\sum_{l=3}^{5}z_{l}$), the observed value
($\sum_{l=3}^{5}y_{l}$), and the filtering distribution as calculated by a
Kalman filter ($\bm{E}_{\pi}(\sum_{l=3}^{5}z_{l}|\bm{y}_{1},...,\bm{y}_{t}$),
for the sum of loci 3-5. The bottom panel shows the the mean of the estimated
filtering distribution for each of the four algorithms I tested — Kalman
filter (analytically correct), bootstrap particle filter, block particle
filter, and Finkelstein particle filter.

Figure 3: Time series from a single run of each algorithm for 10 time steps,
in a model with 90 dimensions. Parameters for each algorithm are given in
text.

In the upper panel of Figure 3, one can see that the observations vary
relatively widely around the truth, while the Kalman filter mean follows those
observations with more conservative moves, thus staying closer to the true
value. In the lower panel, one can see that the bootstrap particle filter
falls to the curse of dimensionality; while the block and Finkelstein particle
filters both manage to approximate the correct Kalman filter distribution
relatively well.

To compare outcomes of the two working algorithms in greater depth, Figure 4
shows the average squared error per locus: that is, the squared difference
between the estimated distribution mean and the correct mean as given by the
Kalman filter, conditional on a single fixed series of observations
$(y_{1},...,y_{10})$. Note that we are measuring error relative to the mean of
the Kalman filter rather than to $z_{l}$; this is because the Kalman filter
result is the ideal filtering distribution that we are trying to capture here.
Though it’s not visible in these graphs, both algorithms do a relatively good
job of reproducing the variance of the filtering distribution; this is within
3% of the true value for both algorithms.

Figure 4: Breakdown of average squared error per locus.

Figure 4 makes two things clear. First, there is a bias/variance tradeoff
between the block particle filter and the Finkelstein particle filter; with
comparable run times, the block filter has a higher bias but almost no
variance in its distribution error. Second, as expected, the performance of
the block particle filter differs for loci that are central to their
neighborhood $Z_{j}$, as opposed to loci that are on the border of their
neighborhood. (The small apparent difference in performance of the Finkelstein
algorithm between the two kinds of loci is largely an artifact of the specific
realization of $(y_{1},...,y_{10})$ that was used to generate this graph. The
Finkelstein algorithm does not use fixed neighborhoods $\\{Z_{j}\\}$, so the
distinction between central and peripheral is simply does not apply, except
for the first and last loci overall.)

The error of the Finkelstein algorithm, like that of the block particle
filter, appears to remain stable over time, as seen in Figure 5. This shows
the time evolution of the KL divergence between the true filtering
distribution, as calculated using the Kalman filter, and a Gaussian with mean
vector and covariance matrix inferred from a Finkelstein particle filter. It
appears that uniform sampling ($g_{\mathrm{uniform}}$) can occasionally be
unsuccessful in sampling good histories, as reflected by the spikes in that
line; log sampling ($g_{\mathrm{bentlog}}$) had superior stability.

Figure 5: Stability of KL divergence across time steps.

##  6 Conclusion

I have introduced the novel Finkelstein particle filtering algorithm for
estimating the filtering distribution of models with high dimensionality due
to large spatial extent. In such models, the simple bootstrap particle
filtering algorithm is unusable. But, as with the previously-proposed block
particle filter, my algorithm relies on the locality of dynamics to resolve
this problem, focusing on a small area at a time.

Using simulations, I have showed that the error of means of my algorithm has
lower bias but higher variance than the block particle filter, given
comparable parameters. All in all, the total squared error of means of the
Finkelstein algorithm is more homogeneous across loci than that of the block
algorithm; lower for loci peripheral to a neighborhood in the block particle
filter, but higher for those which are central. I also give empirical evidence
that the error of this algorithm is stable over time, making it a candidate
for online data assimilation tasks.

It is commonplace to prefer variance over bias when such a tradeoff is
possible, because this allows improving precision with additional computing
power by independent reruns of the algorithm. That improved precision would
certainly be possible in this case with the Finkelstein algorithm. This
picture is slightly complicated by the fact that such computing power might
enable better results from the block particle filter by increasing the
neighborhood size. But there are several problems with just increasing
neighborhood size. Above all, computing power (that is, number of particles)
needed could be up to exponential in neighborhood size, while it’s just
quadratic in number of Finkelstein particles or linear in independent
Finkelstein runs. Second, unlike number of particles, neighborhood size comes
in sizeable discrete intervals; it may not be possible to effectively use a
small additional amount of computing power. And finally, to reduce the bias of
the block particle filter, neighborhood size must be increased up-front, while
the variance Finkelstein particle filter can in be reduced by independent runs
(perhaps even by different scientists).

Thus, I believe that the Finkelstein particle filter algorithm offers
meaningful advantages over prior proposals. In future work, I will extend this
to cover chaotic dynamics in a deterministic or quasi-deterministic model.

## References

  * [1] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso and A.. Stuart  “Importance Sampling: Intrinsic Dimension and Computational Cost”  In _Statistical Science_ 32.3, 2017, pp. 405–431  DOI: 10.1214/17-STS611
  * [2] Christophe Andrieu and Gareth O. Roberts  “The Pseudo-Marginal Approach for Efficient Monte Carlo Computations”  In _The Annals of Statistics_ 37.2, 2009, pp. 697–725  DOI: 10.1214/07-AOS574
  * [3] A. Apte, C…. Jones, A.. Stuart and J. Voss  “Data Assimilation: Mathematical and Statistical Perspectives”  In _International Journal for Numerical Methods in Fluids_ 56.8, 2008, pp. 1033–1046  DOI: 10.1002/fld.1698
  * [4] Thomas Bengtsson, Peter Bickel and Bo Li  “Curse-of-Dimensionality Revisited: Collapse of the Particle Filter in Very Large Scale Systems”  In _Institute of Mathematical Statistics Collections_ Beachwood, Ohio, USA: Institute of Mathematical Statistics, 2008, pp. 316–334  DOI: 10.1214/193940307000000518
  * [5] Peter Bickel, Bo Li and Thomas Bengtsson  “Sharp Failure Rates for the Bootstrap Particle Filter in High Dimensions”  Institute of Mathematical Statistics, 2008  DOI: 10.1214/074921708000000228
  * [6] Olivier Cappe, Eric Moulines and Tobias Ryden  “Inference in Hidden Markov Models”, Springer Series in Statistics  Springer, 2009
  * [7] Alban Farchi and Marc Bocquet  “Review Article: Comparison of Local Particle Filters and New Implementations”  In _Nonlinear Processes in Geophysics_ 25.4, 2018, pp. 765–807  DOI: https://doi.org/10.5194/npg-25-765-2018
  * [8] Walter R Gilks and Carlo Berzuini  “Following a Moving Target – Monte Carlo Inference for Dynamic Bayesian Models”  In _Journal of the Royal Statistical Society_ 63.1, Statistical Methodology, 2001, pp. 127–146  DOI: 1369-7412/01/63127
  * [9] Simon Godsill and Tim Clapp  “Improvement Strategies for Monte Carlo Particle Filters”  In _Sequential Monte Carlo Methods in Practice_ Springer-Verlag, 2000, pp. 139–158
  * [10] M. Morzfeld, D. Hodyss and J. Poterjoy  “Variational Particle Smoothers and Their Localization”  In _Quarterly Journal of the Royal Meteorological Society_ 144.712, 2018, pp. 806–825  DOI: 10.1002/qj.3256
  * [11] Matthias Morzfeld, Xin T. Tong and Youssef M. Marzouk  “Localization for MCMC: Sampling High-Dimensional Posterior Distributions with Local Structure”, 2017  arXiv:1710.07747 [math, stat]
  * [12] Jonathan Poterjoy  “A Localized Particle Filter for High-Dimensional Nonlinear Systems”  In _Monthly Weather Review_ 144.1, 2015, pp. 59–76  DOI: 10.1175/MWR-D-15-0163.1
  * [13] Jonathan Poterjoy, Ryan A. Sobash and Jeffrey L. Anderson  “Convective-Scale Data Assimilation for the Weather Research and Forecasting Model Using the Local Particle Filter”  In _Monthly Weather Review_ 145.5, 2017, pp. 1897–1918  DOI: 10.1175/MWR-D-16-0298.1
  * [14] Patrick Rebeschini and Ramon Handel  “Can Local Particle Filters Beat the Curse of Dimensionality?” 00070 MR: MR3375889  In _The Annals of Applied Probability_ 25.5, 2015, pp. 2809–2866  DOI: 10.1214/14-AAP1061
  * [15] Patrick Rebeschini and Ramon Handel  “Phase Transitions in Nonlinear Filtering”  In _Electronic Journal of Probability_ 20, 2015  DOI: 10.1214/EJP.v20-3281
  * [16] François Septier and Gareth W. Peters  “An Overview of Recent Advances in Monte-Carlo Methods for Bayesian Filtering in High-Dimensional Spaces”  In _Theoretical Aspects of Spatial-Temporal Modeling_ Tokyo: Springer Japan, 2015, pp. 31–61  DOI: 10.1007/978-4-431-55336-6˙2
  * [17] Chris Snyder  “Particle Filters, the “Optimal” Proposal and High-Dimensional Systems”, ECMWF Seminar on Data assimilation for atmosphere and ocean, 2011, pp. 10
  * [18] Chris Snyder, Thomas Bengtsson, Peter Bickel and Jeff Anderson  “Obstacles to High-Dimensional Particle Filtering”  In _Monthly Weather Review_ 136.12, 2008, pp. 4629–4640  DOI: 10.1175/2008MWR2529.1
  * [19] Cyrill Stachniss, John J. Leonard and Sebastian Thrun  “Simultaneous Localization and Mapping”  In _Springer Handbook of Robotics_ , Springer Handbooks  Cham: Springer International Publishing, 2016, pp. 1153–1176  DOI: 10.1007/978-3-319-32552-1˙46
  * [20] Peter Jan Van Leeuwen, Yuan Cheng and Sebastian Reich  “Nonlinear Data Assimilation”, Frontiers in applied dynamical systems reviews and tutorials 2  Cham: Springer, 2015
  * [21] Peter Jan Leeuwen  “Particle Filtering in Geophysical Systems”  In _Monthly Weather Review_ 137.12, 2009, pp. 4089–4114  DOI: 10.1175/2009MWR2835.1
  * [22] Xu Yaxian  “Sequential Monte Carlo Algorithms Fir High-Dimensional Filtering and Smoothing”, 2018  URL: http://scholarbank.nus.sg/handle/10635/146933