Monte Carlo fixed-lag smoothing in state-space models

This paper presents an algorithm for Monte Carlo fixed-lag smoothing in state-space models defined by a diffusion process observed through noisy discrete-time measurements. Based on a particles approximation of the filtering and smoothing distributions, the method relies on a simulation technique of conditioned diffusions. The proposed sequential smoother can be applied to general non linear and multidimensional models, like the ones used in environmental applications. The smoothing of a turbulent flow in a high-dimensional context is given as a practical example.


Introduction
The framework of this paper concerns state-space models described by general diffusions of the 10 form: which are partially observed through noisy measurements at discrete times. Such models can describe many dynamical phenomena in environmental sciences, physics, but also in finance or engineering applications. The main motivation of this work concerns environmental applications, where 15 non linearity and high-dimensionality arise. Indeed, environmental models and data describe non linear phenomena over large domains, with high spatial resolution. The continuous dynamical model (1) is defined from a priori physical laws, while observations are supplied by sensors (satellite data for instance) and can appear with very low time frequency. As an example, in the application presented in the last part of this paper, the dimension of the state and observations is of the order of 20 many thousands, and the model is described by the non linear Navier-Stokes equation. Filtering 1 and smoothing in such state-space models aim at coupling model and observations, which is called data assimilation. The goal of the filtering is to estimate the system state distribution knowing past and present observations. This allows for instance to give proper initial conditions to forecast the future state of a system characterizing atmospheric or oceanographic flows. On the other hand, the 25 smoothing aims at estimating the state distribution using past and future observations, and this retrospective state estimation allows to analyze a spatio-temporal phenomenon over a given time period, for climatology studies for instance. Applications of data assimilation are numerous and the interest is growing in environmental sciences with the increase of available data. However, it is still a challenge to develop filtering and smoothing methods that can be used within a general non linear and 30 high-dimensional context.
Monte Carlo sequential methods, contrary to standard Kalman filters, are able to deal with the filtering problem in non linear state-space models. The particle filtering (Del Moral et al., 2001;Doucet et al., 2000) solves the whole filtering equations through Monte Carlo approximations of the state distribution. On the other hand, ensemble Kalman methods (Evensen, 2003) take into account 35 in some way the non linearities in the system, but are based on a Gaussian assumption. For highdimensional systems, ensemble Kalman methods are preferred in practice to particle filters (Stroud et al., 2010;Van Leeuwen, 2009) since they reach better performance for limited number of particles.
In order to keep this advantage while alleviating the Gaussian assumption, both methods are combined in Papadakis et al. (2010), leading to a particle filter that can be applied to high-dimensional 40 systems. We will use this technique for the filtering step in the high-dimensional application presented in Section 5.
The aim of this paper is to propose a new smoothing method. Within the particle filter framework, the smoothing can be computed backward, reweighting past particles using present observations (Briers et al., 2010;Godsill et al., 2004). There is however one main difficulty for continuous 45 models of type (1). As a matter of fact, it is necessary to know the transition density of the process between observation times, which is not available for general diffusions. This transition density can be approximated through Monte Carlo simulations, as proposed by Durham and Gallant (2002) to solve inference problems for diffusion processes. However, these approximations are based on Brownian bridge (or modified versions of it) simulations, that do not take into account the drift part 50 of the model. For non linear and high-dimensional models with a drift term that dominates, such approximations will be inefficient. It is also possible to obtain an unbiased estimate of the transition density (see Beskos et al. (2006)), but this approach is not adapted to a multi-dimensional context. As a matter of fact, the use of this technique in a multivariate setting imposes constraints on the diffusion drift (in particular the drift function has to be of gradient type). In parallel, within the framework of 55 ensemble Kalman methods, Evensen and van Leeuwen (2000) have proposed to estimate backward the smoothing distribution in a recursive way, based on existing filtering trajectories; Stroud et al. (2010) presented and applied an ensemble Kalman smoothing method, relying on a linearization of 2 the system dynamics.
All previously mentioned smoothing methods require to perform specific assumptions or simpli-60 fications in order to deal with general non linear models of type (1) in a high-dimensional context.
To the best of our knowledge, it remains a challenging problem to develop smoothing methods that can be used in this general setting. In this paper, we deal with this issue sequentially each time a new observation is available, by smoothing the hidden state from this new observation time up to the previous one. This approach, called fixed-lag smoothing, constitutes then a partial answer to the global smoothing problem that would take into account all available observations. Nevertheless, it is reasonable to assume that the distribution of the hidden state depends on future observations through the next observation only, as soon as the time step between measurements is long (which is typically the case in the environmental applications that motivate this work). Under this assumption, a new observation will impact the distribution of the hidden process up to the previous observation 70 only. This point of view justifies the use of a fixed-lag smoothing in our setting as a reasonable approximation of the global smoothing problem.
Such a fixed-lag smoothing may be directly obtained from the particle filtering result, reweighting past trajectories. However, this smoothing will fail in two cases: when the number of particles is too small compared to the size of the system, or when existing trajectories do not correspond to 75 plausible trajectories of the dynamical model. Unfortunately, these two situations have to be faced when smoothing in a high-dimensional state-space model. Firstly, the number of particles has to be reduced for computational reasons. Secondly, particle filters that have been proposed in this highdimensional context require to correct trajectories towards the observations (Papadakis et al., 2010;Van Leeuwen and Ades, 2013). This implies that filtering states are consistent at observation times, 80 but that filtering trajectories may not be plausible realizations of the underlying physical model.
In that case, a smoothing based on existing trajectories will fail. Note that these remarks are not only valid for the fixed-lag smoothing, but also for previously mentioned global techniques relying on existing trajectories. In particular, a genealogical smoothing based on ancestral particle lines (Del Moral, 2004) will be deficient in a high-dimensional setting since many trajectories will share 85 only a few ancestral lines.
In contrast, our method does not rely on existing particles only. It is built on a conditional simulation technique of diffusions proposed by Delyon and Hu (2006) that provides new state trajectories at hidden times between observations. This simulation technique is adapted to a multivariate context where the drift dominates, contrary to techniques based on Brownian bridge sampling (Durham and 90 Gallant, 2002). Moreover, it does not require constraining assumptions for multivariate models, contrary to other techniques based on exact simulation of diffusions (Beskos and Roberts, 2005;Beskos et al., 2006). The proposed smoothing method can then be applied to high-dimensional systems.
Finally, it does not require model linearization nor Gaussian hypotheses, and so is able to deal with general non linear models.

95
3 The remaining of the paper is organized as follows. Section 2 briefly describes sequential Monte Carlo filtering methods in state-space models, and presents the fixed-lag smoothing problem. Section 3 presents the conditional simulation technique of diffusions of Delyon and Hu (2006), and details the construction of the proposed Monte Carlo estimate of smoothing distributions. The method is then experimented on a one-dimensional example in Section 4. Finally, the method is applied in 100 section 5 to a practical non linear and high-dimensional case, similar to the problems that have to be faced in environmental applications. A discussion is given in Section 6.

Monte Carlo filtering and smoothing in state-space models
In this section we recall briefly the particle filtering and smoothing methods for models of type (1), where the hidden state vector x ∈ R n is observed through the observation vector y ∈ R m at discrete 105 times {t 1 , t 2 , . . .}: The drift function f and observation operator g can be non linear. The dynamical model uncertainty is described by a n-dimensional Brownian motion with covariance Σ = σ(x(t))σ(x(t)) T . The functions f , g and σ are assumed to be known, as well as the law of the observation noise γ t k .

110
In particular, we present the standard particle filter and the Weighted Ensemble Kalman filter, that can be used to face the filtering problem in high-dimensional systems.

Particle-based filtering methods
Filtering aims at estimating recursively the distribution p(x t1:t k |y t1:t k ) (and in particular its marginal distribution p(x t k |y t1:t k )) at each observation time t k . This filtering problem can be solved with a

115
Monte Carlo sequential approach, called particle filtering (Del Moral et al., 2001;Doucet et al., 2000). The method relies on a Monte Carlo approximation of the filtering distribution over a set of whose marginal distribution at time t k writes: Particle filters rely on a sequential importance sampling scheme that recursively samples particles, and updates their weights at observation times. The weights corresponds to the ratio between the target distribution and the importance sampling distribution π(x t k |x t0:t k−1 , y t1:t k ). They are recur-4 sively computed as follows: . (5) In practice, a resampling procedure is added in order to avoid degeneracy. This procedure duplicates trajectories with large weights and remove small weighted trajectories.

Standard particle filter
When the proposal distribution π is set to the prior (i.e. π(x t k |x t0:t k−1 , y t1:t k ) = p(x t k |x t k−1 ), the 130 weights updating rule (5) simplifies to the computation of the data likelihood p(y t k |x . This particular instance of the particle filter is called the Bootstrap filter or sequential importance resampling (SIR) filter Gordon et al. (1993). Due to its simplicity it is the most commonly used particle filter. It is however a very inefficient distribution for high dimensional space as it does not take into account the current observation and depends only weakly on the past data through the filtering distribution 135 estimated at the previous instant. This distribution requires a great number of particles to explore meaningful areas of the state space.

Weighted Ensemble Kalman filter (WEnKF)
One way to efficiently incorporate observation within the proposal distribution consists to rely on the ensemble Kalman filtering mechanism to define this distribution. This is the idea proposed in 140 the WEnKF technique (Papadakis et al., 2010). In the WEnKF approach the importance sampling is taken as a Gaussian approximation of p(x t k |x t k−1 , y t k ). This approach is close to the technique proposed in Van Leeuwen (2010). A variation of a similar technique based on a deterministic squareroot formulation is also described in Beyou et al. (2013). In order to make the estimation of the filtering distribution exact (up to the sampling), each member of the ensemble must be weighted at 145 each observation instant t k with appropriate weights w (i) t k , defined from (5). Therefore, the Weighted ensemble Kalman filter (WEnKF) procedure can be simply summarized by Algorithm 1.

Algorithm 1 The WEnKF algorithm
For each t k = t 1 , t 2 , . . .: Note that particle-based filtering techniques update the filtering distribution at observation times only. However, after the estimatep(x t k |y t1:t k ) has been updated at observation time t k , the filtering distribution can be predicted in order to have a continuous estimation ofp(x t |y t1:t k ) for all t ∈ 150 ]t k , t k+1 [ until the next observation time: where, for all i = 1, . . . , N , the state x

155
Contrary to the filtering approach that uses past and present observations, the smoothing in statespace models aims at estimating p(x t |y t1:tend ) for all t ∈ [t 1 , t end ], using all past and future observations over a given time period. As raised in the introduction, existing smoothing methods do not apply directly to a general non linear model of type (1) in a high-dimensional context, since assumptions have to be made that may not be realistic. Instead of solving the global smoothing, we will 160 concentrate in the rest of the paper on a fixed-lag smoothing, which constitutes a partial answer to the global smoothing problem.
The objective of the fixed-lag smoothing will be to replace the predictive distribution (6) by its smoothed version p(x t |y t1:t k+1 ) ∀t ∈]t k , t k+1 ], sequentially each time a new observation y t k+1 arrives. This will allow to reduce the temporal discontinuities inherent to the filtering technique, that 165 successively predicts the distribution of the state between observations, and updates this distribution at observation times.
To achieve this, by construction of the particle filter that weights entire trajectories (see equation

6
(3)), it is known (see for instance Doucet et al. (2000)) that the fixed-lag smoothing distribution p(x t |y t1:t k+1 ) can be directly obtained from the marginal at time t ofp(x t1:t k+1 |y t1:t k+1 ). The 170 empirical smoothing distribution is then given by: However, this approximation is simply a reweighting of past existing particle trajectories, and relies on the support of the filtering distribution at time t k . If the number of particles is too small with respect to the state dimension, the support may be greatly reduced by the correction step (assigning 175 small weights to all particles except a few), leading in practice to a bad estimation of p(x t |y t1:t k+1 ).
Moreover, if particle trajectories have been forced towards observations during the filtering step (like in the WEnKF procedure), a smoothing based on those particles will fail because it will not be able to correct discontinuities. Consequently, since we are interested in smoothing techniques that are efficient in a high-dimensional context, this direct smoothing technique can not be used in its basic 180 form and has to be improved.
In the following, we propose to use a conditional simulation technique of diffusions that will enable the sampling of new smoothed trajectories between times t k and t k+1 . The approximation of the smoothing distribution (7) at each hidden time will then be improved. The conditional simulation technique is presented in the next section, before the resulting smoothing procedure we propose.

Fixed-lag smoothing with conditional simulation
The smoothing method we propose is based on a conditional simulation technique that is presented in section 3.1. We develop then in section 3.2 how this technique can be used to improve the estimation of the smoothing distribution (7).

Conditional simulation 190
Conditional simulation of a diffusion aims at sampling trajectories from a given process: between two times t = 0 and t = T , with the constraints x(0) = u and x(T ) = v. This simulation problem is treated by Delyon and Hu (2006), where the authors show how to obtain the law of the constrained process from a Girsanov theorem. In practice, the proposed algorithms consist in 195 simulating trajectories according to another diffusion process, which is built to respect the constraints and is easy to simulate from. The conditional distribution of the constrained process (8) is absolutely continuous with respect to the distribution of the auxiliary process, with explicitly given density. For instance, in the case where the drift is bounded (a similar algorithm is proposed in Delyon and Hu (2006) for the unbounded case) and for σ invertible, the algorithm is based on the simulation of 200 7 trajectories from the following process: with initial conditionx(0) = u. Note that Lemma 4 in Delyon and Hu (2006) deals with the existence of a unique solution for this equation. This process is a simple modification of (8), where a deterministic part is added to the drift. It is then easy to simulate unconditional trajectories from this 205 process, and all simulated trajectories will satisfyx(T ) = v by construction. For simplicity we will assume in the following that σ is independent of x(t) (note however that this is not an assumption in Delyon and Hu (2006)). The law of the conditioned process is given by: for all measurable function h, where: is the density coming from Girsanov theorem (see Delyon and Hu (2006) Let us note that the presence of the drift part of model (8) in the auxiliary process (9) is crucial to make the simulation efficient. The same process had initially been proposed by Clark (1990) to solve the conditional simulation problem. On the other hand, standard Brownian bridges that could 215 be used as auxiliary processes (Durham and Gallant, 2002) lead in practice to poor approximations of the original constrained diffusion in our high-dimensional setting, since Brownian bridge trajectories are too far away from trajectories of (8).
In the following, the conditional marginal of interest p(x t |x(0) = u, x(T ) = v) will then be approximated as follows: where the M trajectories {x

Proposed fixed-lag smoothing method
We show in the following how the conditional simulation technique can be used to improve the estimation of the local smoothing distribution p(x t |y t1:t k+1 ) for all t ∈]t k , t k+1 ]. 225 We first note that this distribution can be decomposed as: Then, from the state-space model properties, we obtain: 230 Moreover, from the particle filter Monte Carlo approximation described by (3), the joint law p(x t k , x t k+1 |y t1:t k+1 ) can be replaced by: where the w (i) t k+1 are the particle filter importance weights. Plugging (15) into (14) leads then to the following approximation for the fixed-lag smoothing distri-235 bution: The conditional distribution p(x t |x can be estimated using (12) for each pair of initial and end points x The estimation of the smoothing distribution of interest writes finally: The algorithm we propose to compute the fixed-lag smoothing distribution on a given time interval

One-dimensional simulation study
In this section, the smoothing method is first experimented on a one-dimensional state space model.
Since the proposed approach relies on a preliminary particle filtering step, filtering results are first 250 presented in Section 4.2 (either considering a standard particle filter or the WEnKF). The results obtained with the standard fixed-lag smoothing method are then shown in Section 4.3. Finally, Section 4.4 presents the smoothing results obtained with the proposed technique.

State space model
The one-dimensional state space model of interest is a sine diffusion, partially observed with noise 255 (used as an illustration by Fearnhead et al. (2008) for a particle filtering method) : where σ 2 x = 0.5 and γ t k ∼ N (0, σ y ) with σ 2 y = 0.01. One trajectory of the process is first simulated from (19) with an Euler-type discretization scheme of time step ∆t = 0.005. This trajectory will 260 constitute the hidden process, observed through y t k generated according to (20) at every time step t k , with t k − t k−1 = 20∆t. The trajectory is plotted on Figure 1, together with the corresponding discrete observations at times t k .

Particle filtering results
The filtering results are presented for the standard particle filter (denoted PF in the following) and the Weighted Ensemble Kalman filter (WEnKF). Two situations are shown, with reduced (N = 20) and high number (N = 10000) of particles. The case with a high number of particles is shown as the reference for comparison purpose, note however that this ideal situation is not reachable in a 270 high-dimensional context, since the number of particles has to be reduced for evident computational cost reasons.
The results for the two configurations are presented on Figure 2, where the dotted lines represents the filtering mean estimates. The filtering distribution p(x t k |y t1:t k ) is estimated at each observation time t k using (4), and predicted between observation times from (6). The mean is then estimated 275 from weighted particles as Figure 2 (a)-(b) show that the standard particle filter results diverge from the reference solution between observation times, for low or high number of particles. As a matter of fact, when no observation is available, the state distribution is predicted from the dynamics only, so that particles trajectories are not guided towards the next observation. At observation times t k , high weights are given to particles that are close to the 280 observation, so that the estimated mean suddenly gets closer to the solution. These discontinuities between measurement times can also be observed on the WEnKF results (Figure 2 (c)-(d)), because particle trajectories are brutally corrected with the EnKF step at observation times. A smoothing will aim at reducing these temporal discontinuities while providing dynamically consistent solutions.

Standard fixed-lag smoothing results
From the particle filtering results, we present now the results obtained with the direct particles smoothing procedure described in section 2.2. This procedure relies on existing trajectories. The smoothing distributionp(x t |y t1:t k+1 ) is computed backward for all t ∈]t k , t k+1 ] using expression 290 (7), each time a new observation y t k+1 becomes available. The smoothing mean is computed as t for all t ∈]t k , t k+1 ], and the standard deviation is computed in the same way from the weighted particles.
It can be observed on Figure 3(a) that the smoothing based on the standard particle filter is not efficient when the number of particles N is small: Only a few particles are close to the observation 295 at time t k and have nonzero weights, implying that the smoothing distribution is poorly estimated (see for instance between observation times t = 100 and t = 120 where the smoothing distribution is artificially peaked but far from the hidden trajectory). The smoothing result obtained from the reference configuration N = 10000 is plotted on Figure 3(b). In that situation, since many trajectories have high weights at observation times, the estimation of backward smoothing distributions is

Proposed smoothing results
In this section, we show how the proposed method can improve the estimation of backward smooth- Our smoothing is first applied using the filtering output of the standard particle filter with N = 20 315 particles. Figure 4(a) shows the result obtained with a sampling of M = 50 conditional trajectories between each pair {x (18), so the smoothing mean is computed as for all t ∈]t k , t k+1 ], and similarly for the standard deviation. This result highlights the fact that since the proposed method creates new trajectories, it is able to correct the deficiencies of the standard 320 smoothing approach presented on Figure 3(a) when the initial number of filtering particles is too small. On Figure 4(b), the same experiment is presented using M = 500 conditional trajectories. In that case, the result is very similar to the reference particles smoothing result presented on Figure   3(b), obtained from a particle filter with N = 10000.

13
In parallel, the proposed smoothing has been tested using the output of the WEnKF filtering  14

Application to a high-dimensional assimilation problem
This section aims at illustrating the applicability of our method to a high-dimensional and non linear scenario, without extensive study at this stage. The method is applied to a turbulence assimilation problem, where the model of interest is of type (1). The goal is to recover temporal estimates of velocity/vorticity over a given spatial domain of size n = 64 * 64, from a sequence of noisy obser-340 vations and a continuous a priori dynamical model based on a stochastic version of Navier-Stokes equation. Within an environmental framework, a direct application would be the estimation of wind fields or sea surface currents from satellite data.

State space model
Let ξ(x) denote the scalar vorticity at point x = (x, y) T , associated to the 2D velocity w(x) =

345
(w x (x), w y (x)) T through ξ(x) = ∂wy ∂x − ∂wx ∂y . Let ξ ∈ R n be the state vector describing the vorticity over a n = 64 * 64 square domain, and w ∈ R 2n the associated velocity field over the domain. We will focus on incompressible flows such that the divergence of the velocity field is null. A stochastic version of Navier-Stokes equation in its velocity-vorticity form can then be written as:

350
where ℜ denotes the flow Reynolds number (Re = 3000). The uncertainty is modeled by a Brownian motion of size n, with covariance Σ = σσ T , where σ ∈ R n . A velocity field example, generated from the model (21), is shown on Figure 5(a), together with the corresponding vorticity map (b).
We assume the hidden vorticity vector ξ is observed through noisy measurements y t k at discrete times t k , where t k − t k−1 = 100∆t, and ∆t = 0.1 is the time step used to discretize (21). In our 355 experimental setup, measurements correspond to PIV (Particle Image Velocimetry) image sequences used in fluid mechanics applications. Note however that other kind of data can be used similarly within this state space model, like meteorological or oceanographic data for instance. The state and observation are related in our case through y t k = g(ξ t k ) + γ t k , where g is a non linear function linking the vorticity to the image data, and γ t k is a Gaussian noise, uncorrelated in time.

Implementation details
We recall that the smoothing relies first on a particle filter step. Due to the high dimensionality of the state vector, the use of a standard particle filter is not adapted to solve the filtering problem, as dis- t1:t k , w t k } i=1:N , as required by the algorithm proposed in section 3.
The particle filter step requires simulations from the dynamical model (21), and the conditional simulation step requires to sample trajectories from its constrained version, which consists in a similar problem with modified drift (see process (9)). The model is discretized in time with time 375 step ∆t = 0.1; more information about the discretization scheme may be obtained in Papadakis et al. (2010). The random perturbations are assumed to be realizations of Gaussian random fields that are correlated in space with Gaussian covariance function Σ(x i , x j ) = η exp(− ||xi−xj || 2 λ ), where η = 0.01 and λ = 13. In practice, the simulation of these perturbations is performed in Fourier space, with the method described in Evensen (2003). The inverse of the covariance matrix Σ −1 is finally computed as: which only requires the inversion of a diagonal. Note that more efficient procedures could be implemented in our case (homogeneous Gaussian covariance) since the covariance function is separable in x and y directions. This means that the covariance matrix Σ can be written as the Kronecker product of smaller matrices and more easily inverted (Sun et al., 2012). However, the SVD inversion can be applied to any covariance structure, in particular it could deal with a non homogeneous covariance matrix.

Results
In this section, we illustrate the capability of the proposed method to reduce the temporal discontinuities inherently introduced by the filtering in continuous-discrete state-space models.
The smoothing result relies on the output of the WEnKF filtering step, computed with N = 500 particles. Compared to the size of the system, the number of particles is theoretically too small for 400 the filter to be truly efficient. In practice, many filtering trajectories have close to zero weights at observation times. Histograms of filtering weights are given as illustration on Figure 6(a)-(b) at two times t = 400 and t = 500. Note however that the filter is not degenerate and is able to provide results that get close to he hidden vorticity at measurement times. This can be observed on Figure 7, where the mean square error is plotted with full line, averaged at each time over the image domain 405 of size n = 64 * 64. Since the ground truth vorticity sequence is known in our experimental setup, the mean square error is computed between the hidden vorticity and the estimated filtering mean, [. The correction steps lead to successive error decreases at observation times.
The proposed smoothing method has been applied with M = 200. Note that we take benefit 410 from the fact that many filtering have close to zero weights. Indeed, the smoothing method relies in practice on a reduced numberÑ M of sampled conditional trajectories (withÑ << N ), which makes the problem computationally tractable. On this experiment, we have retained around 5% of initial filtering trajectories. The smoothing distributionp(ξ t |y t1:t k+1 ) is computed for all t ∈ ]t k , t k+1 ] from (18), and its mean is computed as . Histograms 415 of conditional simulation weights α(ξ (i)(j) ) are given as illustration on Figure 6(c)-(d) for a given particle (i) at two times t = 400 and t = 500.
The mean square error is computed between the true vorticity and the estimated smoothing mean, and plotted on Figure 7 with dotted line. As expected, the smoothing method reduces the error at hidden times between observations. In addition, we present below a qualitative evaluation of the smoothing result for the same exper-425 iment, over a specific time interval.
The WEnKF result is first presented on Figure 8  t for t = 500. The temporal discontinuity between estimations can be observed when reaching observation time t = 500: the vorticity map is suddenly modified in order to fit to 430 the observations, introducing inconsistencies in the vorticity temporal trajectories. Note that the application of the standard particles smoothing (described in section 2.2) will fail here, and not only because the number of particles is too small. As a matter of fact, we recall that the filtering trajectories have been computed from the method presented in Papadakis et al. (2010), which uses the ensemble Kalman filter step as importance distribution in the particle filter algorithm. The ensemble 435 Kalman filter consists of a prediction step from the dynamical model (21), and a correction step which shifts particles towards the observation. Because of this correction step, the sampled filtering trajectories between two observation times do not correspond to trajectories of the dynamical model. This implies that from such a particle filter, the standard smoothing based on existing trajectories will not be able to reduce the temporal discontinuities observed on Figure 8. This can be 440 observed on Figure 9, where smoothed vorticity maps are computed as  the a priori dynamical model. In order to sample the smoothed trajectories, the method relies on the model and on filtering marginals at observation times, but not on filtering trajectories at hidden times. It is then able to smooth the discontinuities inherent to the particle filtering technique we have used, contrary to the standard smoothing presented on Figure 9.

Conclusion and discussion
In this paper we introduced a smoothing algorithm based on a conditional simulation technique of diffusions. The proposed smoothing is formulated as fixed-lag, in the sense that it is performed 460 sequentially each time a new observation appears, in order to correct the state at hidden times up to the previous observation. Note that a decomposition similar to equations (13) to (18) can be written from an integration up to a previous time t k−h , with h > 1. This implies that the smoother can be formulated with a larger fixed-lag, in order to correct the state backward not only up to the previous observation, but up to further measurement times. Yet, due to the successive resampling steps that 465 have been performed in the filtering steps before time t k , there are in practice only a few distinct filtering trajectories at times t k−h if h is large. Consequently, the estimation of the joint law in (15) will not be reliable anymore for a too large value of h.
We have shown the practical applicability of the method to a high-dimensional problem. Nevertheless, the algorithm remains costly since a second Monte Carlo step is added to the Monte Carlo 470 nature of particle filter algorithms. Yet, from an algorithmic point of view, the sequential nature of the proposed technique allows the smoothing to be implemented with a similar structure as filtering methods (sequential sampling and weighting of model trajectories). It is then easy to couple this smoothing to an operational filtering system and benefit from parallelization strategies for instance.
Since the proposed smoothing uses the filtering result as input, it relies on the success of the 475 21 underlying particle filter. For high-dimensional systems, a standard particle filter is not adapted and it is necessary to use filtering techniques that guide particles towards observations. In this paper, we use the WEnKF algorithm. In practice, any efficient particle filtering technique with such a guiding can be used within our framework. Note however that the construction of such techniques remains an open area of research. 480 We plan to work on the application of the smoothing method to a real high-dimensional case (for the estimation of sea surface currents from satellite image sequences). However, such a work will imply numerous difficulties which are not related to the smoothing technique but to the definition of the state-space model: definition of a suitable physical model, good characterization of state noise structure and model parameters. Therefore, this will be part of a future work.