Nudging the particle filter
Ömer Deniz Akyıldız Joaquín Míguez











Ö. D. Akyildiz 
              University of Warwick

              The Alan Turing Institute

              omer.akyildiz@warwick.ac.uk           

           J. Míguez Universidad Carlos III de Madrid & Instituto de Investigación Sanitaria Gregorio Marañón







Nudging the particle filterThis work was partially supported by Ministerio de Economía y Competitividad of Spain (TEC2015-69868-C2-1-R ADVENTURE), the Office of Naval Research Global (N62909-15-1-2011), and the regional government of Madrid (program CASICAM-CM S2013/ICE-2845).

    Ömer Deniz Akyildiz         Joaquín Míguez 

    Received: date / Accepted: date
======================================================================================================================================================================================================================================================================================




We investigate a new sampling scheme aimed at improving the performance of particle filters whenever (a) there is a significant mismatch between the assumed model dynamics and the actual system, or (b) the posterior probability tends to concentrate in relatively small regions of the state space. The proposed scheme pushes some particles towards specific regions where the likelihood is expected to be high, an operation known as nudging in the geophysics literature. We re-interpret nudging in a form applicable to any particle filtering scheme, as it does not involve any changes in the rest of the algorithm. Since the particles are modified, but the importance weights do not account for this modification, the use of nudging leads to additional bias in the resulting estimators. However, we prove analytically that nudged particle filters can still attain asymptotic convergence with the same error rates as conventional particle methods. Simple analysis also yields an alternative interpretation of the nudging operation that explains its robustness to model errors. Finally, we show numerical results that illustrate the improvements that can be attained using the proposed scheme. In particular, we present nonlinear tracking examples with synthetic data and a model inference example using real-world financial data.  







§ INTRODUCTION







 §.§ Background

State-space models (SSMs) are ubiquitous in many fields of science and engineering, including weather forecasting, mathematical finance, target tracking, machine learning, population dynamics, etc., where inferring the states of dynamical systems from data plays a key role.

A SSM comprises a pair of stochastic processes (x_t)_t≥ 0 and (y_t)_t≥ 1 called signal process and observation process, respectively. The conditional relations between these processes are defined with a transition and an observation model (also called likelihood model) where observations are conditionally independent given the signal process, and the latter is itself a Markov process. Given an observation sequence, y_1:t, the filtering problem in SSMs consists in the estimation of expectations with respect to the posterior probability distribution of the hidden states, conditional on y_1:t, which is also referred to as the filtering distribution.

Apart from a few special cases, neither the filtering distribution nor the integrals (or expectations) with respect to it can be computed exactly; hence, one needs to resort to numerical approximations of these quantities. Particle filters (PFs) have been a classical choice for this task since their introduction by <cit.>; see also <cit.>. The PF constructs an empirical approximation of the posterior probability distribution via a set of Monte Carlo samples (usually termed particles) which are modified or killed sequentially as more data are taken into account. These samples are then used to estimate the relevant expectations. The original form of the PF, often referred to as the bootstrap particle filter (BPF), has received significant attention due to its efficiency in a variety of problems, its intuitive appeal and its straightforward implementation. A large body of theoretical work concerning the BPF has also been compiled. For example, it has been proved that the expectations with respect to the empirical measures constructed by the BPF converge to the expectations with respect to the true posterior distributions when the number of particles is large enough <cit.> or that they converge uniformly over time under additional assumptions related to the stability of the true distributions <cit.>.


Despite the success of PFs in relatively low dimensional settings, their use has been regarded impractical in models where (x_t)_t≥ 0 and (y_t)_t≥ 1 are sequences of high-dimensional random variables. 
In such scenarios, standard PFs have been shown to collapse <cit.>. This problem has received significant attention from the data assimilation community. The high-dimensional models which are common in meteorology and other fields of Geophysics are often dealt with via an operation called nudging <cit.>. Within the particle filtering context, nudging can be defined as a transformation of the particles, which are pushed towards the observations using some observation-dependent map <cit.>. If the dimensions of the observations and the hidden states are different, which is often the case, a gain matrix is computed in order to perform the nudging operation. In <cit.>, nudging is performed after the sampling step of the particle filter. The importance weights are then computed accordingly, so that they remain proper. Hence, nudging in this version amounts to a sophisticated choice of the importance function that generates the particles. It has been shown (numerically) that the schemes proposed by <cit.> can track high-dimensional systems with a low number of particles. However, generating samples from the nudged proposal requires costly computations for each particle and the evaluation of weights becomes heavier as well. It is also unclear how to apply existing nudging schemes when non-Gaussianity and nontrivial nonlinearities are present in the observation model.

A related class of algorithms includes the so-called implicit particle filters (IPFs) <cit.>. Similar to nudging schemes, IPFs rely on the principle of pushing particles to high-probability regions in order to prevent the collapse of the filter in high-dimensional state spaces. In a typical IPF, the region where particles should be generated is determined by solving an algebraic equation. This equation is model dependent, yet it can be solved for a variety of different cases (general procedures for finding solutions are given by <cit.> and <cit.>). The fundamental principle underlying IPFs, moving the particles towards high-probability regions, is similar to nudging. Note, however, that unlike IPFs, nudging-based methods are not designed to guarantee that the resulting particles land on high-probability regions; it can be the case that nudged particles are moved to relatively low probability regions (at least occasionally). Since an IPF requires the solution of a model-dependent algebraic equation for every particle, it can be computationally costly, similar to the nudging methods by <cit.>.  Moreover, it is not straightforward to derive the map for the translation of particles in general models, hence the applicability of IPFs depends heavily on the specific model at hand.









 §.§ Contribution


In this work, we propose a modification of the PF, termed the nudged particle filter (NuPF) and assess its performance in high dimensional settings and with misspecified models. Although we use the same idea for nudging that is presented in the literature, our algorithm has subtle but crucial differences, as summarized below.



  * First, we define the nudging step not just as a relaxation step towards observations but as a step that strictly increases the likelihood of a subset of particles. This definition paves the way for different nudging schemes, such as using the gradients of likelihoods or employing random search schemes to move around the state-space. In particular, classical nudging (relaxation) operations arise as a special case of nudging using gradients when the likelihood is assumed to be Gaussian. Compared to IPFs, the nudging operation we propose is easier to implement as we only demand the likelihood to increase (rather than the posterior density). Indeed, nudging operators can be implemented in relatively straightforward forms, without the need to solve model-dependent equations.

  * Second, unlike the other nudging based PFs, we do not correct the bias induced by the nudging operation during the weighting step. Instead, we compute the weights in the same way they would be computed in a conventional (non-nudged) PF and the nudging step is devised to preserve the convergence rate of the PF, under mild standard assumptions, despite the bias. 


Moreover, computing biased weights  is usually faster than computing proper (unbiased) weights. Depending on the choice of nudging scheme, the proposed algorithm can have an almost negligible computational overhead compared to the conventional PF from which it is derived.


  * Finally, we show that a nudged PF for a given SSM (say _0) is equivalent to a standard BPF running on a modified dynamical model (denoted _1). In particular, model _1 is endowed with the same likelihood function as _0 but the transition kernel is observation-driven in order to match the nudging operation. As a consequence, the implicit model _1 is “adapted to the data” and we have empirically found that, for any sufficiently long sequence y_1, …, y_t, the evidence[Given a data set { y_1, …, y_t }, the evidence in favour of a model  is the joint probability density of y_1, ⋯, y_t conditional on , denoted (y_1:t|).] <cit.> in favour of _1 is greater than the evidence in favour of _0. We can show, for several examples, that this implicit adaptation to the data makes the NuPF robust to mismatches in the state equation of the SSM compared to conventional PFs. In particular, provided that the likelihoods are specified or calibrated reliably, we have found that NuPFs perform reliably under a certain amount of mismatch in the transition kernel of the SSM, while standard PFs degrade clearly in the same scenario.



In order to illustrate the contributions outlined above, we present the results of several computer experiments with both synthetic and real data. In the first example, we assess the performance of the NuPF when applied to a linear-Gaussian SSM. The aim of these computer simulations is to compare the estimation accuracy and the computational cost of the proposed scheme with several other competing algorithms, namely a standard BPF, a PF with optimal proposal function and a NuPF with proper weights. The fact that the underlying SSM is linear-Gaussian enables the computation of the optimal importance function (intractable in a general setting) and proper weights for the NuPF. We implement the latter scheme because of its similarity to standard nudging filters in the literature. This example shows that the NuPF suffers just from a slight performance degradation compared to the PF with optimal importance function or the NuPF with proper weights, while the latter two algorithms are computationally more demanding.  

The second and third examples are aimed at testing the robustness of the NuPF when there is a significant misspecification in the state equation of the SSM. This is helpful in real-world applications because practitioners often have more control over measurement systems, which determine the likelihood, than they have over the state dynamics. We present computer simulation results for a stochastic Lorenz 63 model and a maneuvering target tracking problem. 

In the fourth example, we present numerical results for a  stochastic Lorenz 96 model, in order to show how a relatively high-dimensional system can be tracked without a major increase of the computational effort compared to the standard BPF. For this set of computer simulations we have also compared the NuPF with the Ensemble Kalman filter (EnKF), which is the de facto choice for tackling this type of systems.

Let us remark that, for the two stochastic Lorenz systems, the Markov kernel in the SSM can be sampled in a relatively straightforward way, yet transition probability densities cannot be computed (as they involve a sequence of noise variables mapped by a composition of nonlinear functions). Therefore, computing proper weights for proposal functions other than the Markov kernel itself is, in general, not possible for these examples.

Finally, we demonstrate the practical use of the NuPF on a problem where a real dataset is used to fit a stochastic volatility model using either particle Markov chain Monte Carlo (pMCMC) <cit.> or nested particle filters <cit.>.






 §.§ Organisation


The paper is structured as follows. After a brief note about notation, we describe the SSMs of interest and the BPF in Section <ref>. Then in Section <ref>, we outline the general algorithm and the specific nudging schemes we propose to use within the PF. We prove a convergence result in Section <ref> which shows that the new algorithm has the same asymptotic convergence rate as the BPF. We also provide an alternative interpretation of the nudging operation that explains its robustness in scenarios where there is a mismatch between the observed data and the assumed SSM. We discuss the computer simulation experiments in Section <ref> and present results for real data in Section <ref>. Finally, we make some concluding remarks in Section <ref>. 







 §.§ Notation

We denote the set of real numbers as , while ^d = ×d⋯× is the space of d-dimensional real vectors. We denote the set of positive integers with  and the set of positive reals with _+. We represent the state space with 𝖷⊂^d_x and the observation space with 𝖸⊂^d_y.

In order to denote sequences, we use the shorthand notation x_i_1:i_2 = {x_i_1,…,x_i_2}. For sets of integers, we use [n] = {1,…,n}. The  p-norm of a vector x∈^d is defined by x_p = (x_1^p + ⋯ + x_d^p)^1/p. The L_p norm of a random variable z with probability density function (pdf) p(z) is denoted z_p = (∫ |z|^p p(z)  z)^1/p, for p≥ 1. The Gaussian (normal) probability distribution with mean m and covariance matrix C is denoted (m,C). We denote the identity matrix of dimension d with I_d.

The supremum norm of a real function φ:𝖷→ is denoted φ_∞ = sup_x ∈𝖷 |φ(x)|. A function is bounded if φ_∞ < ∞ and we indicate the space of real bounded functions 𝖷 as B(𝖷). The set of probability measures on  is denoted (), the Borel σ-algebra of subsets of 𝖷 is denoted (𝖷) and the integral of a function φ:𝖷 with respect to a measure μ on the measurable space (𝖷,(𝖷)) is denoted (φ,μ):=∫φμ. The unit Dirac delta measure located at x ∈^d is denoted δ_x( x). The Monte Carlo approximation of a measure μ constructed using N samples is denoted as μ^N. Given a Markov kernel τ(x'|x) and a measure π(x), we define the notation ξ(x') = τπ≜∫τ(x'|x) π(x).













§ BACKGROUND








 §.§ State space models


We consider SSMs of the form

    x_0    ∼π_0( x_0), 
    
    x_t | x_t-1   ∼τ_t( x_t|x_t-1), 
    
    y_t | x_t    ∼ g_t(y_t|x_t),    t ∈,

where x_t ∈𝖷 is the system state at time t, y_t ∈𝖸 is the t-th observation, the measure π_0 describes the prior probability distribution of the initial state, τ_t is a Markov transition kernel on 𝖷, and g_t(y_t|x_t) is the (possibly non-normalised) pdf of the observation y_t conditional on the state x_t. We assume the observation sequence { y_t }_t ∈_+ is arbitrary but fixed. Hence, it is convenient to think of the conditional pdf g_t as a likelihood function and we write g_t(x_t) := g_t(y_t|x_t) for conciseness.

We are interested in the sequence of posterior probability distributions of the states generated by the SSM. To be specific, at each time t=1, 2, ... we aim at computing (or, at least, approximating) the probability measure π_t which describes the probability distribution of the state x_t conditional on the observation of the sequence y_1:t. When it exists, we use π(x_t|y_1:t) to denote the pdf of x_t given y_1:t with respect to the Lebesgue measure, i.e., π_t( x_t) = π(x_t|y_1:t) x_t.

The measure π_t is often termed the optimal filter at time t. It is closely related to the probability measure ξ_t, which describes the probability distribution of the state x_t conditional on y_1:t-1 and it is, therefore, termed the predictive measure at time t. As for the case of the optimal filter, we use ξ(x_t|y_1:t-1) to denote the pdf, with respect to the Lebesgue measure, of x_t given y_1:t-1.








 §.§ Bootstrap particle filter


The BPF <cit.> is a recursive algorithm that produces successive Monte Carlo approximations of ξ_t and π_t for t=1, 2, .... The method can be outlined as shown in Algorithm <ref>. 



After an initialization stage, where a set of independent and identically distributed (i.i.d.) samples from the prior are drawn, it consists of three recursive steps which can be depicted as,

    π_t-1^N →_samplingξ_t^N →_weightingπ̃_t^N →_resamplingπ_t^N.

Given a Monte Carlo approximation π_t-1^N = 1/N∑_i=1^N δ_x_t-1^(i) computed at time t-1, the sampling step yields an approximation of the predictive measure ξ_t of the form 

    ξ_t^N = 1/N∑_i=1^N δ_x̅_t^(i)

by propagating the particles { x_t-1^(i)}_i=1^N via the Markov kernel τ_t(·|x_t-1^(i)). The observation y_t is assimilated via the importance weights w_t^(i)∝ g_t(x_t^(i)), to obtain the approximate filter

    π̃_t^N = ∑_i=1^N w_t^(i)δ_x̅_t^(i),

and the resampling step produces a set of un-weighted particles that completes the recursive loop and yields the approximation

    π_t^N = 1/N∑_i=1^N δ_x_t^(i).

The random measures ξ_t^N, π̃_t^N and π_t^N are commonly used to estimate a posteriori expectations conditional on the available observations. For example, if φ is a function 𝖷, then the expectation of the random variable φ(x_t) conditional on y_1:t-1 is [ φ(x_t) | y_1:t-1] = (φ,ξ_t). The latter integral can be approximated using ξ_t^N, namely,

    (φ,ξ_t)    = ∫φ(x_t) ξ_t( x_t) ≈ (φ,ξ_t^N) 
       = ∫φ(x_t)ξ_t^N( x_t) = 1/N∑_i=1^N φ(x̅_t^(i)).

Similarly, we can have estimators (φ,π̃_t^N) ≈ (φ,π_t) and (φ,π_t^N) ≈ (φ,π_t). Classical convergence results are usually proved for real bounded functions, e.g., if φ∈ B(𝖷) then 

    lim_N∞ | (φ,π_t) - (φ,π_t^N) | = 0

under mild assumptions; see <cit.> and references therein.

The BPF can be generalized by using arbitrary proposal pdf's q_t(x_t|x_t-1^(i),y_t), possibly observation-dependent, instead of the Markov kernel τ_t(·|x_t-1^(i)) in order to generate the particles {x̅_t^(i)}_i=1^N in the sampling step. This can lead to more efficient algorithms, but the weight computation has to account for the new proposal and we obtain <cit.>

    w_t^(i)∝g_t(x̅_t^(i)) τ_t(x̅_t^(i) | x_t^(i))/q_t(x̅_t^(i) | x_t-1^(i), y_t),

which can be more costly to evaluate. This issue is related to the nudged PF to be introduced in Section <ref> below, which can be interpreted as a scheme to choose a certain observation-dependent proposal q_t(x_t|x_t-1^(i),y_t). However, the new method does not require that the weights be computed as in (<ref>) in order to ensure convergence of the estimators.










§ NUDGED PARTICLE FILTER







 §.§ General algorithm


Compared to the standard BPF, the nudged particle filter (NuPF) incorporates one additional step right after the sampling of the particles {x̅_t^(i)}_i=1^N at time t. The schematic depiction of the BPF in (<ref>) now becomes

    π_t-1^N →_samplingξ_t^N →_nudgingξ̃_t^N →_weightingπ̃_t^N →_resamplingπ_t^N,

where the new nudging step intuitively consists in pushing a subset of the generated particles {x̅_t^(i)}_i=1^N towards regions of the state space 𝖷 where the likelihood function g_t(x) takes higher values.

When considered jointly, the sampling and nudging steps in (<ref>) can be seen as sampling from a proposal distribution which is obtained by modifying the kernel τ_t(·|x_t-1) in a way that depends on the observation y_t. Indeed, this is the classical view of nudging in the literature <cit.>. However, unlike in this classical approach, here the weighting step does not account for the effect of nudging. In the proposed NuPF, the weights are kept the same as in the original filter, w_t^(i)∝ g_t(x_t^(i)). In doing so, we save computations but, at the same time, introduce bias in the Monte Carlo estimators. One of the contributions of this paper is to show that this bias can be controlled using simple design rules for the nudging step, while practical performance can be improved at the same time.

In order to provide an explicit description of the NuPF, let us first state a definition for the nudging step.

  A nudging operator α_t^y_t: 𝖷→𝖷 associated with the likelihood function g_t(x) is a map such that

    if      x' = α_t^y_t(x)      then     g_t(x') ≥ g_t(x)

for every x,x'∈𝖷.


Intuitively, we define nudging herein as an operation that increases the likelihood. There are several ways in which this can be achieved and we discuss some examples in Sections <ref> and <ref>. The NuPF with nudging operator α_t^y_t:𝖷𝖷 is outlined in Algorithm <ref>. 





It can be seen that the nudging operation is implemented in two stages.


  * First, we choose a set of indices ℐ_t ⊂ [N] that identifies the particles to be nudged. Let M=|ℐ_t| denote the number of elements in ℐ_t. We prove in Section <ref> that keeping M ≤𝒪(√(N)) allows the NuPF to converge with the same error rates 𝒪(1/√(N)) as the BPF. In Section <ref> we discuss two simple methods to build ℐ_t in practice.

  * Second, we choose an operator α_t^y_t that guarantees an increase of the likelihood of any particle. We discuss different implementations of α_t^y_t in Section <ref>.

We devote the rest of this section to a discussion of how these two steps can be implemented (in several ways).






 §.§ Selection of particles to be nudged
 

The set of indices ℐ_t, that identifies the particles to be nudged in Algorithm <ref>, can be constructed in several different ways, either random or deterministic. In this paper, we describe two simple random procedures with little computational overhead.  








  * Batch nudging: Let the number of nudged particles M be fixed. A simple way to construct ℐ_t is to draw indices i_1, i_2, …, i_M uniformly from [N] without replacement, and then let ℐ_t = i_1:M. We refer to this scheme as batch nudging, referring to selection of the indices at once. One advantage of this scheme is that the number of particles to be nudged, M, is deterministic and can be set a priori.  




  * Independent nudging: The size and the elements of ℐ_t can also be selected randomly in a number of ways. Here, we have studied a procedure in which, for each index i = 1, ..., N, we assign i ∈ℐ_t with probability M/N. In this way, the actual cardinality |ℐ_t| is random, but its expected value is exactly M. This procedure is particularly suitable for parallel implementations, since each index can be assigned to ℐ_t (or not) at the same time as all others.



























 §.§ How to nudge
 

The nudging step is aimed at increasing the likelihood of a subset of individual particles, namely those with indices contained in ℐ_t. Therefore, any map α_t^y_t:𝖷𝖷 such that (g_t ∘α_t^y_t)(x) ≥ g_t(x) when x∈𝖷 is a valid nudging operator. Typical procedures used for optimisation, such as gradient moves or random search schemes, can be easily adapted to implement (relatively) inexpensive nudging steps. Here we briefly describe a few of such techniques.





  * Gradient nudging: If g_t(x_t) is a differentiable function of x_t, one straightforward way to nudge particles is to take gradient steps. In Algorithm <ref> we show a simple procedure with one gradient step alone, and where γ_t>0 is a step-size parameter and ∇_x g_t(x) denotes the vector of partial derivatives of g_t with respect to the state variables, i.e., 

    ∇_x_t g_t = [
    	[   ∂ g_t/∂ x_1,t;   ∂ g_t/∂ x_2,t;               ⋮; ∂ g_t/∂ x_d_x,t;                 ]]      x_t = [
    	[   x_1,t;   x_2,t;       ⋮; x_d_x,t;         ]] ∈𝖷.

Algorithms can obviously be designed where nudging involves several gradient steps. In this work we limit our study to the single-step case, which is shown to be effective and keeps the computational overhead to a minimum. We also note that the performance of gradient nudging can be sensitive to the choice of the step-size parameters γ_t>0, which are, in turn, model dependent[We have found, nevertheless, that fixed step-sizes (i.e., γ_t=γ for all t) work well in practice for the examples of Sections <ref> and <ref>.]. 


  * Random nudging: Gradient-free techniques inherited from the field of global optimisation can also be employed in order to push particles towards regions where they have higher likelihoods. A simple stochastic-search technique adapted to the nudging framework is shown in Algorithm <ref>. We hereafter refer to the latter scheme as random-search nudging.




  * Model specific nudging: Particles can also be nudged using the specific model information. For instance, in some applications the state vector x_t can be split into two subvectors, x_t^obs and x_t^unobs (observed and unobserved, respectively), such that g_t(x_t) = g_t(x_t^obs), i.e., the likelihood depends only on x_t^obs and not on x_t^unobs. If the relationship between x_t^obs and x_t^unobs is tractable, one can first nudge x_t^obs in order to increase the likelihood and then modify x_t^unobs in order to keep it coherent with x_t^obs. A typical example of this kind arises in object tracking problems, where positions and velocities have a special and simple physical relationship but usually only position variables are observed through a linear or nonlinear transformation. In this case, nudging would only affect the position variables. However, using these position variables, one can also nudge velocity variables with simple rules. We discuss this idea and show numerical results in Section <ref>.















 §.§ Nudging general particle filters


In this paper we limit our presentation to BPFs in order to focus on the key concepts of nudging and to ease presentation. It should be apparent, however, that nudging steps can be plugged into general PFs. More specifically, since the nudging step is algorithmically detached from the sampling and weighting steps, it can be easily used within any PF, even if it relies on different proposals and different weighting schemes. We leave for future work the investigation of the performance of nudging within widely used PFs, such as auxiliary particle filters (APFs) <cit.>. 









§ ANALYSIS


The nudging step modifies the random generation of particles in a way that is not compensated by the importance weights. Therefore, we can expect nudging to introduce bias in the resulting estimators in general. However, in Section <ref> we prove that, as long as some basic guidelines are followed, the estimators of integrals with respect to the filtering measure π_t and the predictive measure ξ_t converge in L_p as N∞ with the usual Monte Carlo rate 𝒪(1/√(N)). The analysis is based on a simple induction argument and ensures the consistency of a broad class of estimators. In Section <ref> we briefly comment on the conditions needed to guarantee that convergence is attained uniformly over time. We do not provide a full proof, but this can be done by extending the classical arguments in <cit.> or <cit.> and using the same treatment of the nudging step as in the induction proof of Section <ref>. Finally, in Section <ref>, we provide an interpretation of nudging in a scenario with modelling errors. In particular, we show that the NuPF can be seen as a standard BPF for a modified dynamical model which is “a better fit” for the available data than the original SSM.











 §.§ Convergence in L_p
 

The goal in this section is to provide theoretical guarantees of convergence for the NuPF under mild assumptions. First, we analyze a general NuPF (with arbitrary nudging operator α_t^y_t and an upper bound on the size M of the index set ℐ_t) and then we provide a result for a NuPF with gradient nudging.

Before proceeding with the analysis, let us note that the NuPF produces several approximate measures, depending on the set of particles (and weights) used to construct them. After the sampling step, we have the random probability measure 

    ξ_t^N = 1/N∑_i=1^N δ_x̅_t^(i),

which converts into

    ξ̃_t^N = 1/N∑_i=1^N δ_x̃_t^(i)

after nudging. Once the weights w_t^(i) are computed, we obtain the approximate filter 

    π̃_t^N = ∑_i=1^N w_t^(i)δ_x̃_t^(i),

which finally yields

    π_t^N =  1/N∑_i=1^N δ_x_t^(i)
 
after the resampling step.

Similar to the BPF, the simple Assumption <ref> stated next is sufficient for consistency and to obtain explicit error rates <cit.> for the NuPF, as stated in Theorem <ref> below.


 The likelihood function is positive and bounded, i.e.,

    g_t(x_t) > 0     and    g_t_∞ = sup_x_t∈𝖷 |g_t(x_t)| < ∞

for t = 1,…,T.




 
Let y_1:T be an arbitrary but fixed sequence of observations, with T<∞, and choose any M ≤√(N) and any map α_t^y_t:→. If Assumption <ref> is satisfied and |ℐ_t|=M, then 

    (φ,π_t^N) - (φ, π_t) _p ≤c_t,pφ_∞/√(N)

for every t=1, 2, ..., T, any φ∈ B(𝖷), any p ≥ 1 and some constant c_t,p < ∞ independent of N.



See Appendix <ref> for a proof.

Theorem <ref> is very general; it actually holds for any map α_t^y_t: →, i.e., not necessarily a nudging operator. We can also obtain error rates for specific choices of the nudging scheme. A simple, yet practically appealing, setup is the combination of batch and gradient nudging, as described in Sections <ref> and <ref>, respectively. 


  
The gradient of the likelihood is bounded. In particular, there are constants G_t<∞ such that

    ∇_x g_t(x)_2 ≤ G_t    < ∞

for every x ∈𝖷 and t=1, 2, …, T.




 
Choose the number of nudged particles, M >0, and a sequence of step-sizes, γ_t>0, in such a way that sup_1 ≤ t ≤ Tγ_t M ≤√(N) for some T<0. If Assumption <ref> holds and φ is a Lipschitz test function, then the error introduced by the batch gradient nudging step with |ℐ_t|=M can be bounded as,

    (φ, ξ_t^N) - (φ, ξ̃_t^N) _p ≤LG_t/√(N),

where L is the Lipschitz constant of φ, for every t=1, …, T.


See Appendix <ref> for a proof. 

It is straightforward to apply Lemma <ref> to prove convergence of the NuPF with a batch gradient-nudging step. Specifically, we have the following result.


 
Let y_1:T be an arbitrary but fixed sequence of observations, with T<∞, and choose a sequence of step sizes γ_t>0 and an integer M such that 
    sup_1≤ t≤ Tγ_t M ≤√(N).
 Let π_t^N denote the filter approximation obtained with a NuPF with batch gradient nudging. If Assumptions <ref> and <ref> are satisfied and |ℐ_t|=M, then 

    (φ,π_t^N) - (φ, π_t) _p ≤c_t,pφ_∞/√(N)

for every t=1, 2, ..., T, any bounded Lipschitz function φ, some constant c_t,p < ∞ independent of N for any integer p ≥ 1.



The proof is straightforward (using the same argument as in the proof of Theorem <ref> combined with Lemma <ref>) and we omit it here. We note that Lemma <ref> provides a guideline for the choice of M and γ_t. In particular, one can select M = N^β, where 0 < β < 1, together with γ_t ≤ N^1/2 - β in order to ensure that γ_t M ≤√(N). Actually, it would be sufficient to set γ_t ≤ C N^1/2 - β for some constant C<∞ in order to keep the same error rate (albeit with a different constant in the numerator of the bound). Therefore,  Lemma <ref> provides a heuristic to balance the step size with the number of nudged particles[Note that the step sizes may have to be kept small enough to ensure that g_t(x̅_t^(i) + γ_t ∇_x g_t(x̅_t^(i)) ) ≥ g_t(x̅_t^(i)), so that proper nudging, according to Definition <ref>, is performed.]. We can increase the number of nudged particles but in that case we need to shrink the step size accordingly, so as to keep γ_t M ≤√(N). Similar results can be obtained using the gradient of the log-likelihood, log g_t, if the g_t comes from the exponential family of densities.







 §.§ Uniform convergence
 
Uniform convergence can be proved for the NuPF under the same standard assumptions as for the conventional BPF; see, e.g., <cit.>. The latter can be summarised as follows <cit.>: 


  (i) The likelihood function is bounded and bounded away from zero, i.e., g_t ∈ B(𝖷) and there is some constant a>0 such that inf_t>0, x∈𝖷 g_t(x) ≥ a.

  (ii) The kernel mixes sufficiently well, namely,  for any given integer m there is a constant 0<ε<1 such that

    inf_t>0;(x,x')∈𝖷^2τ_t+m|t(A|x)
    /τ_t+m|t(A|x')
     > ε

for any Borel set A, where τ_t+m|t is the composition of the kernels τ_t+m∘τ_t+m-1∘⋯∘τ_t.

When (i) and (ii) above hold, the sequence of optimal filters {π_t }_t≥ 0 is stable and it can be proved that 

    sup_t>0 (φ,π_t) - (φ,π_t^N) _p ≤c_p/√(N)

for any bounded function φ∈ B(𝖷), where c_p<∞ is constant with respect to N and t and π_t^N is the particle approximation produced by either the NuPF (as in Theorem <ref> or, provided sup_t>0 G_t < ∞, as in Theorem <ref>) or the BPF algorithms. We skip a formal proof as, again, it is straightforward combination of the standard argument by <cit.> (see also, e.g., <cit.> and <cit.>) with the same handling of the nudging operator as in the proofs of Theorem <ref> or Lemma <ref> .







 §.§ Nudging as a modified dynamical model
 
We have found in computer simulation experiments that the NuPF is consistently more robust to model errors than the conventional BPF. In order to obtain some analytical insight of this scenario, in this section we reinterpret the NuPF as a standard BPF for a modified, observation-driven dynamical model and discuss why this modified model can be expected to be a better fit for the given data than the original SSM. In this way, the NuPF can be seen as an automatic adaptation of the underlying model to the available data.

The dynamic models of interest in stochastic filtering can be defined by a prior measure τ_0, the transition kernels τ_t and the likelihood functions g_t(x)=g_t(y_t|x), for t≥ 1. In this section we write the latter as g_t^y_t(x) = g_t(y_t|x), in order to emphasize that g_t is parametrised by the observation y_t, and we also assume that every g_t^y_t is a normalised pdf in y_t for the sake of clarity. Hence, we can formally represent the SSM defined by (<ref>), (<ref>) and (<ref>) as _0={τ_0,τ_t,g_t^y_t}.







Now, let us assume y_1:T to be fixed and construct the alternative dynamical model _1 = {τ_0, τ̃_t^y_t, g_t^y_t}, where

    τ̃_t^y_t( x_t | x_t-1 ) :=    (1-ε_M)τ_t( x_t|x_t-1)   + 
       ε_M ∫δ_α_t^y_t(x̅_t) ( x_t) τ_t(x̅_t | x_t-1)

is an observation-driven transition kernel, ε_M = M/N and the nudging operator α_t^y_t is a one-to-one map that depends on the (fixed) observation y_t. We note that the kernel τ̃_t^y_t jointly represents the Markov transition induced by the original kernel τ_t followed by an independent nudging transformation (namely, each particle is independently nudged with probability ε_M). As a consequence, the standard BPF for model _1 coincides exactly with a NuPF for model _0 with independent nudging and operator α_t^y_t. Indeed, according to the definition of τ̃_t^y_t in (<ref>), generating a sample x̃_t^(i) from τ̃_t^y_t( x_t | x_t-1^(i)) is a three-step process where 


  * we first draw x̅_t^(i) from τ_t( x_t|x_t-1^(i)), 

  * then generate a sample u_t^(i) from the uniform distribution (0,1), and 

  * if u_t^(i) < ε_M then we set x̃_t^(i) = α_t^y_t(x̅_t^(i)), else we set x̃_t^(i) = x̅_t^(i).

After sampling, the importance weight for the BPF applied to model _1 is w_t^(i)∝ g_t^y_t(x̃_t^(i)). This is exactly the same procedure as in the NuPF applied to the original SSM _0 (see Algorithm <ref>).

Intuitively, one can expect that the observation-driven _1 is a better fit for the data sequence y_1:T than the original model _0. Within the Bayesian methodology, a common approach to compare two competing probabilistic models (_0 and _1 in this case) for a given data set y_1:t is to evaluate the so-called model evidence <cit.> for both _0 and _1. 

The evidence (or likelihood) of a probabilistic model  for a given data set y_1:t is the probability density of the data conditional on the model, that we denote as (y_1:t|).

We say that _1 is a better fit than _0 for the data set y_1:t when (y_1:t|_1) > (y_1:t|_0). Since

    (y_1:t|_0) = ∫⋯∫∏_l=1^t g_l(x_l)τ_l( x_l | x_l-1) τ_0( x_0),

and

    (y_1:t|_1) = ∫⋯∫∏_l=1^t g_l(x_l)τ̃_l^y_l( x_l | x_l-1) τ_0( x_0),

the difference between the evidence of _0 and the evidence of _1 depends on the difference between the transition kernels τ̃_t and τ̃_t^y_t.

We have empirically observed in several computer experiments that (y_1:t|_1) > (y_1:t|_0) and we argue that the observation-driven kernel τ̃_t^y_t implicit to the NuPF is the reason why the latter filter is robust to modelling errors in the state equation, compared to standard PFs. This claim is supported by the numerical results in Sections <ref> and <ref>, which show how the NuPF attains a significant better performance than the standard BPF, the auxiliary PF <cit.> or the extended Kalman filter <cit.> in scenarios where the filters are built upon a transition kernel different from the one used to generate the actual observations.


While it is hard to show that (y_1:t|_1) > (y_1:t|_0) for every NuPF, it is indeed possible to guarantee that the latter inequality holds for specific nudging schemes. An example is provided in Appendix <ref>, where we describe a certain nudging operator α_t^y_t and then proceed to prove that (y_1:t|_1) > (y_1:t|_0), for that particular scheme, under some regularity conditions on the likelihoods and transition kernels.







§ COMPUTER SIMULATIONS


In this section, we present the results of several computer experiments. In the first one, we address the tracking of a linear-Gaussian system. This is a very simple model which enables a clearcut comparison of the NuPF and other competing schemes, including a conventional PF with optimal importance function (which is intractable for all other examples) and a PF with nudging and proper importance weights. Then, we study three nonlinear tracking problems: 


  * a stochastic Lorenz 63 model with misspecified parameters, 

  * a maneuvering target monitored by a network of sensors collecting nonlinear observations corrupted with heavy-tailed noise, 

  * and, finally,  a high-dimensional stochastic Lorenz 96 model[For the experiments involving Lorenz 96 model, simulation from the model is implemented in C++ and integrated into Matlab. The rest of the simulations are fully implemented in Matlab.].


We have used gradient nudging in all experiments, with either M ≤√(N) (deterministically, with batch nudging) or [M]≤√(N) (with independent nudging). This ensures that the assumptions of Theorem <ref> hold. For simplicity, the gradient steps are computed with fixed step sizes, i.e., γ_t=γ for all t. For the object tracking experiment, we have used a large step-size, but this choice does not affect the convergence rate of the NuPF algorithm either.








 §.§ A high-dimensional, inhomogeneous Linear-Gaussian state-space model


In this experiment we compare different PFs implemented to track a high-dimensional linear Gaussian SSM. In particular, the model under consideration is

    x_0    ∼(0,I_d_x), 
    
    x_t|x_t-1   ∼(x_t-1,Q),
    
    y_t|x_t    ∼(C_t x_t, R),

where {x_t}_t≥ 0 are hidden states, {y_t}_t≥ 1 are observations, and Q and R are the process and the observation noise covariance matrices, respectively. The latter are diagonal matrices, namely Q = q I_d_x and R = I_d_y, where q = 0.1, d_x = 100 and d_y = 20. The sequence {C_t}_t≥ 1 defines a time-varying observation model. The elements of this sequence are chosen as random binary matrices, i.e., C_t ∈{0,1}^d_y × d_x where each entry is generated as an independent Bernoulli random variable with p = 0.5. Once generated, they are fixed and fed into all algorithms we describe below for each independent Monte Carlo run.

We compare the NuPF with three alternative PFs. The first method we implement is the PF with the optimal proposal pdf p(x_t|x_t-1,y_t) ∝ g_t(y_t|x_t) τ_t(x_t|x_t-1), abbreviated as Optimal PF. The pdf p(x_t|x_t-1,y_t) leads to an analytically tractable Gaussian density for the model (<ref>)–(<ref>) <cit.> but not in the nonlinear tracking examples below. Note, however, that at each time step, the mean and covariance matrix of this proposal have to be explicitly evaluated in order to compute the importance weights. 

The second filter is a nudged PF with proper importance weights (NuPF-PW). In this case, we treat the generation of the nudged particles as a proposal function to be accounted for during the weighting step. To be specific, the proposal distribution resulting from the NuPF has the form

    τ̃_t(x_t | x_t-1) = (1-ϵ_N) τ_t(x_t|x_t-1) + ϵ_N τ̅_t(x_t|x_t-1),

where ϵ_N = 1/√(N) and

    τ̅_t(x_t|x_t-1) = ∫δ_α_t^y_t(x̅_t) (x_t) τ_t(x̅_t|x_t-1).

The latter conditional distribution admits an explicit representation as a Gaussian for model (<ref>)-(<ref>) when the α_t operator is designed as a gradient step, but this approach is intractable for the examples in Section <ref> and Section <ref>. Note that τ̃_t is a mixture of two time-varying Gaussians and this fact adds to the cost of the sampling and weighting steps. Specifically, computing weights for the NuPF-PW is significantly more costly, compared to the BPF or the NuPF, because the mixture (<ref>) has to be evaluated together with the likelihood and the transition pdf. 

The third tracking algorithm implemented for model (<ref>)–(<ref>) is the conventional BPF.

For all filters, we have set the number of particles as[When N is increased the results are similar for the NuPF, the optimal PF and the NuPF-PW larger number particles, as they already perform close to optimally for N=100, and only the BPF improves significantly.] N = 100 . In order to implement the NuPF and NuPF-PW schemes, we have selected the step size γ = 2 × 10^-2. We have run 1,000 independent Monte Carlo runs for this experiment. To evaluate different methods, we have computed the empirical normalised mean squared errors (NMSEs). Specifically, the NMSE for the j-th simulation is 

    (j) =  ∑_t=1^t_fx̅_t - x̂_t(j) _2^2/∑_t=1^t_f x_t _2^2,

where x̅_t = 𝔼[x_t|y_1:t] is the exact posterior mean of the state x_t conditioned on the observations up to time t and x̂_t(j) is the estimate of the state vector in the j-th simulation run. Therefore, the notation  implies the normalised mean squared error is computed with respect to x̅_t. In the figures, we usually plot the mean and the standard deviation of the sample of errors, (1), …, (1000).

The results  are shown in Fig. <ref>. In particular, in Fig. <ref>(a), we observe that the  performance of the NuPF compared to the optimal PF and NuPF-PW (which is similar to a classical PF with nudging) is comparable. However, Fig. <ref>(b) reveals that the NuPF is significantly less demanding compared to the optimal PF and the NuPF-PW method. Indeed, the run-times of the NuPF are almost identical to the those of the plain BPF. As a result, the plot of the s multiplied by the running times displayed in Fig. <ref>(b) reveals that the proposed algorithm is as favorable as the optimal PF, which can be implemented for this model, but not for general models unlike the NuPF.








 §.§ Stochastic Lorenz 63 model with misspecified parameters
 

In this experiment, we demonstrate the performance of the NuPF when tracking a misspecified stochastic Lorenz 63 model. The dynamics of the system is described by a stochastic differential equation (SDE) in three dimensions,

    x_1    = - 𝖺(x_1 - x_2) s +  w_1, 
     x_2    = ( 𝗋x_1 - x_2 - x_1 x_3 )  s +  w_2, 
     x_3    = ( x_1 x_2 - 𝖻x_3 )  s +  w_3,

where s denotes continuous time, {w_i(s)}_s∈(0,∞) for i = 1,2,3 are 1-dimensional independent Wiener processes and 𝖺,𝗋,𝖻∈ are fixed model parameters. We discretise the model using the Euler-Maruyama scheme with integration step 𝖳 > 0 and obtain the system of difference equations

    x_1,t   =    x_1,t-1 - 𝖳𝖺 (x_1,t-1 - x_2,t-1) 
           + √(𝖳) u_1,t, 
    
    x_2,t   =    x_2,t-1 + 𝖳 (𝗋 x_1,t-1 -
    x_2,t-1 - x_1,t-1 x_3,t-1) 
           + √(𝖳) u_2,t, 
    
    x_3,t   =    x_3,t-1 + 𝖳 (x_1,t-1 x_2,t-1 - 𝖻x_3,t-1) + √(𝖳) u_3,t,

where {u_i,t}_t∈, i = 1,2,3 are i.i.d. Gaussian random variables with zero mean and unit variance. We assume that we can only observe the variable x_1,t, contaminated by additive noise, every t_s>1 discrete time steps. To be specific, we collect the sequence of observations

    y_n = k_o x_1,n t_s + v_n,    n = 1, 2, ...,

where {v_n}_n∈ is a sequence of i.i.d. Gaussian random variables with zero mean and unit variance and the scale parameter k_o = 0.8 is assumed known. 

In order to simulate both the state signal and the synthetic observations from this model, we choose the so-called standard parameter values 

    (𝖺,𝗋,𝖻) = (10,28,8/3),

which make the system dynamics chaotic. The initial condition is set as 

    x_0 = [-5.91652, -5.52332, 24.5723]^⊤.

The latter value has been chosen from a deterministic trajectory of the system (i.e., with no state noise) with the same parameter set  (𝖺,𝗋,𝖻) = (10,28,8/3) to ensure that the model is started at a sensible point. We assume that the system is observed every t_s = 40 discrete time steps and for each simulation we simulate the system for t=0, 1, …, t_f, with t_f = 20,000. Since t_s = 40, we have a sequence of t_f/t_s=500 observations overall. 

Let us note here that the Markov kernel which takes the state from time n-1 to time n (i.e., from the time of one observation to the time of the next observation) is straightforward to simulate using the Euler-Maruyama scheme (<ref>), however the associated transition probability density cannot be evaluated because it involves the mapping of both the state and a sequence of t_s noise samples through a composition of nonlinear functions. This precludes the use of importance sampling schemes that require the evaluation of this density when computing the weights.

We run the BPF and NuPF algorithms for the model described above, except that the parameter 𝖻 is replaced by 𝖻_ϵ = 𝖻+ϵ, with ϵ = 0.75 (hence 𝖻_ϵ≈ 3.417 versus 𝖻≈ 2.667 for the actual system). As the system underlying dynamics is chaotic, this mismatch affects the predictability of the system significantly.



We have implemented the NuPF with independent gradient nudging. Each particle is nudged with probability 1/√(N), where N is the number of particles (hence [M]=√(N)), and the size of the gradient steps is set to γ = 0.75 (see Algorithm <ref>). 

As a figure of merit, we evaluate the NMSE for the 3-dimensional state vector, averaged over 1,000 independent Monte Carlo simulations. For this example (as well as in the rest of this section), it is not possible to compute the exact posterior mean of the state variables. Therefore, the NMSE values are computed with respect to the ground truth, i.e.,

    (j) =  ∑_t=1^t_fx_t - x̂_t(j) _2^2/∑_t=1^t_f x_t _2^2,

where (x_t)_t≥ 1 is the ground truth signal.

Fig. <ref> (a) displays the , attained for varying number of particles N, for the standard BPF and the NuPF. It is seen that the NuPF outperforms the BPF for the whole range of values of N in the experiment, both in terms of the mean and the standard deviation of the errors, although the NMSE values become closer for larger N. The plot on the right displays the values of x_2,t and its estimates for a typical simulation. In general, the experiment shows that the NuPF can track the actual system using the misspecified model and a small number of particles, whereas the BPF requires a higher computational effort to attain a similar performance.

As a final experiment with this model, we have tested the robustness of the algorithms with respect to the choice of parameters in the nudging step. In particular, we have tested the NuPF with independent gradient nudging for a wide range of step-sizes γ. Also, we have tested the NuPF with random search nudging using a wide range of covariances of the form C = σ^2 I by varying σ^2.

The results can be seen in Fig. <ref>. This figure shows that the algorithm is robust to the choice of parameters for a range of step-sizes and variances of the random search step. As expected, random search nudging takes longer running time compared to gradient steps. This difference in run-times is expected to be larger in higher-dimensional models since random search is expected to be harder in such scenarios.









 §.§ Object tracking with a misspecified model
 

In this experiment, we consider a tracking scenario where a target is observed through sensors collecting radio signal strength (RSS) measurements contaminated with additive heavy-tailed noise. The target dynamics are described by the model,

    x_t    = A x_t-1 + B L (x_t-1 - x_∙) + u_t,

where x_t∈^4 denotes the target state, consisting of its position r_t ∈^2 and its velocity, v_t ∈^2, hence x_t = [ [ r_t; v_t;     ]]∈^4. The vector x_∙ is a deterministic, pre-chosen state to be attained by the target. Each element in the sequence {u_t}_t∈ is a zero-mean Gaussian random vector with covariance matrix Q. The parameters A,B,Q are selected as

    A = [      I_2    κ I_2;        0 0.99 I_2 ],   
    B = [   0 I_2 ]^⊤,

and

    Q = [ κ^3/3 I_2 κ^2/2 I_2; κ^2/2 I_2     κ I_2 ],

where I_2 is the 2× 2 identity matrix and κ = 0.04. The policy matrix L ∈^2×4 determines the trajectory of the target from an initial position x_0 = [140,140,50,0]^⊤ to a final state x_T = [140,-140,0,0]^⊤ and it is computed by solving a Riccati equation (see <cit.> for details), which yields

    L = [
    	[ -0.0134       0 -0.0381       0;       0 -0.0134       0 -0.0381;         ]].

This policy results in a highly maneuvering trajectory. In order to design the NuPF, however, we assume the simpler dynamical model 

    x_t    = A x_t-1 + u_t,

hence there is a considerable model mismatch.

The observations are nonlinear and coming from 10 sensors placed in the region where the target moves. The measurement collected at the i-th sensor, time t, is modelled as 

    y_t,i = 10 log_10(P_0/r_t - s_i^2 + η) + w_t,i,

where r_t ∈^2 is the location vector of the target, s_i is the position of the ith sensor and w_t,i∼𝒯(0,1,ν) is an independent t-distributed random variable for each i = 1,…,10. Intuitively, the closer the parameter ν to 1, the more explosive the observations become. In particular, we set ν = 1.01 to make the observations explosive and heavy-tailed. As for the sensor parameters, we set the transmitted RSS as P_0 = 1 and the sensitivity parameter as η = 10^-9. The latter yields a soft lower bound of -90 decibels (dB) for the RSS measurements.

We have implemented the NuPF with batch gradient nudging, with a large-step size γ = 5.5 and M = ⌊√(N)⌋. Since the observations depend on the position vector r_t only, an additional model-specific nudging step is needed for the velocity vector v_t. In particular, after nudging the r_t^(i) = [ x^(i)_1,t, x^(i)_2,t]^⊤, we update the velocity variables as

    v_t^(i) = 1/κ (r^(i)_t - r^(i)_t-1),      v_t^(i) =   [ x^(i)_3,t, x^(i)_4,t]^⊤,

where κ= 0.04 as defined for the model. The motivation for this additional transformation comes from the physical relationship between position and velocity. We note, however, that the NuPF also works without nudging the velocities.

We have run 10,000 Monte Carlo runs with N = 500 particles in the auxiliary particle filter (APF) <cit.>, the BPF <cit.> and the NuPF. We have also implemented the extended Kalman filter (EKF), which uses the gradient of the observation model.


Fig. <ref> shows a typical simulation run with each one of the four algorithms (on the left side, plots (a)–(d)) and a box-plot of the NMSEs obtained for the 10,000 simulations (on the right, plot (e)). Plots (a)–(d) show that, while the EKF also uses the gradient of the observation model, it fails to handle the heavy-tailed noise, as it relies on Gaussian approximations. The BPF and the APF collapse due to the model mismatch in the state equation. Plot (d) shows that the NMSE of the NuPF is just slightly smaller in the mean than the NMSE of the EKF, but much more stable.








 §.§ High-dimensional stochastic Lorenz 96 model
 

In this computer experiment, we compare the NuPF with the ensemble Kalman filter (EnKF) for the tracking of a stochastic Lorenz 96 system. The latter is described by the set of stochastic differential equations (SDEs)

    x_i    =   [ (x_i+1 - x_i-2) x_i-1 - x_i + F ]  s +  w_i, 
           i = 1, …, d,

where s denotes continuous time, { w_i(s) }_s ∈ (0,∞), 1 ≤ i ≤ d, are independent Wiener processes, d is the system dimension and the forcing parameter is set to F = 8, which ensures a chaotic regime. The model is assumed to have a circular spatial structure, so that x_-1 = x_d-1, x_0 = x_d, and x_d+1 = x_1. Note that each x_i, i=1, ..., d, denotes a time-varying state associated to a different space location. In order to simulate data from this model, we apply the Euler-Maruyama discretisation scheme and obtain the difference equations,

    x_i,t   =    x_i,t-1 + 𝖳[ (x_i+1,t-1 - x_i-2,t-1) x_i-1,t-1. 
          . - x_i,t-1 + F ] + √(𝖳) u_i,t,

where u_i,t are zero-mean, unit-variance Gaussian random variables. We initialise this system by generating a vector from the uniform distribution on (0,1)^d and running the system for a small number of iterations and set x_0 as the output of this short run.




We assume that the system is only partially observed. In particular, half of the state variables are observed, in Gaussian noise, every t_s = 10 time steps, namely,

    y_j,n = x_2j-1,nt_s + u_j,n,

where n = 1, 2, …, j = 1, 2, …, ⌊ d/2 ⌋, and u_j,n is a normal random variable with zero mean and unit variance. The same as in the stochastic Lorenz 63 example of Section <ref>, the transition pdf that takes the state from time (n-1)t_s to time nt_s is simple to simulate but hard to evaluate, since it involves mapping a sequence of noise variables through a composition of nonlinearities.

In all the simulations for this system we run the NuPF with batch gradient nudging (with M = ⌊√(N)⌋ nudged particles and step-size γ = 0.075). In the first computer experiment, we fixed the dimension d = 40 and run the BPF and the NuPF with increasing number of particles. The results can be seen in Fig. <ref>, which shows how the NuPF performs better than the BPF in terms of NMSE (plot (a)). Since the run-times of both algorithms are nearly identical, it can be seen that, when considered jointly with NMSEs, the NuPF attains a significantly better performance (plot (b)).

In a second computer experiment, we compared the NuPF with the EnKF. Fig. <ref>(a) shows how the NMSE of the two algorithms grows as the model dimension d increases and the number of particles N is kept fixed. In particular, the EnKF attains a better performance for smaller dimensions (up to d=10^3), however its NMSE blows up for d > 10^3 while the performance of the NuPF remains stable. The running time of the EnKF was also higher than the running time of the NuPF in the range of higher dimensions (d≥ 10^3).






 §.§ Assessment of bias


In this section, we numerically quantify the bias of the proposed algorithm on a low-dimensional linear-Gaussian state-space model. To assess the bias, we compute the marginal likelihood estimates given by the BPF and the NuPF. The reason for this choice is that the BPF is known to yield unbiased estimates of the marginal likelihood[Note that the estimates of integrals (φ,π_t) computed using the self-normalised importance sampling approximations (i.e., (φ,π_t^N) ≈ (φ,π_t)) produced by the BPF and the NuPF methods are biased and the bias vanishes with the same rate for both algorithms as a result of Theorem <ref>. The same is true for the approximate predictive measures ξ_t^N.] <cit.>. The NuPF leads to biased (typically overestimated) marginal likelihood estimates, hence it is of interest to compare them with those of the BPF. To this end, we choose a simple linear-Gaussian state space model for which the marginal likelihood can be exactly computed as a byproduct of the Kalman filter. We then compare this exact marginal likelihood to the estimates given by the BPF and the NuPF.

Particularly, we define the state-space model,

    x_0    ∼(x_0;μ_0,P_0),
    
    x_t | x_t-1   ∼(x_t; x_t-1,Q), 
    
    y_t | x_t    ∼(y_t; C_t x_t,R),

where (C_t)_t≥ 0∈ [0,1]^1 × 2 is a sequence of observation matrices where each entry is generated as a realisation of a Bernoulli random variable with p = 0.5, μ_0 is a zero vector, and x_t ∈^2 and y_t ∈. The state variables are cross-correlated, namely,

    Q = [   2.7 -0.48; -0.48  2.05 ],

and R = 1. We have chosen the prior covariance as P_0 = I_d_x. We have simulated the system for T=100 time steps. Given a fixed observation sequence y_1:T, the marginal likelihood for the system given in Eqs. (<ref>)-(<ref>) is

    Z^⋆ = 𝗉(y_1:T),

which can be exactly computed via the Kalman filter.


We denote the estimate of Z^⋆ given by the BPF and the NuPF as Z^N_BPF and Z^N_NuPF, respectively. It is well-known that the BPF estimator is unbiased <cit.>,

    [Z^N_BPF] = Z^⋆,

where [·] denotes the expectation with respect to the randomness of the particles. Numerically, this suggests that as one runs identical, independent Monte Carlo simulations to obtain {Z^N,k_BPF}_k=1^K and compute the average

    Z̅^N_BPF = 1/K∑_k=1^K Z^N,k_BPF,

then it follows from the unbiasedness property (<ref>) that the ratio of the average in (<ref>) and the true value Z^⋆ should satisfy

    Z̅^N_BPF/Z^⋆→ 1   as   K →∞.

Since the marginal likelihood estimates provided by the NuPF are not unbiased for the original SSM (and tend to attain higher values), if we define

    Z̅^N_NuPF = 1/K∑_k=1^K Z^N,k_NuPF,

then as K→∞, we should see

    Z̅^N_NuPF/Z^⋆→ 1+ϵ  as   K →∞,

for some ϵ > 0.

We have conducted an experiment aimed at quantifying the bias ϵ>0 above. In particular, we have run 20,000 independent simulations for the BPF and the NuPF with N=100, N=1,000 and N=10,000. For each value of N, we have computed running empirical means as in (<ref>) and (<ref>) for K=1, …, 20,000. The variance of Z̅_BPF^N increases with T, hence the estimators for small K display a relatively large variance and we need K>>1 to clearly observe the bias. The NuPF filter performs independent gradient nudging with step size γ = 0.1. 

The results of the experiment are displayed in Figure <ref>, which shows how, as expected, the NuPF overestimates Z^⋆. We can also see how the bias becomes smaller as N increases (because only and average of √(N) particles are nudged per time step).






§ EXPERIMENTAL RESULTS ON MODEL INFERENCE









In this section, we illustrate the application of the NuPF to estimate the parameters of a financial time-series model. In particular, we adopt a stochastic-volatility SSM and we aim at estimating its unknown parameters (and track its state variables) using the EURUSD log-return data from 2014-12-31 to 2016-12-31 (obtained from ). For this task, we apply two recently proposed Monte Carlo schemes: the nested particle filter (NPF) <cit.> (a purely recursive, particle-filter style Monte Carlo method) and the particle Metropolis-Hastings (pMH) algorithm <cit.> (a batch Markov chain Monte Carlo procedure). In their original forms, both algorithms use the marginal likelihood estimators given by the BPF to construct a Monte Carlo approximation of the posterior distribution of the unknown model parameters. Here, we compare the performance of these algorithms when the marginal likelihoods are computed using either the BPF or the proposed NuPF. 

We assume the stochastic volatility SSM <cit.>,

    x_0    ∼(μ, σ_v^2/1-ϕ^2), 
    
    x_t|x_t-1   ∼(μ + ϕ (x_t-1 - μ), σ_v^2), 
    
    y_t|x_t    ∼(0,exp(x_t)),

where μ∈, σ_v ∈_+, and ϕ∈ [-1,1] are fixed but unknown parameters. The states {x_t}_1 ≤ t ≤ T are log-volatilities and the observations {y_t}_1≤ t ≤ T are log-returns. We follow the same procedure as <cit.> to pre-process the observations. Given the historical price sequence s_0, …, s_T, the log-return at time t is calculated as

    y_t = 100 log(s_t/s_t-1)

for 1≤ t ≤ T. Then, given y_1:T, we tackle the joint Bayesian estimation of x_1:T and the unknown parameters θ = (μ,σ_v,ϕ). In the next two subsections we compare the conventional BPF and the NuPF as building blocks of the NPF and the pMH algorithms.



 §.§ Nudging the nested particle filter


The NPF in <cit.> consists of two layers of particle filters which are used to jointly approximate the posterior distributions of the parameters and the states. The filter in the first layer builds a particle approximation of the marginal posterior distribution of the parameters. Then, for each particle in the parameter space, say θ^(i), there is an inner filter that approximates the posterior distribution of the states conditional on the parameter vector θ^(i). The inner filters are classical particle filters, which are essentially used to compute the importance weights (marginal likelihoods) of the particles in the parameter space. In the implementation of <cit.>, the inner filters are conventional BPFs. We have compared this conventional implementation with an alternative one where the BPFs are replaced by the NuPFs. For a detailed description of the NPF, see <cit.>.

In order to assess the performances of the nudged and classical versions of the NPF, we compute the model evidence estimate given by the nested filter by integrating out both the parameters and the states. In particular, if the set of particles in the parameter space at time t is {θ_t^(i)}_i=1^K and for each particle θ_t^(i) we have a set of particles in the state space {x_t^(i,j)}_j=1^N, we compute 


    p(y_1:T) = ∏_t=1^T [ 1/KN∑_i=1^K ∑_j=1^N g_t(x_t^(i,j)) ].

The model evidence quantifies the fitness of the stochastic volatility model for the given dataset, hence we expect to see a higher value when a method attains a better performance (the intuition is that if we have better estimates of the parameters and the states, then the model will fit better). For this experiment, we compute the model evidence for the nudged NPF before the nudging step, so as to make the comparison with the conventional algorithm fair.

We have conducted 1,000 independent Monte Carlo runs for each algorithm and computed the model evidence estimates. We have used the same parameters and the same setup for the two versions of the NPF (nudged and conventional). In particular, each unknown parameter is jittered independently. The parameter μ is jittered with a zero-mean Gaussian kernel variance σ_μ^2 = 10^-3, the parameter σ_v is jittered with a truncated Gaussian kernel on (0,∞) with variance σ_σ_v^2 = 10^-4, and the parameter ϕ is jittered with a zero-mean truncated Gaussian kernel on [-1,1], with variance σ_ϕ^2 = 10^-4. We have chosen a large step-size for the nudging step, γ = 4,  and we have used batch nudging with M = ⌊√(N)⌋.

The results in Fig. <ref> demonstrate empirically that the use of the nudging step within the NPF reduces the variance of the model evidence estimators, hence it improves the numerical stability of the NPF.








 §.§ Nudging the particle Metropolis-Hastings


The pMH algorithm is a Markov chain Monte Carlo (MCMC) method for inferring parameters of general SSMs <cit.>. The pMH uses PFs as auxiliary devices to estimate parameter likelihoods in a similar way as the NPF uses them to compute importance weights. In the case of the pMH, these estimates should be unbiased and they are needed to determine the acceptance probability for each element of the Markov chain. For the details of the algorithm, see <cit.> (or <cit.> for a tutorial-style introduction). Let us note that the use of NuPF does not lead to an unbiased estimate of the likelihood with respect to the assumed SSM. However, as discussed in Section <ref>, one can view the use of nudging in this context as an implementation of pMH with an implicit dynamical model _1 derived from the original SSM _0.



We have carried out a computer experiment to compare the performance of the pMH scheme using either BPFs or NuPFs to compute acceptance probabilities. The two algorithms are labeled pMH-BPF and pMH-NuPF, respectively, hereafter. The parameter priors in the experiment are

    p(μ) = (0,1),          p(σ_v) = 𝒢(2,0.1),          p(ϕ) = ℬ(120,2),

where 𝒢(a,θ) denotes the Gamma pdf with shape parameter a and scale parameter θ, and ℬ(α,β) denotes the Beta pdf with shape parameters (α,β). Unlike <cit.>, who use a truncated Gaussian prior centered on 0.95 with a small variance for ϕ, we use the Beta pdf, which is defined on [0,1], with mean α/(α + β) = 0.9836, which puts a significant probability mass on the interval [0.9,1].

We have compared the pMH-BPF algorithm and the pMH-NuPF scheme (using a batch nudging procedure with γ = 0.1 and M = ⌊√(N)⌋) by running 1,000 independent Monte Carlo trials. We have computed the marginal likelihood estimates in the NuPF after the nudging step. 



First, in order to illustrate the impact of the nudging on the parameter posteriors, we have run the pMH-NuPF and the pMH-BPF and obtained a long Markov chain (2× 10^6 iterations) from both algorithms. Figure <ref> displays the two-dimensional marginals of the resulting posterior distribution. It can be seen from Fig. <ref> that the bias of the NuPF yields a perturbation compared to the posterior distribution approximated with the pMH-BPF. The discrepancy is small but noticeable for small N (see Fig. <ref>(a) for N=100) and vanishes as we increase N (see Fig. <ref>(b) and (c), for N=500 and N=1,000, respectively). We observe that for a moderate number of particles, such as N=500 in Fig. <ref>(b), the error in the posterior distribution due to the bias in the NuPF is very slight.

Two common figures of merit for MCMC algorithms are the acceptance rate of the Markov kernel (desirably high) and the autocorrelation function of the chain (desirably low). Figure <ref> shows the acceptance rates for the pMH-NuPF and the pMH-BPF algorithms with N=100, N=500 and N=1,000 particles in both PFs. It is observed that the use of nudging leads to noticeably higher acceptance rates, although the difference becomes smaller as N increases.

Figure <ref> displays the average autocorrelation functions (ACFs) of the chains obtained in the 1,000 independent simulations.We see that the autocorrelation of the chains produced by the pMH-NuPF method decays more quickly than the autocorrelation of the chains output by the conventional pMH-BPF, especially for lower values of N. Even for N=1,000 (which ensures an almost negligible perturbation of the posterior distribution, as shown in Figure <ref>(c)) there is an improvement in the ACFs of the parameters ϕ and σ_v using the NuPF.   Less correlation can be expected to translate into better estimates as well for a fixed length of the chain.






§ CONCLUSIONS


We have proposed a simple modification of the particle filter which, according to our computer experiments, can improve the performance of the algorithm (e.g., when tracking high-dimensional systems) or enhance its robustness to model mismatches in the state equation of a SSM. The modification of the standard particle filtering scheme consists of an additional step, which we term nudging, in which a subset of particles are pushed towards regions of the state space with a higher likelihood. In this way, the state space can be explored more efficiently while keeping the computational effort at nearly the same level as in a standard particle filter. We refer to the new algorithm as the “nudged particle filter” (NuPF). While, for clarity and simplicity, we have kept the discussion and the numerical comparisons restricted to the modification (nudging) of the conventional BPF, the new step can be naturally incorporated to most known particle filtering methods. 

We have presented a basic analysis of the NuPF which indicates that the algorithm converges (in L_p) with the same error rate as the standard particle filter. In addition, we have also provided a simple reinterpretation of nudging that illustrates why the NuPF tends to outperform the BPF when there is some mismatch in the state equation of the SSM. To be specific, we have shown that, given a fixed sequence of observations, the NuPF amounts to a standard PF for a modified dynamical model which empirically leads to a higher model evidence (i.e., a higher likelihood) compared to the original SSM.

The analytical results have been supplemented with a number of computer experiments, both with synthetic and real data. In the latter case, we have tackled the fitting of a stochastic volatility SSM using Bayesian methods for model inference and a time-series dataset consisting of euro-to-US-dollar exchange rates over a period of two years. We have shown how different figures of merit (model evidence, acceptance probabilities or autocorrelation functions) improve when using the NuPF, instead of a standard BPF, in order to implement a nested particle filter <cit.> and a particle Metropolis-Hastings <cit.> algorithm.
 
Since the nudging step is fairly general, it can be used with a wide range of differentiable or non-differentiable likelihoods. Besides, the new operation does not require any modification of the well-defined steps of the PF so it can be plugged into a variety of common particle filtering methods. Therefore, it can be adopted by a practitioner with hardly any additional effort. In particular, gradient nudging steps (for differentiable log-likelihoods) can be implemented using automatic differentiation tools, currently available in many software packages, hence relieving the user from explicitly calculating the gradient of the likelihood. 

Similar to the resampling step, which is routinely employed for numerical stability, we believe the nudging step can be systematically used for improving the performance and robustness of particle filters.

























§ PROOF OF THEOREM <REF>
 

In order to prove Theorem <ref>, we need a preliminary lemma, which can be found, e.g., in <cit.>.
 Let α,β,α̅,β̅∈() be probability measures and f,h ∈ B() be two real bounded functions on  such that (h,α̅) > 0 and (h,β̅)>0. If the identities,

    (f,α) = (fh,α̅)/(h,α̅)  and   (f,β) = (fh,β̅)/(h,β̅)

hold, then we have,

    | (f,α)-(f,β) | ≤
    	1
    /
    	(h,α̅)
    |
    	(fh,α̅) - (fh,β̅)
    | +  f _∞/
    	(h,α̅)
    |
    	(h,α̅) - (h,β̅)
    |.





Now we proceed with the proof of Theorem  <ref>. We follow the same kind of induction argument as in, e.g., <cit.> and <cit.>.

For the base case, i.e. t = 0, we draw x_0^(i), i=1, ..., N, i.i.d. from π_0 and obtain,

    (φ,π_0^N) - (φ,π_0)_p = [|1/N∑_i=1^N (φ(x_0^(i)) - (φ,π_0))|^p]^1/p.

We define S_0^(i) = φ(x_0^(i)) - (φ,π_0) and note that S_0^(i), i = 1,…,N are zero-mean and independent random variables. Using the Marcinkiewicz-Zygmund inequality <cit.>, we arrive at,

    [|1/N∑_i=1^N S_0^(i)|^p]    ≤B_0,p/N^p[(∑_i=1^N | S_0^(i)|^2)^p/2] 
       ≤B_0,p/N^p(N 4 φ_∞^2)^p/2,

where B_0,p is a constant independent of N and the last inequality follows from |S_0^(i)| = | φ(x_0^(i)) - (φ,π_0)| ≤ 2 φ_∞. Therefore, we have proved that Eq. (<ref>) holds for the base case,

    (φ,π_0^N) - (φ,π_0)_p ≤c_0,pφ_∞/√(N),

where c_0,p = 2 B_0,p^1/p is a constant independent of N.

The induction hypothesis is that, at time t-1, 

    (φ,π_t-1^N) - (φ,π_t-1) _p ≤c_t-1,pφ_∞/√(N)

for some constant c_t-1,p < ∞ independent of N. 

We start analyzing the predictive measure ξ_t^N,

    ξ_t^N(x) = 1/N∑_i = 1^N δ_x̅_t^(i)(x),

where x̅_t^(i), i = 1,…,N are the particles sampled from the transition kernels τ_t^x_t-1^(i)(x_t) ≜τ_t( x_t |x_t-1^(i)). Since we have ξ_t = τ_t π_t-1 (see Sec. 1.4), a simple triangle inequality yields,

    (φ,ξ_t^N) - (φ,ξ_t)_p    = (φ,ξ_t^N) - (φ,τ_t π_t-1)_p 
       ≤(φ,ξ_t^N) - (φ,τ_t π_t-1^N)_p 
       + (φ,τ_t π_t-1^N)-(φ,τ_t π_t-1)_p,

where,

    (φ,τ_t π_t-1^N) = 1/N∑_i=1^N (φ,τ_t^x_t-1^(i)).

For the sampling step, we aim at bounding the two terms on the rhs of (<ref>).

For the first term, we introduce the σ-algebra generated by the random variables x_0:t^(i) and x̅_1:t^(i), i = 1,…,N, denoted _t = σ(x_0:t^(i),x̅_1:t^(i),i = 1,…,N). Since π_t-1^N is measurable w.r.t. _t-1, we can write

    [(φ,ξ_t^N) | _t-1] = 1/N∑_i=1^N (φ,τ_t^x_t-1^(i)) = (φ,τ_t π_t-1^N).

Next, we define the random variables S_t^(i) = φ(x̅_t^(i)) - (φ,τ_t π_t-1^N) and note that, conditional on _t-1, S_t^(i), i = 1,…,N are zero-mean and independent. Then, the approximation error of ξ_t^N can be written as,

    [| (φ,ξ^N_t.   .) - (φ,τ_t π_t-1^N)|^p | _t-1]= [ | 1/N∑_i=1^N S_t^(i)|^p |	 _t-1].



Resorting again to the Marcinkiewicz-Zygmund inequality, we can write,

    [ | 1/N∑_i=1^N S_t^(i)|^p |	 _t-1] ≤B_t,p/N^p[ (∑_i=1^N |S_t^(i)|^2  )^p/2|	 _t-1],

where B_t,p< ∞ is a constant independent of N. Moreover, since |S_t^(i)| = |φ(x̅_t^(i)) - (φ,τ_t π_t-1^N)| ≤ 2 φ_∞, we have,

    [ | 1/N∑_i=1^N S_t^(i)|^p |	 _t-1] ≤B_t,p/N^p( N 4 φ_∞^2)^p/2 = B_t,p/N^p/2 2^p φ_∞^p.

If we take unconditional expectations on both sides of the equation above, then we arrive at

    (φ,ξ^N_t   ) - (φ,τ_t π_t-1^N)_p ≤c_1,pφ_∞/√(N),

where c_1,p = 2 B_t,p^1/p < ∞ is a constant independent of N.

To handle the second term in the rhs of (<ref>), we define (φ̅,π_t-1) = (φ,τ_t π_t-1) where φ̅∈ B() and given by,

    φ̅(x) = (φ,τ_t^x).

We also write (φ̅,π_t-1^N) = (φ,τ_tπ_t-1^N). Since φ̅_∞≤φ_∞, the induction hypothesis leads,

    (φ,τ_tπ_t-1^N) - (φ,τ_tπ_t-1)_p    = (φ̅,π_t-1^N) - (φ̅,π_t-1)_p 
       ≤c_t-1,pφ_∞/√(N),

where c_t-1,p is a constant independent of N. Combining (<ref>) and (<ref>) yields,

    (φ, ξ_t^N) - (φ,ξ_t) _p ≤c_1,t,pφ_∞/√(N)

where c_1,t,p = c_t-1,p + c_1,p < ∞ is a constant independent of N.

Next, we have to bound the error between the predictive measure ξ_t^N and the nudged measure ξ̃_t^N. As the sets of samples {x̅_t^(i)}_i=1^N, used to construct ξ_t^N, and {x̃_t^(i)}_i=1^N, used to construct ξ̃_t^N as shown in (<ref>), differ exactly in M particles, namely x̃_t^(j_1), …, x̃_t^(j_M), where { j_1, …, j_M } = ℐ_t, we readily obtain the relationship

    (φ, ξ_t^N) - (φ,ξ̃_t^N) _p    =   1/N∑_i ∈ℐ_t( φ(x̅_t^(i)) - φ(x̃_t^(i)) ) _p 
       ≤   2φ_∞ M/N
       ≤   2 φ_∞/√(N)

where the first inequality holds trivially (since |φ(x) - φ(x')| ≤ 2 φ_∞ for every (x,x') ∈𝖷^2) and the second inequality follows from the assumption M ≤√(N). Combining (<ref>) and (<ref>) we arrive at

    (φ,ξ_t) - (φ,ξ̃_t^N) _p ≤c̃_1,tφ_∞/√(N),

where the constant c̃_1,t,p = 2 + c_1,t,p < ∞ is independent of N. 

Next, we aim at bounding (φ,π_t) - (φ,π̃_t^N)_p using (<ref>). We note that, after the computation of weights, we define the weighted random measure,

    π̃_t^N = ∑_i=1^N w_t^(i)δ_x̃_t^(i)  where   w_t^(i) = g_t(x̃_t^(i))/∑_i=1^N g_t(x̃_t^(i)).

The integrals computed with respect to the weighted measure π̃_t^N takes the form,

    (φ,π̃_t^N) = (φ g_t,ξ̃^N)/(g_t,ξ̃_t^N).

On the other hand, using Bayes theorem, integrals with respect to the optimal filter can also be written in a similar form as,

    (φ,π_t) = (φ g_t,ξ_t)/(g_t,ξ_t).

Using Lemma <ref> together with (<ref>) and (<ref>), we can readily obtain,

    |(φ,π̃_t^N) - (φ,π_t)| ≤   1/(g_t,ξ_t)(
    	φ_∞| (g_t,ξ_t) - (g_t,ξ_t^N)| . 
       . + |(φ g_t,ξ_t) - (φ g_t,ξ_t^N)|),

where (g_t,ξ_t) > 0 by assumption. Using Minkowski's inequality, we can deduce from (<ref>) that

    (φ,π̃_t^N) - (φ,π_t)_p ≤   1/(g_t,ξ_t)(
    	φ_∞ (g_t,ξ_t) - (g_t,ξ_t^N)_p . 
       . + (φ g_t,ξ_t) - (φ g_t,ξ_t^N)_p).

Noting that we have φ g_t_∞≤φ_∞g_t_∞, (<ref>) and (<ref>) together yield,

    (φ,π_t) - (φ, π̃_t^N) _p ≤c_2,t,pφ_∞/√(N),

where 

    c_2,t,p= 2 g_t_∞c̃_1,t,p/(g_t, ξ_t) < ∞

is a finite constant independent of N (the denominator is positive and the numerator is finite as a consequence of Assumption <ref>).

Finally, the analysis of the multinomial resampling step is also standard. We denote the resampled measure as π_t^N. Since the random variables which are used to construct π_t^N are sampled i.i.d from π̃_t^N, the argument for the base case can also be applied here to yield,

    (φ, π̃_t^N) - (φ, π_t^N) _p ≤c_3,t,pφ_∞/√(N),

where c_3,t,p < ∞ is a constant independent of N. Combining bounds (<ref>) and (<ref>) to obtain the inequality (<ref>), with c_t,p = c_2,t,p + c_3,t,p < ∞, concludes the proof.












§ PROOF OF LEMMA <REF>
 



Since x̃_t^(i) = x̅_t^(i) + γ_t ∇_x_t g_t(x̅^(i)_t), we readily obtain the relationships

    | φ(x̃_t^(i)) - φ(x̅_t^(i))|    ≤ Lx̃_t^(i) - x̅_t^(i)_2 
       = L γ_t ‖∇_x_t g(x̅^(i)_t)‖_2 
       ≤γ_t L G_t

where the first inequality follows from the Lipschitz assumption, the identity is due to the implementation of the gradient-nudging step and the second inequality follows from Assumption <ref>. Then we bound the error (φ, ξ_t^N) - (φ, ξ̃_t^N) _p as 

    (φ, ξ_t^N) - (φ,ξ̃_t^N) _p    =   1/N∑_i∈ℐ_t( φ(x̅_t^(i)) - φ(x̃_t^(i)) ) _p 
       ≤   1/N∑_i∈ℐφ(x̅_t^(i)) - φ(x̃_t^(i)) _p 
       ≤   M/Nγ_t L G_t

where the identity is a consequence of the construction of ℐ_t and we apply Minkowski's inequality, (<ref>) and the assumption |ℐ_t|=M to obtain (<ref>). However, we have assumed that sup_1≤ t ≤ Tγ_t M ≤√(N), hence

    (φ, ξ_t^N) - (φ, ξ̃_t^N) _p ≤LG_t/√(N).













§ NUDGING SCHEME THAT INCREASES THE MODEL EVIDENCE
 

Consider the SSM _0={τ_0,τ_t,g_t} where the likelihoods (g_t)_t≥ 1 and the Markov kernels (τ_t)_t≥ 1 satisfy the regularity assumptions below.

 
The functions log g_t(x), t =1, 2, ..., are differentiable and the gradients ∇log g_t(x) are Lipschitz with constant L_t^g < ∞. To be specific,

    ∇log g_t(x) - ∇log g_t (x') _2 ≤𝖫_t^g  x - x'_2.



 
The Markov kernels τ_t( x|x') are absolutely continuous with respect to the Lebesgue measure, hence there are conditional pdf's m_t(x|x') such that τ_t( x|x') = m_t(x|x')  x for any x'∈. Moreover, the log-pdf's log m_t(x|x') are uniformly Lipschitz in x', i.e., there are non-negative bounded functions L_t(x) such that

    | log m_t(x|x') - log m_t(x|x”) | ≤ L_t(x)  x' - x”_2

and
𝖫^τ_t = sup_x∈ L_t(x) < ∞.




For any subset A ⊆, let us introduce the indicator function 

    1_A(x) := {[ 1,  ,; 0,   ;    ].

We construct the nudging operator α_t^y_t : ↦ of the form

    α_t^y_t(x) := (
    	x + γ_t ∇log g_t(x)
    )  1_S_y_t(x) + x  1_S_y_t(x),

where S_y_t := { 
	x ∈: ∇log g_t(x_t) ≥ 2 𝖫_t^τ} (recall that g_t(x)=g_t(y_t|x)), S_y_t = \ S_y_t is the complement of the set S_y_t and γ_t>0 is small enough to guarantee that g_t(x+γ_t∇log g_t(x)) ≥ g_t(x). Intuitively, this nudging scheme only takes a gradient step when the slope of the likelihood g_t is sufficient to insure an improvement of the likelihood with a small move of the state x.

Assume, for simplicity, that we apply the nudging operator (<ref>) to every particle at every time step. The recursive step of the resulting NuPF can be outlined as follows:


  * For i=1, ..., N,
	
	
  * draw x̅_t^(i)∼τ_t( x_t|x_t-1^(i)),
	
  * nudge every particle, i.e., x̃_t^(i) = α_t^y_t(x̅_t^(i)),
	
  * and compute weights w_t^(i)∝ g_t(x̃_t^(i)).
	

  * Resample to obtain { x_t^(i)}_i=1, ..., N.

The asymptotic convergence of this algorithm can be insured whenever the step sizes γ_t are selected small enough to guarantee that sup_x  x - α_t^y_t(x)≤1/√(N) (this is a consequence of Corollary 1 in <cit.>). The implicit model for this NuPF is _1={τ_0,τ̃_t,g_t} where the transition kernel is

    τ̃_t^y_t( x_t|x_t-1) = ∫δ_α_t^y_t(x̅)( x_t) τ_t(x̅|x_t-1).

Note that for any integrable function f:↦ we have

    (f,τ̃_t^y_t)(x_t-1)    =   ∫∫ f(x_t) δ_α_t^y_t(x̅)( x_t) τ_t(x̅ | x_t-1) 
       =   ∫ (f ∘α_t^y_t)(x̅) τ_t(x̅ | x_t-1) 
       =    (f ∘α_t^y_t, τ_t)(x_t-1),

where (f ∘α_t^y_t)(x)=f(α_t^y_t(x)) is the composition of f and α_t^y_t. In particular, we note that (g_t,τ̃_t^y_t)(x)= (g_t∘α_t^y_t,τ_t)(x) and ( 1_,τ̃_t^y_t)(x) = ∫ 1_(α_t^y_t(x)) τ_t(x̅ | x) = 1 because α_t^y_t is ↦ and, therefore, 1_∘α_t^g_t = 1.

Let β_t(x) := x + γ_t ∇log g_t(x).  Assumptions <ref> and <ref> entail the following result.









 
If Assumptions <ref> and <ref> hold and the inequalities

    γ_t ≤1/𝖫_t^g    ∇log g_t(x_t) _2 ≥ 2 𝖫_t^τ,
 
are satisfied, then 

    (g_t ∘β_t)(x_t)
    /
    	g_t(x_t)
    ≥sup_x_t+1∈
    	m_t+1( x_t+1 | x_t)
    /
    	m_t+1(x_t+1 | β_t(x_t))
    .





Proof. Assumption <ref> implies that, for any pair x,x'∈,

    log g_t(x) ≥log g_t(x') + ⟨∇log g_t(x'), x - x'⟩ - 𝖫_t^g/2x - x'^2_2

see, e.g., <cit.> for a proof. We can readily use (<ref>)  to obtain a lower bound for log g_t(β_t(x_t)). Indeed,

    log g_t(β_t(x_t))    ≥   log g_t(x_t) + ⟨∇log g_t(x_t), γ_t ∇log g_t(x_t)⟩
           - 𝖫_t^gγ_t^2/2∇log g_t(x_t)^2_2 
    
       =   log g_t(x_t) + (γ_t -  𝖫_t^gγ_t^2/2) ∇log g_t(x_t)^2_2 
       ≥   log g_t(x_t) + γ_t/2∇log g_t(x_t)^2_2

where the last inequality follows from the assumption γ_t ≤1/𝖫_t^g. In turn, (<ref>) implies

    g_t(β_t(x_t))/g_t(x_t)≥exp{γ_t/2∇log g_t(x_t)^2_2}.


We now turn to the problem of upper bounding the ratio of transition pdf's. From Assumption <ref> and the definition of β_t(x_t) we readily obtain that

    log m_t+1(x_t+1 | x_t) - log m_t+1(x_t+1 | β_t(x_t)) ≤𝖫_t^τγ_t ∇log g_t(x_t)_2

holds for any x_t+1∈. Taking exponentials on both sides of (<ref>) we arrive at

    sup_x_t+1∈m_t+1(x_t+1 | x_t)/m_t+1(x_t+1 | β_t(x_t))≤exp{𝖫_t^τγ_t ∇log g_t(x_t)_2 }.


If ∇log g_t(x_t) _2=0 then expressions (<ref>) and (<ref>) together readily yield the desired relationship (<ref>).

If ∇log g_t(x_t) _2 > 0, then the assumption ∇log g_t(x_t) ≥ 2𝖫_t^τ implies that

    γ_t/2∇log g_t(x_t)_2^2 ≥γ_t 𝖫_t^τ∇log g_t(x_t)_2

which, together with (<ref>) and (<ref>), again yield the desired inequality (<ref>). 

Finally, we prove that the evidence in favour of _1 is greater than the evidence in favour of  _0.



Let the nudging scheme be defined as in (<ref>). If Assumptions <ref> and <ref> hold and the inequality

    γ_t ≤1/𝖫_t^g
 
is satisfied for t=1, …,  T<∞, then (y_1:T|_1) ≥(y_1:T|_0).



Proof. From the definition of τ̃_t in (<ref>) and the ensuing relationship (<ref>), the evidence of model _1 can be readily written down as

    (y_1:T|_1)    =   ∫⋯∫ g_T(α_t^y_T(x_T)) ×
          ×∏_t=1^T-1 m_t+1(x_t+1|α_t^y_t(x_t)) g_t(α_t^y_t(x_t)) 
          × m_1(x_1|x_0)m_0(x_0)  x_0 ⋯ x_T.

It is apparent that g_T(α_T^y_T(x_T)) ≥ g_T(x_T) for every x_T. Moreover, for any x_t ∈, if x_t ∈ S_y_t then 

    ∇log g_t(x_t) ≥ 2 𝖫_t^τ    α_t^y_t(x_t) = β_t(x_t),
 
hence we can apply Lemma <ref>, which yields 

    m_t+1(x_t+1|α_t^y_t(x_t)) g_t(α_t^y_t(x_t)) ≥ m_t+1(x_t+1|x_t) g_t(x_t)

for every x_t+1. Alternatively, if x_t ∉ S_y_t, then α_t^y_t(x_t)=x_t and, trivially,

    m_t+1(x_t+1|α_t^y_t(x_t)) g_t(α_t^y_t(x_t)) = m_t+1(x_t+1|x_t) g_t(x_t).

Therefore,

    (y_1:T|_1)    ≥   ∫⋯∫ g_T(x_T) ∏_t=1^T-1 m_t+1(x_t+1|x_t) g_t(x_t) 
          × m_1(x_1|x_0) m_0(x_0)  x_0 ⋯ x_T 
       =   (y_1:T|_0).








spbasic













