A critical analysis of resampling strategies for the regularized particle filter

TL;DR

This paper critically examines resampling strategies in regularized particle filters, revealing that resampling can hinder convergence and proposing methods like kernel bandwidth modulation to improve performance, supported by numerical experiments.

Contribution

It provides a detailed analysis of resampling effects in regularized particle filters and introduces techniques to mitigate convergence issues, outperforming traditional bootstrap methods.

Findings

Resampling can prevent convergence to the true posterior.

Kernel bandwidth modulation improves filter performance.

Regularized particle filters outperform bootstrap filters with proposed methods.

Abstract

We analyze the performance of different resampling strategies for the regularized particle filter regarding parameter estimation. We show in particular, building on analytical insight obtained in the linear Gaussian case, that resampling systematically can prevent the filtered density from converging towards the true posterior distribution. We discuss several means to overcome this limitation, including kernel bandwidth modulation, and provide evidence that the resulting particle filter clearly outperforms traditional bootstrap particle filters. Our results are supported by numerical simulations on a linear textbook example, the logistic map and a non-linear plant growth model.

Figures

Tables

Table 1. Table 1: Comparison of particle filter predictions for different resampling strategies at n = 100 𝑛 100 n=100 averaged over 10 runs.

	Benchmark	RPF (ESS${}_{\text{crit}}$ = 0.5)	RPF (${\alpha }_n={\left(n+{\alpha }_h^{-1}\right)}^{-1}$)	SIR (ESS${}_{\text{crit}}$ = 0.5)	SIS (ESS${}_{\text{crit}}$ = 0)
${\mu }_{\text{RUE}}$	3.521	3.479 (0.0096)	3.505 (0.0291)	3.554 (0.0431)	3.503 (0.0569)
${\sigma }_{\text{RUE}}$	0.0325	0.0487 (0.0041)	0.0341 (0.0083)	0.0126 (0.0161)	0.0060 (0.0049)
${\mu }_{\gamma }$	0.632	0.639 (0.0026)	0.636 (0.0090)	0.623 (0.0095)	0.634 (0.0089)
${\sigma }_{\gamma }$	0.0070	0.0146 (0.0016)	0.0074 (0.0032)	0.0024 (0.0030)	0.0007 (0.0005)
${\mu }_{{\mu }_a}$	542.6	537.1 (2.37)	540.4 (7.55)	547.2 (18.08)	540.4 (10.50)
${\sigma }_{{\mu }_a}$	8.297	14.43 (1.69)	9.996 (1.74)	4.764 (5.89)	2.477 (3.09)

Equations

$$\begin{array}{rl}& {\int }_{\mathcal{X}}\varphi ({\mathbf{x}}_n){\pi }_n(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n})\\ & =\frac{1}{Z_n}{\int }_{{\mathcal{X}}^{n+1}}\varphi ({\mathbf{x}}_n){\pi }_0(d{\mathbf{x}}_0)\prod _{j=1}^n{g}_j({\mathbf{y}}_j|{\mathbf{x}}_j)f_j(d{\mathbf{x}}_j|{\mathbf{x}}_{j-1}),\end{array}$$

$$Z_n={\int }_{{\mathcal{X}}^{n+1}}{\pi }_0(d{\mathbf{x}}_0)\prod _{j=1}^n{g}_j({\mathbf{y}}_j|{\mathbf{x}}_j)f_j(d{\mathbf{x}}_j|{\mathbf{x}}_{j-1}).$$

$${\pi }_{1|0}^N(d\mathbf{x})=\frac{1}{N}\sum _{i=1}^N{\delta }_{{\mathbf{x}}_{1|0}^{(i)}}(d\mathbf{x})$$

$${\pi }_1^N(d{\mathbf{x}}_1|{\mathbf{y}}_1)=\sum _{i=1}^N{w}_1^{(i)}{\delta }_{{\mathbf{x}}_{1|0}^{(i)}}(d{\mathbf{x}}_1)$$

$$w_1^{(i)}=\frac{g_1({\mathbf{y}}_1|{\mathbf{x}}_{1|0}^{(i)})}{{\sum }_{j=1}^N{g}_1({\mathbf{y}}_1|{\mathbf{x}}_{1|0}^{(j)})}$$

$${\int }_{\mathcal{X}}\varphi ({\mathbf{x}}_1){\pi }_1^N(d{\mathbf{x}}_1|{\mathbf{y}}_1)=\sum _{i=1}^N{w}_1^{(i)}\varphi ({\mathbf{x}}_{1|0}^{(i)}),$$

$$\underset{N\to \infty }{\lim }{\int }_{\mathcal{X}}\varphi ({\mathbf{x}}_1){\pi }_1^N(d{\mathbf{x}}_1|{\mathbf{y}}_1)={\int }_{\mathcal{X}}\varphi ({\mathbf{x}}_1){\pi }_1(d{\mathbf{x}}_1|{\mathbf{y}}_1).$$

$$Z_n^N=\frac{1}{N}\sum _{i=1}^N\prod _{j=1}^n{g}_j({\mathbf{y}}_j|{\mathbf{x}}_{j|j-1}^{(i)}).$$

$$\frac{1}{N^2}\sum _{i=1}^N{\left(\prod _{j=1}^n{g}_j({\mathbf{y}}_j|{\mathbf{x}}_{j|j-1}^{(i)})-Z_n^N\right)}^2=\frac{(Z_n^N{)}^2}{{\text{ESS}}_n}\left(1-\frac{{\text{ESS}}_n}{N}\right),$$

$${\text{ESS}}_n=\frac{1}{{\sum }_{i=1}^N(w_n^{(i)}{)}^2}$$

$${\tilde{\pi }}_n^N(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n})=\frac{1}{N}\sum _{i=1}^N{\delta }_{{\mathbf{x}}_n^{(i)}}(d{\mathbf{x}}_n).$$

$$\frac{1}{\lambda ({\mathcal{A}}_{ϵ})}{\pi }_n^N({\mathcal{A}}_{ϵ})=\frac{1}{\lambda ({\mathcal{A}}_{ϵ})}\sum _{i=1}^N{w}_n^{(i)}{\delta }_{{\mathbf{x}}_{n|n-1}^{(i)}}({\mathcal{A}}_{ϵ}),$$

$$p_n^N({\mathbf{x}}_n|{\mathbf{y}}_{1:n})=\sum _{i=1}^N{w}_n^{(i)}{\mathcal{K}}_{h_n}({\mathbf{x}}_n-{\mathbf{x}}_{n|n-1}^{(i)}),$$

$${\alpha }_h={\left(\frac{4}{N(d_x+2)}\right)}^{\frac{2}{d_x+4}}.$$

$${\Sigma }_n^N=\frac{N}{N-1}\sum _{i=1}^N{w}_n^{(i)}({\mathbf{x}}_{n|n-1}^{(i)}-{\mu }_n^N{)}^2$$

$$\mathbb{E}\left[||p_n^N-p_n|{|}_{L^2}\right]\underset{N\to \infty }{=}\mathcal{O}\left(N^{-\frac{2}{d_x+4}}\right).$$

$$\begin{array}{rl}\text{RMSE}({\hat{\mathbf{x}}}_n,{\mathbf{x}}_n)& =\mathbb{E}{\left[||{\hat{\mathbf{x}}}_n-{\mathbf{x}}_n|{|}^2\right]}^{1/2}\\ & =\sqrt{||{\mu }_n-{\mathbf{x}}_n|{|}^2+\text{Tr}[{\Sigma }_n]},\end{array}$$

$$\left|\begin{array}[]{l}{\bf x}_{n}=A_{n}{\bf x}_{n-1}+\epsilon_{n}\vspace{0.1cm}\\ {\bf y}_{n}=B_{n}{\bf x}_{n}+\eta_{n}\end{array}\right.$$

$$\left|\begin{array}[]{l}\Sigma_{n}^{-1}=\Sigma_{n|n-1}^{-1}+B_{n}^{T}R_{n}^{-1}B_{n}\vspace{0.1cm}\\ \Sigma_{n}^{-1}{\bf\mu}_{n}=\Sigma_{n|n-1}^{-1}{\bf\mu}_{n|n-1}+B_{n}^{T}R_{n}^{-1}{\bf y}_{n}\end{array}\;,\right.$$

$$\begin{array}{rl}{\pi }_{n|n-1}(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n-1})={\int }_{\mathcal{X}}& f_n(d{\mathbf{x}}_n|{\mathbf{x}}_{n-1})\\ & \times {\pi }_{n-1}(d{\mathbf{x}}_{n-1}|{\mathbf{y}}_{1:n-1}),\end{array}$$

$${\pi }_n(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n})=\frac{g_n({\mathbf{y}}_n|{\mathbf{x}}_n){\pi }_{n|n-1}(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n-1})}{{\int }_{\mathcal{X}}{g}_n({\mathbf{y}}_n|{\mathbf{x}}_n){\pi }_{n|n-1}(d{\mathbf{x}}_n|{\mathbf{y}}_{1:n-1})},$$

$$\left|\begin{array}[]{l}\Sigma_{n}=R(R+n\Sigma_{0})^{-1}\Sigma_{0}\vspace{0.1cm}\\ {\bf\mu}_{n}=\Sigma_{n}\Sigma_{0}^{-1}{\bf\mu}_{0}+(\mathbb{I}-\Sigma_{n}\Sigma_{0}^{-1})\langle{\bf y}\rangle_{n}\end{array}\right.$$

$$\begin{array}{rl}& p_n^N({\mathbf{x}}_n|{\mathbf{y}}_{1:n})=\frac{\frac{1}{N}{\sum }_{i=1}^N{g}_n({\mathbf{y}}_n|{\mathbf{x}}_{n|n-1}^{(i)}){\mathcal{K}}_{h_n}({\mathbf{x}}_n-{\mathbf{x}}_{n|n-1}^{(i)})}{\frac{1}{N}{\sum }_{j=1}^N{g}_n({\mathbf{y}}_n|{\mathbf{x}}_{n|n-1}^{(j)})}\\ & \underset{N\to \infty }{\approx }\frac{{\int }_{\mathcal{X}}\lambda (d{\mathbf{x}}_n^{\prime })p_{n|n-1}^N({\mathbf{x}}_n^{\prime }|{\mathbf{y}}_{1:n-1})g_n({\mathbf{y}}_n|{\mathbf{x}}_n^{\prime }){\mathcal{K}}_{h_n}({\mathbf{x}}_n-{\mathbf{x}}_n^{\prime })}{{\int }_{\mathcal{X}}\lambda (d{\mathbf{x}}_n^{\prime })p_{n|n-1}^N({\mathbf{x}}_n^{\prime }|{\mathbf{y}}_{1:n-1})g_n({\mathbf{y}}_n|{\mathbf{x}}_n^{\prime })},\end{array}$$

$$\left(({\Sigma }_{n|n-1}^N{)}^{-1}+B_n^T{R}_n^{-1}{B}_n\right){\mu }_n^N=({\Sigma }_{n|n-1}^N{)}^{-1}{\mu }_{n|n-1}^N+B_n^T{R}_n^{-1}{\mathbf{y}}_n$$

$$\begin{array}{rl}{\Sigma }_n^N& =h_n^2+{\left(({\Sigma }_{n|n-1}^N{)}^{-1}+B_n^T{R}_n^{-1}{B}_n\right)}^{-1}\\ & =(1+{\alpha }_h){\left(({\Sigma }_{n|n-1}^N{)}^{-1}+B_n^T{R}_n^{-1}{B}_n\right)}^{-1},\end{array}$$

$$\begin{array}[]{l}\left(\Sigma_{n-1}^{-1}+R^{-1}\right){\bf\mu}_{n}=\Sigma_{n-1}^{-1}{\bf\mu}_{n-1}+R^{-1}{\bf y}_{n}\vspace{0.1cm}\\ \Sigma_{n}=(1+\alpha_{n})R\left(R+\Sigma_{n-1}\right)^{-1}\Sigma_{n-1}\end{array}$$

$$\left|\begin{array}[]{l}\Sigma_{n}=\prod_{j=1}^{n}(1+\alpha_{j})R\left(R+\sum_{j=1}^{n}\prod_{k=0}^{j-1}(1+\alpha_{k})\Sigma_{0}\right)^{-1}\Sigma_{0}\vspace{0.1cm}\\ \begin{split}&{\bf\mu}_{n}=\frac{1}{\prod_{j=1}^{n}(1+\alpha_{j})}\Sigma_{n}\Sigma_{0}^{-1}{\bf\mu}_{0}\\ &+\left(\mathbb{I}-\frac{1}{\prod_{j=1}^{n}(1+\alpha_{j})}\Sigma_{n}\Sigma_{0}^{-1}\right)\frac{\sum_{j=1}^{n}\prod_{k=0}^{j-1}(1+\alpha_{k}){\bf y}_{j}}{\sum_{j=1}^{n}\prod_{k=0}^{j-1}(1+\alpha_{k})}\end{split}\end{array}\right.$$

$$\left|\begin{array}[]{l}\Sigma_{n}^{N}=(1+\alpha_{h})^{n}R\left(R+\sum_{j=1}^{n}(1+\alpha_{h})^{j-1}\Sigma_{0}\right)^{-1}\Sigma_{0}\vspace{0.1cm}\\ \begin{split}{\bf\mu}_{n}^{N}&=\frac{1}{(1+\alpha_{h})^{n}}\Sigma_{n}^{N}\Sigma_{0}^{-1}{\bf\mu}_{0}\\ &+\left(\mathbb{I}-\frac{1}{(1+\alpha_{h})^{n}}\Sigma_{n}^{N}\Sigma_{0}^{-1}\right)\frac{\sum_{j=1}^{n}(1+\alpha_{h})^{j-1}{\bf y}_{j}}{\sum_{j=1}^{n}(1+\alpha_{h})^{j-1}}\end{split}\end{array}\right.$$

$$\underset{n\to \infty }{\lim }{\Sigma }_n^N={\alpha }_{h}R(1+\mathcal{O}(e^{-n{\alpha }_h})).$$

$${\mu }_n^N\underset{n\to \infty }{\approx }\frac{{\sum }_{j=1}^n(1+{\alpha }_h{)}^{j-1}{\mathbf{y}}_j}{{\sum }_{j=1}^n(1+{\alpha }_h{)}^{j-1}}\left(1+\mathcal{O}(e^{-n{\alpha }_h})\right)$$

Full text

††journal: Statistics and Computing

1 Introduction

Sequential Monte Carlo (SMC) methods DelMoral04 ; SMC are powerful algorithms which allow performing Bayesian inference in dynamic systems. These methods, also known under the heading of particle filters in the context of dynamic systems, have become extremely popular way beyond the borders of the statistics community and are nowadays widely used by engineers and scientists from other fields faced with non-trivial inference tasks which naturally abound with the emergence of data science. The purpose of this paper is to provide a critical analysis of the performance of regularized particle filter (RPF) algorithms LeGland98 ; Musso12 which represent an interesting variation of traditional particle filters. They rely on kernel-based Parzen62 posterior density estimation at each filtering step, with particle resampling performed from the reconstructed density. In many situations, this method was proven to help prevent sample impoverishment SMC and thereby enhance the robustness of particle filtering. Our main achievement will be to demonstrate in this paper that the approach which consists in resampling systematically can lead to seriously over-estimate the posterior mean squared error, motivating us to exhibit alternative resampling strategies which allow fully resolving this issue.

Before going any further, let us introduce the general framework we shall be working with in the rest of this paper. To begin with, we shall assume that the processes from which we seek to infer some information belong to the family of hidden Markov models (HMM) HMM , which describe stochastic processes characterizing the evolution in discrete time $n$ of two sets of random variables $\{{\bf x}_{n}\}_{n\geq 0}$ and $\{{\bf y}_{n}\}_{n\geq 1}$, respectively defined on some space ${\cal X}\subseteq\mathbb{R}^{d_{x}}$ and ${\cal Y}\subseteq\mathbb{R}^{d_{y}}$. The former are hidden random variables referred to as states, while the latter are noisy observations of the former. Starting from a prior probability measure $\pi_{0}(d{\bf x}_{0})$ on the initial state, hidden variables follow a Markovian evolution characterized by a sequence of Markov kernels $f_{n}(d{\bf x}_{n}| {\bf x}_{n-1})$ also named transition functions. Observations are considered conditionally independent and described by a conditional probability density function $g_{n}({\bf y}_{n}| {\bf x}_{n})$, the so-called measurement function, with respect to the Lebesgue measure $\lambda$. Parameters $\{\Theta\}$, with a compact support set, can easily be accounted for in that framework by augmenting the size of the state space, as can the presence of environmental variables $\{{\bf u}_{n}\}_{n\geq 0}$ defined on some space ${\cal U}\subseteq\mathbb{R}^{d_{u}}$.

A problem of broad interest in this context is to determine the posterior probability measure $\pi_{k}(d{\bf x}_{k}| {\bf y}_{1:n})$ given a sequence of observations ${\bf y}_{1:n}=\{{\bf y}_{1},...,{\bf y}_{n}\}$. If $k<n$, this is called a smoothing problem, while if $k>n$ this is a prediction problem. In the following we shall be interested in the situation where $k=n$, which is called the filtering problem. Thus, given any Borel set ${\cal A}\subset{\cal X}$, $\pi_{n}({\cal A})=\int_{\cal A}\pi_{n}(d{\bf x}_{n}| {\bf y}_{1:n})$ shall be the posterior probability of event ${\cal A}$. More generally, for any test function $\phi$ defined on ${\cal X}$, applying Bayes rule recursively allows expressing the integral of $\phi$ with respect to the posterior probability measure as a function of past measurement and transition functions HMM

[TABLE]

introducing the marginal likelihood

[TABLE]

The purpose of SMC methods Liu98 ; DelMoral06 ; SMC is to provide an estimation of such probability measures in the general case of nonlinear state-space models, where analytical attempts are essentially hopeless. The idea of such methods is to represent the probability distributions by particles, i.e. independent identically distributed samples of the distributions, and to propagate and update them using the HMM equations, hence the equivalent nomenclature of particle filter methods. The estimation of parameters in this context has received acute attention from the community these last years (see Kantas15 for a review). One of the main challenges is to circumvent the celebrated particle degeneracy problem Gordon93 and to ensure that parameter space is sufficiently explored to enable accurate estimation. Several methods have been proposed to address this challenge, among which iterated filtering Ionides11 and particle Markov chain Monte Carlo Andrieu10 in the specific context of joint parameter / hidden state estimation. They constitute batch or offline methods, however, in that they are unable to take into account new observations dynamically, as opposed to online algorithms such as SMC2 Chopin13 ; Jacob15 .

A different idea, introduced some time ago in West93 and discussed in more detail in Liu10 ; Musso12 , is to allow freedom in parameter space by regularizing the posterior density using a mixture. In this paper, we focus on the RPF introduced in Ref. LeGland98 and further developed and discussed in Refs Musso12 ; Campillo09 ; ChenPhD ; Gauchi13 . We shall highlight some shortcomings of this approach when resampling is performed systematically (as routinely done), which can be made explicit by computing the posterior mean squared error. In particular, building upon some well-known results obtained for linear Gaussian state-space models, we prove that the mean squared error is lower bounded by a positive constant when the number of particles is finite, independently of the number of observations. We shall provide a detailed understanding of this phenomenon on a very simple linear textbook example, offering evidence that it is related to the particle filter losing memory of past observations. We will then show that long-ranged memory can be restored by appropriately modulating the kernel bandwidth and that the efficiency of this suitably modified particle filter clearly outperforms that of the original bootstrap particle filter Gordon93 . The rest of the paper is organized as follows. In section 2, we recall basic principles of SMC methods and introduce the bootstrap particle filter and the RPF, before illustrating what can potentially go wrong when resampling is performed on a systematic basis. We then provide analytical insight regarding the effect of different resampling strategies on the behavior of the RPF in section 3, using a very simple linear textbook example for illustrative purposes. That our results qualitatively hold for more general, non-linear, models is discussed in section 4. Conclusions are drawn in section 5, while a series of technical details are gathered in appendices.

2 Sequential Monte Carlo methods

Let us begin by recalling the basic principles of particle filter methods. Starting from the initial prior probability measure from which $N$ particles are sampled, ${\bf x}_{0}^{(i)}\sim\pi_{0}(d{\bf x}_{0})$, particles are propagated until the next time step using a Markov kernel generally chosen Gordon93 ; Kong94 as ${\bf x}_{1|0}^{(i)}\sim f_{1}(d{\bf x}_{1}|{\bf x}_{0}^{(i)})$ in nonlinear state space models for which $f_{1}(d{\bf x}_{1}|{\bf x}_{0}^{(i)},{\bf y}_{1})$ is not readily available Liu98 . The resulting empirical estimator for the predicted probability measure can be expressed as

[TABLE]

with $\delta$ the Dirac measure. The estimator for the updated posterior probability measure is obtained by weighting each particle according to the observation ${\bf y}_{1}$:

[TABLE]

with weights

[TABLE]

suitably normalized to ensure that total weight is conserved. Given any bounded function $\phi$ defined on ${\cal X}$, its integral with respect to the empirical posterior probability measure is then given by

[TABLE]

and the estimator can be shown to converge asymptotically, for an infinite number of particles, towards the true posterior probability measure given by Bayes rule Chopin04 :

[TABLE]

Practically speaking, however, one always works with a finite number of particles and life is not as simple. A well-known problem which can occur with the above estimation procedure is the so-called weight degeneracy issue SMC . If the likelihood is relatively peaked around the observation value, or in other words if the observation noise is rather small, particle weights will differ by orders of magnitude depending on the distance between the particle and the observation. Actually, this will happen almost regardless of the magnitude of the observation noise for a sufficient number of observations, since iterating the above procedure yields $w_{n}^{(i)}\propto\prod_{j=1}^{n}g_{j}({\bf y}_{j}|{\bf x}_{j|j-1}^{(i)})$, using the recursive formulation ${\bf x}_{j|j-1}^{(i)}\sim f_{j}(d{\bf x}_{j}|{\bf x}_{j-1}^{(i)})$. This implies that, eventually, most of the weight will be carried by a small number of particles, with a dramatic effect on the precision of the resulting estimation as can be seen by looking at the unbiased estimator of the marginal likelihood at step $n$:

[TABLE]

Its variance can be estimated using the standard empirical estimator

[TABLE]

where we introduced the effective sample size

[TABLE]

which is a measure of the effective number of particles Kong94 ; Carpenter99 : if all particles carry the same weight, $\text{ESS}_{n}=N$, while if a single particle carries all the weight, $\text{ESS}_{n}=1$. This result therefore illustrates that the accuracy of the estimation by the empirical measure relies on a relatively homogeneneous distribution of particle weights.

A natural strategy to counter the effect of particle degeneracy is to ’unweight’ particles by resampling from the empirical posterior probability measure, ${\bf x}_{n}^{(i)}\sim\pi_{n}^{N}(d{\bf x}_{n}|{\bf y}_{1:n})$, such that the new approximation to the posterior is

[TABLE]

This is the idea behind the bootstrap particle filter Gordon93 , also known as the sequential importance resampling (SIR) algorithm. Its main drawback is that resampling according to a probability measure such as Eq. (4) essentially duplicates particles (particles are selected following a multinomial law), leading to so-called sample impoverishment SMC .

2.1 Regularized particle filter

In Ref. LeGland98 , LeGland, Musso and Oujdane proposed to address these shortcomings by regularizing the probability measure using a Parzen-Rosenblatt Parzen62 kernel density estimator. The idea is to replace the Dirac measure by a measure able to smoothly interpolate between all weighted particles. In other words, observing that the coarse-grained density on Borel set ${\cal A}_{\epsilon}\subset{\cal X}$

[TABLE]

where $\lambda({\cal A}_{\epsilon})=\int_{{\cal A}_{\epsilon}}\lambda(d{\bf x}_{n})$, the regularization procedure consists in replacing $\delta_{{\bf x}}({\cal A}_{\epsilon})/\lambda({\cal A}_{\epsilon})$ as $\lambda({\cal A}_{\epsilon})\to 0$ by a kernel density function with respect to the Lebesgue measure Tsybakov , thanks to which the regularized estimator for the posterior probability density function reads

[TABLE]

with

• $\bullet$

$w_{n}^{(i)}=g_{n}({\bf y}_{n}| {\bf x}_{n|n-1}^{(i)})/\sum_{j=1}^{N}g_{n}({\bf y}_{n}| {\bf x}_{n|n-1}^{(j)})$

• $\bullet$

${\bf x}_{n|n-1}^{(i)}\sim f_{n}(d{\bf x}_{n}|{\bf x}_{n-1}^{(i)})$

• $\bullet$

${\bf x}_{n-1}^{(i)}\sim p_{n-1}^{N}({\bf x}_{n-1}|{\bf y}_{1:n-1})$   .

In the following, we shall consider a Gaussian kernel ${\cal K}_{h_{n}}={\cal N}(0,h_{n}^{2})$: Eq. (13) then effectively describes a Gaussian mixture and resampling according to this distribution now consists in a two-step procedure, with the particle selection (performed using the systematic resampling algorithm Hol06 ) complemented by a perturbation controlled by the kernel bandwidth $h_{n}$. Each particle generated in the process is thus unique. Note that, if the likelihood function $g_{n}$ is unknown, it can itself be emulated by a kernel density estimator Campillo09 ; ChenPhD .

There is a whole literature on how the kernel bandwidth $h_{n}$ can be chosen. For now, we will follow Ref. Musso12 and pick $h_{n}^{2}=\alpha_{h}\Sigma_{n}$, with

[TABLE]

This choice minimizes the mean integrated squared error $\mathbb{E}\left[||p_{n}^{N}-p_{n}||_{L^{2}}^{2}\right]$, provided the true probability density function $p_{n}$ is Gaussian (see Appendix A for the proof in one dimension with uniform weights). In practice, the true covariance matrix $\Sigma_{n}$ is unknown and is therefore replaced by an estimator $\Sigma_{n}^{N}$ based on the weighted samples $({\bf x}_{n|n-1}^{(i)},w_{n}^{(i)})$:

[TABLE]

with ${\bf\mu}_{n}^{N}=\sum_{i=1}^{N}w_{n}^{(i)}{\bf x}_{n|n-1}^{(i)}$. As an important corollary, because $\lim_{N\to\infty}h_{n}=0$, we retrieve asymptotically that, just as for the empirical probability measure, the regularized estimator converges towards the true posterior density:

[TABLE]

2.2 Problem statement

While the general strategy underlying the RPF seems very principled, we proceed to highlight a problem which can be encountered when resampling is performed on a systematic basis, i.e. at every step, as generally done. Recall that, given a prior probability measure $\pi_{0}(d{\bf x}_{0})$ and a series of observations $\{{\bf y}_{1:n}\}$, we seek to determine the true posterior distribution $\hat{{\bf x}}_{n}\sim\pi_{n}(d{\bf x}_{n}|{\bf y}_{1:n})$. In particular, we wish to estimate how the RPF estimator $\hat{{\bf x}}_{n}^{N}\sim p_{n}^{N}({\bf x}_{n}|{\bf y}_{1:n})$ fares with respect to $\hat{{\bf x}}_{n}$. In what follows, we shall provide evidence that, for a finite number of particles, the RPF estimator features a systematic deviation from the true posterior density which increases with the number of observations. This deviation can be made explicit by comparing the root mean squared error of the RPF with that of the true posterior density:

[TABLE]

where $||.||$ indicates the usual Euclidean norm, and introducing the mean value ${\bf\mu}_{n}:=\mathbb{E}[\hat{{\bf x}}_{n}| {\bf y}_{1:n}]$ and the covariance matrix $\Sigma_{n}:=\mathbb{E}[(\hat{{\bf x}}_{n}-{\bf\mu}_{n})(\hat{{\bf x}}_{n}-{\bf\mu}_{n})^{T}]$. For illustrative purposes, consider the very simple case of a one-dimensional stationary state $x_{n+1}=x_{n}$ with time-independent additive observation noise $\eta\sim{\cal N}(0,R)$: $y_{n}=x_{n}+\eta$. This is of course an elementary textbook setting, which is in fact a special case of the Kalman filter framework which shall be reviewed in the following section. The advantage of this setup is to provide us with the exact result with respect to which particle filter expectations can be benchmarked. In the left frame of Fig. 1, we plot the analytical solution and compare it with the outcome of the RPF when resampling is performed at each step.

As can be clearly seen, while the analytical expectation steadily decreases with each new observation, the RMSE computed using the RPF rather quickly saturates and remains lower bounded by a positive constant no matter how many observations are considered. The explanation for this behavior essentially stems from the fact that resampling from the regularized estimator, while allowing to explore particle space, also acts as a perturbation which increases the variance of the estimator. This perturbation becomes detrimental for the estimation process in the long run, as is manifest in the left frame of Fig. 1. While this problem was already acknowledged in a slightly different context in the early days of particle filters West93 , some of its implications may not have been fully grasped at the time. In particular, we claim that the over-estimation of the variance prevents the estimator from remaining asymptotically efficient as the number of observations increases, in a sense which shall be defined shortly. Mathematical arguments supporting this claim for linear state-space models will be provided in the following section, while numerical simulations indicating that our results extend to non-linear models will be presented in section 4.

3 Theoretical results in the linear Gaussian case

We now address the origin of the saturation effect in the estimated RMSE in more detail and claim that it can be understood using results which we obtain in the setting of linear Gaussian state-space models. Before providing mathematical evidence supporting our claim, we begin by recalling standard results for the derivation of the filtered probability density function in the linear Gaussian case, which is the framework of the famous Kalman filter.

3.1 Kalman filter framework

The framework of the Kalman filter Kalman60 ; Kalman61 is obtained by considering the special case where transition and measurement functions are linear, and assuming that noise and prior probability density functions are Gaussian distributed with respect to the Lebesgue measure. The HMM equations then read

[TABLE]

where $\{A_{n}\}_{n\geq 1}$ and $\{B_{n}\}_{n\geq 1}$ are referred to as transition and measurement matrices, and with $\epsilon_{n}\sim{\cal N}(0,Q_{n})$ and $\eta_{n}\sim{\cal N}(0,R_{n})$. In this situation, given a prior $\hat{{\bf x}}_{0}\sim{\cal N}(\mu_{0},\Sigma_{0})$, the filtering problem can be solved exactly: analytical expressions can be obtained for the posterior $\hat{{\bf x}}_{n}\sim{\cal N}(\mu_{n},\Sigma_{n})$ which is fully characterized by its mean value ${\bf\mu}_{n}$ and (symmetric positive definite) covariance matrix $\Sigma_{n}$. Their expressions are provided by the recursive formulae

[TABLE]

with ${\bf\mu}_{n|n-1}=A_{n}{\bf\mu}_{n-1}$ and $\Sigma_{n|n-1}=A_{n}\Sigma_{n-1}A_{n}^{T}+Q_{n}$ (see Appendix B for a derivation). A simple way of obtaining these results is to make use of the Bayesian conjugacy property of Gaussian distributions: in other words, given a Gaussian prior and a Gaussian likelihood, the posterior is also Gaussian. Applying Bayes rule to the predicted and posterior probability measures, respectively as

[TABLE]

and

[TABLE]

thus allows for a straightforward derivation of Eq. (19). We next focus on the special case where all states are assumed stationary. This corresponds to taking $A_{k}=\mathbb{I}$ and $Q_{k}=0$ for $k\geq 1$ in Eq. (18). If we further fix the measurement covariance matrices to unity $B_{k}=\mathbb{I}$ and assume that observation noise covariance matrices are time-independent, one can show using the matrix inversion lemma that

[TABLE]

with $\langle{\bf y}\rangle_{n}:=\sum_{j=1}^{n}{\bf y}_{j}/n$ the maximum likelihood estimator. This situation corresponds to a multi-dimensional generalization of the stationary model we considered earlier. Note that $\Sigma_{n}$ reaches the Cramer-Rao bound $R/n$ as $n\to\infty$, while $\lim_{n\to\infty}{\bf\mu}_{n}={\bf x}_{0}$: we thus retrieve that, in accordance with the Bernstein-von Mises theorem, the Bayes estimator $\hat{{\bf x}}_{n}$ is asymptotically efficient and coincides in that limit with the maximum likelihood estimator, as the dependence on the choice of prior vanishes. The corresponding $\text{RMSE}=\mathcal{O}(n^{-1/2})$ also vanishes asymptotically. This behavior might seem like a natural byproduct in the limit of a large number of observations, since the latter might be expected to eventually provide sufficient information for an arbitrarily accurate estimation of the true state. There are, however, elementary counterexamples to this fact which can be very easily proven using the Kalman filter equations we have derived. For example, if we restore a simple time-independent transition matrix $A=a\mathbb{I}$ to the above estimation paradigm, one can easily show that the asymptotic behavior of $\Sigma_{n}$ depends on the magnitude of $a$: while if $a<1$, the covariance matrix decays exponentially fast, if $a>1$ however, then the covariance matrix remains lower bounded as $\Sigma_{n}\geq(1-a^{-2})R$. Note that if $a=1$, we recover the algebraic dependence as stated in Eq. (22). Another enlightening case is when state noise $Q\neq 0$ is included. Given time-independent state and observation noises, one can show that $\Sigma_{n}$ necessarily remains lower bounded by a positive matrix HMM .

3.2 Regularized RMSE asymptotics

Let us now investigate how the above results can be transposed in the framework of the RPF. For a Gaussian regularizing kernel, Bayesian conjugacy ensures that the posterior distribution Eq. (13) is asymptotically Gaussian:

[TABLE]

with mean

[TABLE]

and covariance matrix

[TABLE]

with ${\bf\mu}_{n|n-1}^{N}=A_{n}{\bf\mu}_{n-1}^{N}$ and $\Sigma_{n|n-1}^{N}=A_{n}\Sigma_{n-1}^{N}A_{n}^{T}+Q_{n}$, and where $\alpha_{h}$ is defined in Eq. (14). Note that, as $N\to\infty$, we recover the standard Kalman result from Eq. (19) as we should, since the kernel acts as a Dirac delta-function in Eq. (23) as $\alpha_{h}\to 0$. Additionnally, while we work in the limit of $N\to\infty$ and choose to conserve the $N$-dependence of the kernel bandwidth, we explicitly assume the Monte Carlo error $\propto N^{-1/2}$ is negligible when replacing the sum by the integral in Eq. (23). We now proceed to illustrate the consequences of Eqs. (24,25) by returning for simplicity to the case of stationary states, where exact expressions for the posterior moments can be obtained by solving Eqs. (24,25) recursively. For reasons which will be clear later, we shall allow the possibility for $\alpha_{h}$ to be time-dependent. Note that our analysis should carry over to any dynamical system described by Eqs. (18) for which the Bayes estimator is asymptotically efficient: this includes systems with noiseless and compact supported state dynamics, and for which the observation model is regular in some sense CarmierKalman .

Lemma 1

Consider the recursively defined set of equations for $n\geq 1$ on the couple of variables $(\mu_{n},\Sigma_{n})$,

[TABLE]

where $\{\alpha_{n}\}_{n\geq 1}$ is a sequence of positive real numbers, and $R$ and $\{{\bf y}\}_{n\geq 1}$ are defined as in Eq. (22). Starting from the initial condition $(\mu_{0},\Sigma_{0})$, the solution to the above set of equations can be expressed as

[TABLE]

with $\alpha_{0}=0$.

Proof

This can be straightforwardly proven by recurrence. ∎

As a first application of lemma 1, consider the case where the sequence $\{\alpha_{n}\}_{n\geq 1}$ is time-independent. This is our original motivation, corresponding to the choice of bandwidth indicated in Eq. (14): $\alpha_{n}=\alpha_{h}\;\forall n\geq 1$. Expressions for the posterior moments simplify as

[TABLE]

We obtain as a direct consequence that, if $\alpha_{h}>0$ (for finite $N$), the covariance matrix does not converge to zero as the number of observations increases.

Theorem 3.1 (RPF fixed point)

The RPF estimator for the posterior covariance matrix converges exponentially fast towards a non-vanishing fixed point:

[TABLE]

Likewise, the estimated mean value converges exponentially fast towards a random variable

[TABLE]

centered around the observation mean value ${\bf x}_{n}={\bf x}_{0}$, with residual covariance

[TABLE]

As a result, the asymptotic $\text{RMSE}=\mathcal{O}(\sqrt{\alpha_{h}})$ remains lower bounded by a positive constant as soon as the number of particles is finite, no matter how many observations.

As we have already mentioned, this situation can also occur in the original Kalman filter (retrieved in the limit of an infinite number of particles) as soon as process noise $Q\neq 0$ is included. In the present case, the reason for this phenomenon is that the kernel bandwidth used for resampling acts as a perturbation which increases the variance and competes with the variance reduction due to each new observation. The problem is that the variance reduction rate with $\Sigma_{n}^{N}$ is faster (quadratic) than that of the bandwidth (linear): $\Sigma_{n+1}^{N}-\Sigma_{n}^{N}=(\alpha_{h}R-\Sigma_{n}^{N})(R+\Sigma_{n}^{N})^{-1}\Sigma_{n}^{N}$. Arrives therefore a moment when the ever smaller piece of information provided by the new observation is exactly compensated by the bandwidth perturbation.

A complementary and enlightening interpretation of our results can be put forward by noting that each observation in Eq. (30) contributes differently to the estimated posterior mean value, i.e. each ${\bf y}_{j}$ carries a different weight. In particular, the weighting ratio between first (oldest) and last (most recent) observations roughly scales like $e^{-n\alpha_{h}}$, which signifies that the RPF estimator loses memory of previous observations exponentially fast. This loss of memory is directly responsible for the inability of the estimator to process the information provided by the full history of observations as exhaustively as the Kalman estimator, since it acts as an effective temporal resolution. The lower bound on the RMSE obtained in theorem 3.1 is thus simply the translation of this temporal resolution in state space.

Returning to Fig. 1, we check the validity of our results by comparing them with numerical simulations of the one-dimensional stationary HMM introduced in section 2.2. Interestingly, although conclusions of theorem 3.1 based on Eqs. (24,25) were derived in the asymptotic limit $N\to\infty$, numerical simulations support their well-foundedness already for $N=1000$ particles. The late time fluctuations in the RMSE of the RPF correspond to those of mean value in Eq. (31). Next, because $\lim_{N\to\infty}\alpha_{h}=0$, we expect according to theorem 3.1 the lower bound to decrease as $N$ increases. This is clearly illustrated by the blue curve in the left frame of Fig. 1. Increasing the observation noise variance obviously leads to the opposite effect (green curve). We plotted as a dashed black line the analytical prediction based on the Kalman filter from Eq. (28), where the observation noise variance $R/\Sigma_{0}$ was quenched from the value 0.25 to 1 at $n=300$. As one can see, the computed posterior variance on the state reacts to the quench by increasing and converging towards the new asymptotic fixed point. Notice that increasing $d_{x}$ is also detrimental (not shown), since $\lim_{d_{x}\to\infty}\alpha_{h}=1$ according to Eq. (14). Next, we show in the right frame that other bandwidth selection criteria such as Silverman and Scott’s rule of thumbs or the direct plug-in method by Sheather and Jones Sheather04 lead to the same problem when resampling is performed at each step. Another natural strategy is to perform resampling only occasionally. However, resampling periodically every $p$ step only allows reducing the location of the fixed point roughly a factor $1/p$, as is shown in the following corollary.

Corollary 1

Given a resampling period $p\in\mathbb{N}^{*}$ and an integer $0\leq q<p$, the RPF estimator for the posterior covariance matrix $\Sigma_{np+q}^{N}$ converges exponentially fast towards the non-vanishing q-dependent solution:

[TABLE]

As a result, the asymptotic covariance matrix fluctuates between the upper bound $\alpha_{h}R/p$ and the lower bound $\alpha_{h}R/(p+(p-1)\alpha_{h})$. Likewise, the estimated mean value converges exponentially fast towards a random variable centered around the observation mean value ${\bf x}_{n}$, with a q-dependent residual covariance

[TABLE]

This translates in an asymptotic $\text{RMSE}=\mathcal{O}(\sqrt{\alpha_{h}/p})$.

Proof

This can be proven by deriving from Eq. (27) the new expression for the posterior covariance matrix when $\alpha_{np+q}=\alpha_{h}\delta_{q,0}$:

[TABLE]

where we defined $[j-1]_{p}$ as the integer part of $(j-1)/p$ . Similar calculations allow obtaining the corresponding expression for the mean value:

[TABLE]

∎

That the RMSE remains lower bounded by a positive constant when resampling is performed periodically is illustrated numerically in the left frame of Fig. 2, where the analytical lower bound (32) is indicated by the dotted line.

3.3 Kernel bandwidth modulation

One is thus left with two possibilities in order to allow the estimator to keep up with the Kalman filter prediction: either resample less and less often with time (which seems compatible with the idea of resampling only if the ESS falls lower than a critical threshold, see subsection 3.4), or choose the bandwidth parameter differently so that consistency of the estimator is maintained despite performing resampling at each time. To begin with, it is clear from theorem 3.1 that choosing $\{\alpha_{n}\}_{n\geq 1}$ such that $\lim_{n\to\infty}\alpha_{n}=0$ enforces that $\lim_{n\to\infty}\text{RMSE}(\hat{{\bf x}}_{n},{\bf x}_{n})=0$ (this can easily be proven by contradiction). This idea was actually hinted at in ChenPhD . We proceed one step further by characterizing in theorem 3.2 which sequences ensure that the rate of convergence of the RPF estimator is optimal in the sense defined below.

Definition 1 (Optimal convergence rate)

The rate of convergence of an estimator $\hat{{\bf x}}_{n}$ towards the true state ${\bf x}_{n}$ as the number of observations $n\to\infty$ will be called optimal if it is comparable with that given by the inverse of the Fisher information matrix: $\text{RMSE}(\hat{{\bf x}}_{n},{\bf x}_{n})=\mathcal{O}(n^{-1/2})$.

Since the Kalman filter estimator is asymptotically efficient, this implies in particular that the optimal rate of convergence is that given by the Kalman filter.

Theorem 3.2

Let $\{\alpha_{n}\}_{n\geq 1}$ be a sequence of positive numbers as in Eq. (27). If $\alpha_{n}=\mathcal{O}(n^{-1})$, the RPF estimator converges at an optimal rate. In contrast, given any $\epsilon\in]0,1/2[$, if $\alpha_{n}\approx n^{\epsilon-1}$, then $\exists\;c\in\mathbb{R}^{+}\;\forall n\geq 1\;\text{RMSE}\geq cn^{(\epsilon-1)/2}$.

Proof

Suppose that $\exists\;c\in\mathbb{R}^{+}\;\forall n\geq 1\;n\alpha_{n}\leq c$. Then, using that $x-x^{2}/2<\log{(1+x)}<x$ for $x>0$, one has that $\prod_{j=1}^{n}(1+\alpha_{j})=\mathcal{O}(n^{c})$, and thus also that $\sum_{j=1}^{n}\prod_{j=0}^{k-1}(1+\alpha_{j})=\mathcal{O}(n^{c+1})$. Plugging this in Eq. (27) implies that $\Sigma_{n}^{N}=\mathcal{O}(R/n)$ and therefore that the RPF estimator converges at an optimal rate. In contrast, if $\alpha_{n}\approx n^{\epsilon-1}$, then the same reasoning implies that $\prod_{j=1}^{n}(1+\alpha_{j})=\mathcal{O}(\exp{(n^{\epsilon}/\epsilon)})$. Making use of the following integral inequality

[TABLE]

we deduce that $\exists\;c\in\mathbb{R}^{+}\;\forall n\geq 1\;\Sigma_{n}^{N}\geq cn^{\epsilon-1}R$. Finally, noting that $\sum_{j=0}^{n-1}\exp{(2j^{\epsilon}/\epsilon)}\geq\exp{(2(n^{\epsilon}-1)/\epsilon)}$, we conclude that $\exists\;c\in\mathbb{R}^{+}\;\forall n\geq 1\;\text{RMSE}\geq cn^{(\epsilon-1)/2}$. ∎

This result still leaves us with a wide range of sequences at our disposal. Finding a suitable statistical criterion both justifying this type of kernel bandwidth modulation (and therefore different from the mean integrated squared error used to derive Eq. (14)) and allowing to identify a sequence which is optimal in some sense remains an open question. Let us look at two examples for the sequence $\{\alpha_{n}\}_{n\geq 1}$.

• $\bullet$

Example 1: $\alpha_{n}=(n+\alpha_{h}^{-1})^{-1}$.

For this choice, products and sums in Eq. (27) can be computed exactly and we get for the covariance matrix $\Sigma_{n}^{N}\underset{n\to\infty}{\approx}2R/n$, and for the mean value

[TABLE]

centered around the observation mean value ${\bf x}_{n}$ and with residual covariance $\mathbb{E}[({\bf\mu}_{n}^{N}-{\bf x}_{n})({\bf\mu}_{n}^{N}-{\bf x}_{n})^{T}]\underset{n\to\infty}{\approx}4R/(3n)$. Thus, optimal convergence rate is recovered, as expected from theorem 3.2. Interestingly, as compared to Eq. (30), the dependence of Eq. (37) on the observations is long-ranged: the weighting ratio between the last and the first observations scales linearly with $n$ instead of exponentially as in Eq. (30). This illustrates that consistency of the regularized estimator requires that memory of past observations should be kept long enough. Note that such a requirement, while understandable for parameter estimation, may be irrelevant for tracking applications where memory of the past might be considered relatively useless.

• $\bullet$

Example 2: $\alpha_{n}=\alpha_{h}e^{-n\alpha_{h}}$.

In this case, one cannot compute sums and products in Eq. (27) exactly but one can easily show that $\prod_{j=1}^{n}(1+\alpha_{h}e^{-j\alpha_{h}})$ converges to a non-zero limit using that it is a strictly increasing bounded sequence:

[TABLE]

This result allows immediately deducing that $\sum_{j=1}^{n}\prod_{k=0}^{j-1}(1+\alpha_{h}e^{-k\alpha_{h}})=\mathcal{O}(n)$, and thus that $\Sigma_{n}^{N}=\mathcal{O}(R/n)$ and that the mean value

[TABLE]

is centered around the observation mean value ${\bf x}_{n}$ with residual covariance $\mathbb{E}[({\bf\mu}_{n}^{N}-{\bf x}_{n})({\bf\mu}_{n}^{N}-{\bf x}_{n})^{T}]=\mathcal{O}(R/n)$, using that $\prod_{k=0}^{j-1}(1+\alpha_{h}e^{-k\alpha_{h}})^{2}$ is also bounded. Optimal convergence rate is thus once more recovered, this time with a weighting ratio between different observations bounded by $\prod_{j=1}^{\infty}(1+\alpha_{h}e^{-j\alpha_{h}})$.

This strategy of modulating the bandwidth by enforcing it to decrease fast enough is positively supported by numerical simulations in the right frame of Fig. 1. It thus constitutes a functional method to overcome the limitations in the RPF we discussed earlier. Note that a different strategy discussed in Ref. West93 consists in counterbalancing the variance originating from the bandwidth perturbation by rescaling the particle values towards their center of mass. Centering the kernel densities around the shrunken values $a{\bf x}_{n|n-1}^{(i)}+(1-a){\bf\mu}_{n}^{N}$, and taking $a=\sqrt{1-\alpha_{h}}$, one can indeed verify that the posterior density in the RPF asymptotically coincides with that expected from the Kalman filter approach.

3.4 Strategies based on the resampling ratio

We now examine another possibility, which is to perform resampling only if the ESS introduced in Eq. (10) becomes lower than a determined threshold which we shall refer to as the resampling ratio $\text{ESS}_{\text{crit}}$ DelMoral12 . In the left frame of Fig. 2, we compare the posterior RMSE for different values of this threshold, both for the RPF and the SIR algorithm. As one can clearly see, as soon as the resampling ratio is smaller than 1, the RPF yields results in agreement with the exact solution. The SIR algorithm provides comparable results, whatever the value of the resampling ratio. However, if the location $\mu_{0}$ of the prior density is taken farther away from the correct value, the situation changes since estimating the posterior mean value $\mu_{n}$ accurately will increasingly depend on the ability of particles to explore phase space well enough to diffuse towards the correct value. This is illustrated in the right frame of Fig. 2, where $\mu_{0}$ was taken to be $3\sqrt{\Sigma_{0}}$ away from the observed mean value, instead of $\sqrt{\Sigma_{0}}$ until now. In this case, the SIR algorithm almost always fails to correctly estimate the posterior variance due to sample impoverishment arising as a direct consequence of severe weight degeneracy. In fact, the estimated variance often goes to zero, reflecting the fact that only a single particle eventually survives the resampling process. On the other hand, the RPF proves very robust in this situation and provides correct results as before, despite the fact that its RMSE has a tendency to underestimate the optimal Kalman expectation. Note that the curves displayed on Fig. 2 are outputs of single particle filter simulations (no averaging is performed).

Resampling conditional on a sub-unit value of $\text{ESS}_{\text{crit}}$ thus also appears to resolve the RPF shortcomings discussed earlier. One can actually show that the resampling frequency for a sub-unit resampling ratio decreases exponentially fast with the number of observations, thereby justifying that the posterior RMSE should no longer be lower bounded by a positive constant. We provide below a formal proof of this statement in the case of the stationary HMM.

Lemma 2

Assume that $\text{ESS}_{\text{crit}}<1$. Then the spacing between consecutive resampling times increases exponentially fast with the number of observations for a stationary HMM.

Proof

Our starting point is the observation that

[TABLE]

For simplicity, we consider the case of the one-dimensional stationary state, where integration with respect to the Markov kernels reduces to $\prod_{j=1}^{n}g_{j}(y_{j}|x_{0})$. It is then a simple task to compute the effective sample size:

[TABLE]

We see that $\text{ESS}_{n}$ depends sensitively on how close the ML $\langle y\rangle_{n}$ is from the typical particle value $\mu_{0}$ and on the quantity $\gamma_{n}:=R/(n\Sigma_{0})$. The number of observations essentially reduces the effective observation noise variance. If resampling is performed at step $n$, then one can express the ESS at step $n+m$ asymptotically as

[TABLE]

The value of $m$ required for the ESS to fall below a critical threshold $\text{ESS}_{\text{crit}}N$, with $\text{ESS}_{\text{crit}}\in[0,1]$, can then easily be calculated and we get

[TABLE]

with

[TABLE]

Iterating, we thus see that the number of steps $m_{k}$ asymptotically required to perform resampling for the $k^{\text{th}}$ time grows exponentially fast:

[TABLE]

∎

Note that similar conclusions were reached in Chopin02 . We provide numerical simulations supporting this claim in Fig. 3.

Interestingly, we compared the results we obtain for noisy observations (inset) with that for noiseless observations (symbols in main frame). The latter can be interpreted as representing the most favourable scenario, in which observations are always equal to the noiseless signal, and shall therefore be referred to from now on as the oracle limit, by analogy with the influential decision-theoretic concept Candes06 . As one can notice in the inset, the exponential growth of the resampling step sizes can be occasionally violated if observations cause large fluctuations in the exponential term of Eq. (41). These fluctuations also explain why the ESS is not necessarily monotonous between consecutive resampling times for noisy observations.

Equipped with lemma 2, we can now conclude our demonstration that the RPF estimator is well-behaved as soon as $\text{ESS}_{\text{crit}}<1$.

Theorem 3.3

Assume that $\text{ESS}_{\text{crit}}<1$. Then the rate of convergence of the RPF estimator is optimal in the sense of definition 1.

Proof

According to the previous lemma, we know that the spacing between consecutive resampling times increases exponentially fast, with rate $1+\beta_{\text{crit}}$. As a consequence,

[TABLE]

and thus $\exists\;c\in\mathbb{R}^{+}\;\forall n\geq 1\;\prod_{j=1}^{n}(1+\alpha_{j})=\mathcal{O}(n^{c})$. The same reasoning as that followed in the proof of theorem 3.2 allows us to conclude. ∎

4 Extension to non-linear models

We now provide numerical evidence that the above results qualitatively hold in the case of non-linear models. We shall consider two examples: the logistic map, which is notably famous for possessing regions of parameter space displaying chaos, and a multi-dimensional plant growth model.

4.1 Logistic map

The logistic map, popularized by Robert May May76 , is a discretized version of the logistic equation first introduced by Verhulst in the 19${}^{\text{th}}$ century to describe population dynamics with a finite carrying capacity. It is described by the deceptively simple equation $x_{n+1}=ax_{n}(1-x_{n})$, with ${\cal X}=[0,1]$ for $a\in[0,4]$. The nature of the state dynamics can change drastically depending on the value of $a$: it converges to a stationary solution for $0\leq a\leq 3$, oscillates between multiple values for $3<a\leq a_{c}$, and displays chaotic behavior for $a>a_{c}$, where $a_{c}\simeq 3.57$, with islands of stability found amid. The effect of the initial condition, on the other hand, is often of litte relevance and shall be assumed known to us: $x_{0}=0.5$. In the following, we seek to estimate $a$ given multiplicative noisy observations $y_{n}=x_{n}\eta$, with $\eta\sim\log{\cal N}(1,R)$. Starting from a prior distribution $a^{(i)}\sim{\cal N}(3.0,0.3)$, the resulting inference of parameter $a$ given observations generated with the value $a^{*}=3.33$ and $R=0.1$ is displayed in Fig. 4. Saturation of $\text{RMSE}(\hat{a}_{n},a^{*})$ when resampling is performed at each step is clearly apparent (solid red). This effect disappears using either of the resampling strategies we discussed earlier, however, as the comparison with a benchmark simulation for $N=10^{5}$ particles (solid black) testifies. Interestingly, we compared the results we obtain for noisy observations (thick lines) with that in the oracle limit (dotted lines). The RMSE obtained in the latter case acts as a lower bound with respect to that obtained for noisy observations, with numerical simulations even indicating that $\lim_{n\to\infty}\text{RMSE}(\hat{a}_{n},a^{*})=0$ in this case for $\text{ESS}_{\text{crit}}<1$.

If instead we choose to work in the chaotic parameter region and take $a^{*}>a_{c}$, our numerical results (not shown) seemingly point to the fact that correctly inferring the value of $a^{*}$ is always challenging, even in the oracle limit. The same can be told for the estimation of the initial condition $x_{0}$, for a fixed known value of $a$. The reason for this behavior seems to stem from the fact that, in such situations, the information provided by observations on the quantity of interest is a (rapidly) decreasing function of time. This loss of memory on the initial condition ultimately limits the resolution with which inference can be carried, as concentration results recently obtained in a similar context may testify Paulin16 .

4.2 LNAS model

Next, we look at a more elaborate model, which describes the growth of leaf mass $l_{n}$ and root mass $r_{n}$ of a sugarbeet plant as a function of daily temperature $T_{n}$ and radiation $\phi_{n}$. The latter constitute environmental variables, while the former are the hidden states. This model is a simplified version of the LNAS model Cournede13 . In particular, we consider deterministic transition functions. Leaf and root masses are updated as

[TABLE]

where $q_{n}$ is the daily produced biomass and $a_{n}$ is an allocation function. Biomass production is governed by the Beer-Lambert law

[TABLE]

where RUE is the radiation use efficiency and $\rho^{-1}$ the leaf mass density. Allocation depends on the thermal time $\tau_{n}$ as

[TABLE]

where $\gamma$ controls early biomass partition, erf is the error function and $(\mu_{a},\sigma_{a})$ are parameters governing the shape of the allocation function. Thermal time increases monotonically as a function of temperature as $\tau_{n}=\tau_{n-1}+\text{max}[0,T_{n}]$. The system is initialized as $\tau_{0}=r_{0}=0$, and $l_{0}\neq 0$ such that $\rho^{-1}l_{0}\ll 1$. For environmental variables, we use real data obtained from a weather station in France. The filtering problem in the present case then consists in estimating the parameters (RUE, $\gamma$, $\mu_{a}$), given independent prior densities $\text{RUE}^{(i)}\sim{\cal N}(3.8,0.3)$, $\gamma^{(i)}\sim{\cal N}(0.7,0.05)$ and $\mu_{a}^{(i)}\sim{\cal N}(500,50)$, and given multiplicative noisy observations of the hidden states $l_{n}\eta$ and $r_{n}\eta$ with $\eta\sim\log{\cal N}(1,R)$. We shall assume $\rho$ and $\sigma_{a}$ known beforehand.

Because we no longer have access to the optimal posterior density, as opposed to the case of the stationary model, we plot the evolution of posterior standard deviations and mean values of each parameter (instead of the RMSE) and use results obtained by the RPF for $N=10^{6}$ particles and ESS${}_{\text{crit}}=0.5$ as a benchmark. We checked that this corresponds to the optimal posterior density by plotting the results obtained for $N=10^{6}$ particles using the SIR algorithm by comparison, and indeed observed that convergence as a function of $N$ is achieved (see Fig. 5). These simulations were carried by running our algorithms on the computing mesocenter Fusion. Given observations generated with values $(\text{RUE}^{*},\gamma^{*},\mu_{a}^{*})=(3.56,0.625,550)$ and $R=0.1$, the resulting inference for various numbers of particles and resampling strategies is displayed in Fig. 6. We see that the performance of the RPF for $N=1000$ particles and ESS${}_{\text{crit}}=0.5$ is already close to that obtained in the benchmark scenario. Using our kernel bandwidth modulation strategy, we additionally demonstrate that equivalent results can be obtained even when ESS${}_{\text{crit}}=1$, whereas the traditional choice of Eq. (14) leads to an unstable exploration of parameter space (not shown). This holds independently of the chosen parameter. In contrast, expectations based on the bootstrap particle filter or in the absence of any resampling are evidently unable to capture the correct behavior for the posterior variance. As a consequence, these particle filters will fail to provide appropriate credibility regions for the parameters and thereby make any reliable prediction impossible. This is also illustrated in table 1 where we compared the performances of our various resampling strategies, each time averaged over 10 independent runs. These results thus clearly demonstrate the ability of the RPF to outperform such traditional particle filters for a moderate number of particles.

5 Conclusions

To conclude, we demonstrated that the choice of resampling strategy in the regularized particle filter can be paramount, as far as estimation of the posterior probability density is concerned. In particular, using analytical arguments inspired from linear Gaussian state-space models, we proved that resampling on a systematic basis causes the RMSE to be lower bounded by a constant which depends on the number of particles and on the dimensionality of the state space. Numerical evidence supporting our claims, including for non-linear state-space models, was also highlighted. Finally, we showed that if one introduces a resampling ratio, or if the kernel bandwidth is appropriately modulated, then the regularized particle filter clearly outperforms traditional bootstrap particle filters.

As a closing remark, we note that the idea we highlighted in this paper of progressively reducing the bandwidth perturbation seems related to recent approaches known as tempering Svensson17 , which consist in building a probability measure which depends on a tunable parameter allowing one to smoothly transport particles from the prior $\pi_{0}$ to the posterior $\pi_{n}$. Similar ideas have also been invoked to perform Bayesian update progressively Daum08 , and investigating the ramifications of such connections in more detail would certainly be worthy of interest. We also leave for future research the question of which statistical criterion (Kullback-Leibler divergence, Wasserstein metric, etc.) could justify this type of kernel bandwidth modulation.

Acknowledgements.

The authors wish to thank Gautier Viaud for useful conversations at various stages of this work.

Appendix A Optimal bandwidth selection

There are various ways Sheather04 of choosing the bandwidth parameter $h_{n}$ in Eq. (13). The most popular choice relies on minimizing the mean integrated squared error (MISE) between the true posterior probability density function $p$ and the kernel density $p_{h}^{N}(x)=\sum_{i=1}^{N}{\cal K}_{h}(x-x^{(i)})/N$, with $x^{(i)}\sim p(x)$:

[TABLE]

where we introduced the $L^{2}(\mathbb{R})$ norm $p\mapsto||p||_{L^{2}}=\sqrt{\int dz\;p(z)^{2}}$. As explicited above, the MISE can allegedly be decomposed into the integrated squared bias and the integrated variance. Observing that $\mathbb{E}[p_{h}^{N}(x)]=\int dz\;p(z){\cal K}_{h}(x-z)$ and developing $p$ around $x$ as a Taylor series, one obtains:

[TABLE]

Likewise, the variance can be computed by decomposing $\mathbb{E}[p_{h}^{N}(x)^{2}]$ into diagonal and non-diagonal contributions:

[TABLE]

such that $\mathbb{V}[p_{h}^{N}(x)]=(||{\cal K}_{1}||_{L^{2}}^{2}p(x)+o(1))/(Nh)$. Injecting these results in Eq. (50) and minimizing with respect to $h$ is straightforward and provides the optimal choice:

[TABLE]

Note that this implies that $\text{MISE}\propto N^{-4/5}$. The only problem left is that this choice depends on the second moment of the underlying true probability density function which is generally unknown. There are different strategies to estimate this quantity, which we will not delve into Sheather04 . If we assume that both $p$ and ${\cal K}$ are Gaussian distributed, the expression for $h$ reduces to the simple $h_{\text{MISE}}^{5}=4\Sigma_{p}^{5/2}/(3N)$ which coincides with Eq. (14) for $d_{x}=1$.

Appendix B Kalman filter equations

The Kalman filter equations (19) can be derived recursively in two steps, starting from the probability density function at time $n-1$. The first step is the prediction step, in which we seek to determine $p(\hat{{\bf x}}_{n}| {\bf y}_{1:n-1})$, in other words how the probability density function is propagated between two time steps. Using Eq. (18), one can easily show that

[TABLE]

using that state noise has null expectancy, and

[TABLE]

using that state noise is uncorrelated with the dynamics. The second step is the update step, in which we seek to compute the posterior moments ${\bf\mu}_{n}$ and $\Sigma_{n}$, expressions for which are traditionally obtained by enforcing unbiasedness of the estimator and minimizing the mean squared error, as we will now show. We begin by writing the posterior mean as a linear combination between the prior mean and the observation, that is ${\bf\mu}_{n}=K^{\prime}_{n}{\bf\mu}_{n-1}+K_{n}{\bf y}_{n}$. Enforcing ${\mathbb{E}}[\hat{{\bf x}}_{n}-{\bf\mu}_{n}| {\bf y}_{1:n-1}]=0$ leads to $K^{\prime}_{n}=(\mathbb{I}-K_{n}B_{n})A_{n}$, and thus:

[TABLE]

which expresses the posterior mean as a linear combination between the predicted mean and the observation, with a tradeoff controlled by the so-called Kalman gain matrix $K_{n}$, an expression for which will be derived below. We begin by computing $\Sigma_{n}$, which follows directly from the previous equation:

[TABLE]

Here we defined $\text{Cov}[\hat{{\bf x}}]=\mathbb{E}[\hat{{\bf x}}\hat{{\bf x}}^{T}]$ as a shorthand and used the fact that observation noise is uncorrelated with the dynamics. The Kalman gain is then obtained by minimizing the mean squared error, which corresponds to the trace of the posterior covariance matrix:

[TABLE]

where derivation with respect to a matrix is performed element-wise. Solving Eq. (58) yields the Kalman gain expression

[TABLE]

where matrix inversion is justified by positive definiteness. Injecting this in Eq. (57), we obtain a much simpler expression for the posterior covariance matrix:

[TABLE]

From Eqs. (59,60), we also derive the useful and synthetic relation $\Sigma_{n}B_{n}^{T}=K_{n}R_{n}$. Finally, using this and rearranging Eqs. (56,60), we obtain the central result Eq. (19). Notice that we made no real use of the assumption of Gaussian distributions. In fact, the above results hold true irrespective of the underlying probability density function, provided noises are centered. The distinctive feature of the Gaussian case is that the knowledge of the first two moments allows recovering the entire distribution. Specifically, for Gaussian probability density functions, the Kalman filter equations can actually also be derived without any particular assumptions using Bayes rule.