Ensemble transform algorithms for nonlinear smoothing problems

TL;DR

This paper explores ensemble transform algorithms for nonlinear smoothing, introducing a new ensemble transform particle smoother that is consistent and flexible, with practical techniques tested on complex problems.

Contribution

It introduces the ensemble transform particle smoother, a novel method that is consistent in the particle limit and adaptable for complex nonlinear smoothing tasks.

Findings

The ensemble transform particle smoother is consistent in the particle limit.

Adaptive spread corrections improve smoothing accuracy.

Localization techniques facilitate application to complex problems.

Abstract

Several numerical tools designed to overcome the challenges of smoothing in a nonlinear and non-Gaussian setting are investigated for a class of particle smoothers. The considered family of smoothers is induced by the class of linear ensemble transform filters which contains classical filters such as the stochastic ensemble Kalman filter, the ensemble square root filter and the recently introduced nonlinear ensemble transform filter. Further the ensemble transform particle smoother is introduced and particularly highlighted as it is consistent in the particle limit and does not require assumptions with respect to the family of the posterior distribution. The linear update pattern of the considered class of linear ensemble transform smoothers allows one to implement important supplementary techniques such as adaptive spread corrections, hybrid formulations, and localization in order to facilitate their application to complex estimation problems. These additional features are derived and numerically investigated for a sequence of increasingly challenging test problems.

Equations

$${\mathbf{x}}_k=f({\mathbf{x}}_{k-1})+{ϵ}_k,{ϵ}_k\sim \mathcal{N}(0,\mathbf{Q})$$

$${\mathbf{y}}_k=h({\mathbf{x}}_k)+{\nu }_k,{\nu }_k\sim \mathcal{N}(0,\mathbf{R})$$

$${\mathbf{y}}_{1:k}=[{\mathbf{y}}_1^{\mathrm{T}},{\mathbf{y}}_2^{\mathrm{T}},...,{\mathbf{y}}_k^{\mathrm{T}}{]}^{\mathrm{T}}\in {\mathbb{R}}^{k{N}_y\times 1}$$

$${\mathbf{x}}_k=f({\mathbf{x}}_{k-n},\dots ,{\mathbf{x}}_{k-1})+{ϵ}_k,{ϵ}_k\sim \mathcal{N}(0,\mathbf{Q}),$$

$$p({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k})=\frac{p({\mathbf{y}}_k|{\mathbf{x}}_k)p({\mathbf{x}}_k|{\mathbf{x}}_{k-1})}{p({\mathbf{y}}_k|{\mathbf{y}}_{1:k-1})}p({\mathbf{x}}_{0:k-1}|{\mathbf{y}}_{1:k-1})=\left[\prod _{l=1}^k\frac{p({\mathbf{y}}_l|{\mathbf{x}}_l)}{p({\mathbf{y}}_l|{\mathbf{y}}_{1:l-1})}\cdot p({\mathbf{x}}_l|{\mathbf{x}}_{l-1})\right]p({\mathbf{x}}_0).$$

$$p({\bf y}_{k}|{\bf x}_{k})\propto\exp\Big{(}-\frac{1}{2}(h({\bf x}_{k})-{\bf y}_{k})^{\rm T}{\bf R}^{-1}(h({\bf x}_{k})-{\bf y}_{k})\Big{)}$$

$${\mathbf{x}}_{l|k}^{(i)}$$

$${\mathbf{X}}_{l|k}=\left({\mathbf{x}}_{l|k}^{(1)},...,{\mathbf{x}}_{l|k}^{(M)}\right).$$

$${\mathbf{x}}_{k-L:k|k}^{(i)}\sim p({\mathbf{x}}_{k-L:k}|{\mathbf{y}}_{1:k}).$$

$${\mathbf{X}}_{k-L:k|k}=\left({\mathbf{x}}_{k-L:k|k}^{(1)},...,{\mathbf{x}}_{k-L:k|k}^{(M)}\right).$$

$$p({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k-1})\approx \frac{1}{M}\sum _{i=1}^M\delta ({\mathbf{x}}_{0:k}-{\mathbf{x}}_{0:k|k-1}^{(i)}).$$

$$w_{k|k}^{(i)}=\frac{p({\mathbf{y}}_k|{\mathbf{x}}_{k|k-1}^{(i)})}{\sum _{l=1}^{M}p({\mathbf{y}}_k|{\mathbf{x}}_{k|k-1}^{(l)})}.$$

$$p({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k})\approx \sum _{i=1}^M{w}_{k|k}^{(i)}\delta ({\mathbf{x}}_{0:k}-{\mathbf{x}}_{0:k|k-1}^{(i)})$$

$${\mathbf{X}}_{l|k}={\mathbf{X}}_{l|k-1}{\mathbf{D}}_{k|k}.$$

$$\frac{1}{M}\mathbb{E}[{\mathbf{D}}_{k|k}1]={\mathbf{w}}_{k|k},$$

$${\mathbf{w}}_{k|k}={\left(w_{k|k}^{(1)},\dots ,w_{k|k}^{(M)}\right)}^{\mathrm{T}}\in {\mathbb{R}}^{M\times 1}$$

$${\mathbf{X}}_{0|k}={\mathbf{X}}_{0|0}{\mathbf{D}}_{1|1}{\mathbf{D}}_{2|2}\cdots {\mathbf{D}}_{k|k}$$

$${\mathbf{X}}_{0:k|k}={\mathbf{X}}_{0:k|k-1}{\mathbf{D}}_{k|k}$$

$${\mathbf{X}}_{k-L:k|k}={\mathbf{X}}_{k-L:k|k-1}{\mathbf{D}}_{k|k}$$

$${\mathbf{D}}_{k|k}^{\mathrm{ESRS}}:={\hat{\mathbf{w}}}_{k|k}{1}^{\mathrm{T}}+{\text{S}}_{k|k},$$

$${\hat{\mathbf{w}}}_{k|k}=\frac{1}{M-1}{\text{S}}_{k|k}^2(\mathbf{H}{\text{A}}_{k|k-1}{)}^{\mathrm{T}}{\text{R}}^{-1}({\mathbf{y}}_k-\mathbf{H}{\mathbf{m}}_{k|k-1})$$

$${\text{S}}_{k|k}={\left\{\mathbf{I}+\frac{1}{M-1}(\mathbf{H}{\text{A}}_{k|k-1}{)}^{\mathrm{T}}{\text{R}}^{-1}\mathbf{H}{\text{A}}_{k|k-1}\right\}}^{-1/2},$$

$$\textbf{A}_{k|k-1}=\Big{(}{\bf x}^{(1)}_{k|k-1}-{{\bf m}}_{k|k-1},\dots,{\bf x}^{(M)}_{k|k-1}-{{\bf m}}_{k|k-1}\Big{)}$$

$${\mathbf{m}}_{k|k-1}=\frac{1}{M}\sum _{i=1}^M{\mathbf{x}}_{k|k-1}^{(i)},$$

$${\mathbf{D}}_{k|k}^{\mathrm{NETS}}={\mathbf{w}}_{k|k}{1}^{\mathrm{T}}+\sqrt{M}{\left[{\mathbf{W}}_{k|k}-{\mathbf{w}}_{k|k}({\mathbf{w}}_{k|k}{)}^{\mathrm{T}}\right]}^{1/2}{\mathbf{\Omega }}_{k|k}$$

$${\mathbf{W}}_{k|k}=\text{mbox}diag({\mathbf{w}}_{k|k})$$

$${\mathbf{D}}_{k|k}^{\mathrm{ETPF}}=\arg \min \sum _{i,j=1}^M{d}_{ij}||{\mathbf{x}}_{k|k-1}^{(i)}-{\mathbf{x}}_{k|k-1}^{(j)}|{|}^2$$

$$x_k={ϵ}_k,{ϵ}_k\sim \mathcal{N}(0,1)$$

$$y_k=h(x_k)+{\nu }_k,{\nu }_k\sim \mathcal{N}(0,1)$$

$$p(x_1|y_1)=\mathcal{N}(0,0.5)$$

Full text

Keywords. Data assimilation, particle smoother, optimal transport, localisation, hybrid, second order correction

AMS(MOS) subject classifications. 65C05, 62M20, 93E11, 62F15, 86A22

1 Introduction

In many applicational areas, e.g., in meteorology [Shu57], neuroscience [SQJZ18] and robotics [Sär13], there is a need to smooth model states using data extending further into the future. The underlying setting of such state estimation problems consists of two information sources: a model describing the system of interest and observations of the model state. Here we consider a discrete-time Markov state space model of the general form

[TABLE]

where ${\bf x}_{k}\in\mathbb{R}^{N_{x}}$ denotes the state of the system at time $k$, ${\bf Q}\in\mathbb{R}^{N_{x}\times N_{x}}$ is the positive semi-definite covariance matrix of model errors and the initial states ${\bf x}_{0}$ are drawn from a given distribution $p({\bf x}_{0})$. Here $\mathcal{N}({\bf 0},{\bf Q})$ stands for the Gaussian distribution with mean zero and covariance matrix ${\bf Q}$ and we also allow for the case ${\bf Q}={\bf 0}$ in which case the model (1) becomes deterministic. Additionally, partial and noisy measurements of the system states are available according to the following process:

[TABLE]

where ${\bf y}_{k}\in\mathbb{R}^{N_{y}}$ is the measurement at time $k$ for $k=1,\ldots,K$, and ${\bf R}\in\mathbb{R}^{N_{y}\times N_{y}}$ is a positive definite matrix of measurement errors. The notation

[TABLE]

is used to indicate the collection of observations from time one to $k\leq K$. The same short-hand notation will be used for collection of state vectors. The overall goal in the resulting estimation problems is usually to either approximate the marginal filtering $p({\bf x}_{k}|{\bf y}_{1:k})$ or smoothing distributions $p({\bf x}_{l}|{\bf y}_{1:k})$ for $l=0,\ldots,k-1$.111With a slight abuse of notation we use the same symbol $p$ to denote different distributions. It should be clear from the arguments in $p$ which specific distribution is referred to. Whilst filtering is useful for real-time prediction where the data comes in sequentially, smoothers are more desirable when one has the possibility to asses a whole time interval of measurements at once and can therefore obtain more informative posteriors by considering the full data set. This is particularly effective for state estimation in the context of models with long memory or long temporal correlations, e.g. ocean circulation modeling [CBV*+*10a, CV06]. In some applications, one might also be interested in the joint smoothing distribution $p({\bf x}_{k-L:k}|{\bf y}_{1:k})$ over some time window of length $L\leq k$. This is, for example, the case when the discrete-time Markov model (1) is replaced by a model of order $n>1$ of the general form

[TABLE]

in which case sequential state estimation requires knowledge of the joint smoothing distribution $p({\bf x}_{k-n:k-1}|{\bf y}_{1:k-1})$ at time $k-1$, that is $L=n$, in order to proceed to time $k$ using the model (4) and the data (2).

Due to these advantages, smoothing is of fundamental importance in several applicational areas. For example variational approaches [FLK05] such as 4DVAR [DT86] are well established in Numerical Weather Prediction and randomized maximum likelihood and iterative ensemble Kalman smoothers [BS14, CO17, Eve18] are key methods for oil reservoir modelling. Our focus is strictly on derivative-free ensemble based smoothing techniques as opposed to the large family of methods that rely on minimisation of certain cost functionals such as the aforementioned methods.

In the linear-Gaussian case, the smoothing distribution remains Gaussian with mean and covariance that can be computed by the Kalman smoother, which is a direct generalization of the Kalman filter [Kal60, Jaz70]. In the nonlinear setting, the popular family of ensemble Kalman filters (EnKFs) [Eve06] has inspired a number of different smoothers [EvL00, CBV*+*10a, WC02]. The reader is referred to [CBBE18] and references therein for a recent review on the various existing approaches.

Alternatively there are several smoothing extensions of particle filters [Sär13]. The canonical ansatz for a particle smoother is to build on the classical bootstrap filter to obtain smoothed state space estimates with respect to observations extending beyond the current time [DGA00, KS01]. However particle degeneracy is even more prominent in the smoothing case than in filtering. If the smoother is only applied over a relatively small time interval, degeneracy can be minimised significantly [CG98]. Further, a number of different implementations, e.g., using variance reduction techniques, have been proposed and a nice overview can be found in [DdFe01, Sär13]. Further solutions to the smoothing problem can also be derived using optimal control techniques by assuming a linear evolution model for the trajectories [DGA00, DdFe01]. Moreover smoothing approaches have recently been used to improve popular filtering techniques such as the auxiliary particle filter [YP11] or the implicit particle method in [WMS13]. Conceptually the idea is to utilize proposal densities that factor in observations at future time steps (typically one data point beyond the current time instance is already sufficient to improve the underlying filter).

Here we propose a unifying framework, referred to as linear ensemble transform smoothers (LETSs), to highlight the connection between various smoothing algorithms. As we will show, the family of LETSs includes a range of popular smoothers such as the class of ensemble square root smoothers (ESRS) [EvL00, CBV*+*10b, WC02] and the nonlinear ensemble transform smoother (NETS) [TA15]. We propose for the first time, a smoother extension of the ensemble transform particle filter (ETPF) [RC15], which requires a non-trivial modification of existing LETS formulations. The advantage of the ensemble transform particle smoother (ETPS) is that it has no underlying Gaussianity assumptions and is formally consistent with the true posterior in the ensemble limit (i.e., $M\rightarrow\infty$).

The other main contribution of this manuscript is to present a series of tools that can be used to address computational difficulties associated with smoothing in moderate dimensional nonlinear problems, particularly for small to moderate ensemble sizes. This is motivated by the fact that in many applications, particularly numerical weather prediction, ensemble methods with very small ensemble sizes relative to the state dimension are utilised [BTM08, MK11, NAV13]. The first tool we discuss is the so-called hybrid formulation, which allows one to develop superior smoothers by combining very robust but inconsistent smoothers such as the ESRS with consistent smoothers such as the ETPS. We demonstrate how such hybrid formulations can be developed in the context of LETSs, taking inspiration from hybrid formulations recently introduced in the context of filters in [CRR16]. Another important practical issue is that one is typically forced to use a small ensemble in many applications. A common problem associated with this in filtering is that the posterior spread is often underestimated. In practice this issue is approached via a random increase of the ensemble variance. A more judicious strategy is to match the second moment of ensemble to a spread one trusts to be better than the current approximation. With this in mind, we extend a moment matching technique described in [AdWR16] to the class of smoothers in this paper. Localisation is another popular strategy to deal with effects of spurious correlations due to small ensemble sizes [HKS07, HM01], and can also be seen as a way to handle the curse of dimensionality in the particle filtering context [FB18]. This strategy has also been utilized in the context of smoothing, e.g., see recent work in [MHP18, Boc16, CO17]. Whilst in filtering one is typically concerned with localising in space, in the more general smoothing problem it is pertinent to also consider localisation in time. In this paper, we propose a number of spatio-temporal localisation techniques specifically for LETSs. Finally, we demonstrate the benefits of these approaches through a range of numerical experiments.

The remainder of the paper is structured as follows. First an overview of the mathematical background and notation is given. In Section 3 we first summarize the key idea behind the ESRS and the NETS and demonstrate why their straightforward extension to the ETPF fails. We then introduce a generalised family of LETSs in the form of (34) which overcomes this limitation. This is followed by a discussion on numerical tools such as localisation in time, ensemble spread correction and hybrid formulations that can be employed in the context of the proposed LETSs. The presented class of smoothers is numerically investigated for a range of systems in Section 5: in particular we consider a higher order Markov system, as the effects of smoothing are particularly prominent when the underlying model has long memory. We also consider a classical moderate dimensional toy system, since finite ensemble sizes are particularly problematic for high dimensional applications. Lastly in Section 6, we discuss how spatio-temporal localisation techniques can be utilized to make the smoothing algorithms feasible in that setting.

2 Mathematical background and notation

As the smoothing algorithms we will consider throughout the paper are implemented in a sequential manner it is important to understand the following equalities: The joint smoothing distribution $p({\bf x}_{0:k}|{\bf y}_{1:k})$ can simply be linked to the previous smoothing distribution $p({\bf x}_{0:k-1}|{\bf y}_{1:k-1})$ using the assumption that model and measurement errors in (1) and (2), respectively, are independent:

[TABLE]

Here $p({\bf x}_{k}|{\bf x}_{k-1})$ denotes the transition probability associated with the Markov model (1) and

[TABLE]

denotes the likelihood function associated with the measurement model (2). The joint smoothing distribution (5) gives rise to the associated marginal distributions $p({\bf x}_{l}|{\bf y}_{1:k})$ in the state variable ${\bf x}_{l}$ for $l=k$ (filtering), $l<k$ (smoothing), or the fixed-lag smoothing distribution $p({\bf x}_{k-L:k}|{\bf y}_{1:k})$. Computing expectation values with respect to the joint smoothing distribution (5) or any of its marginal distributions is generally intractable except for the linear Gaussian case. We therefore consider sample-based Monte Carlo methods222This paper being at the crossroad of the ensemble Kalman and the sequential Monte Carlo/particle filter/smoothing literature, we use the notion of particles, samples, and ensemble members synonymously throughout this paper. in this paper and use the following notation to indicate $M$ samples ${\bf x}_{l|k}^{(i)}\in\mathbb{R}^{N_{x}}$, $i=1,\ldots,M$, representing the marginal prediction ($l>k$), filtering ($l=k$), and smoothing ($l<k$) distributions respectively, which we write with a slight abuse of notation as

[TABLE]

where $0\leq l\leq k+1$ for $k=0,\ldots,K$. Further, the full ensemble of samples at any time $k$ is collected in an $N_{x}\times M$ matrix such as

[TABLE]

One can also collect past and current state values of the $i$th ensemble member in a column vector ${\bf x}_{k-L:k|k}^{(i)}$ of dimension $(L+1)N_{x}$ for any $1\leq L\leq k$ in order to produce samples from the associated fixed-lag $L$ smoothing distributions, which we write again with a slight abuse of notation as

[TABLE]

Again the full ensemble of such trajectory samples can be collected in an $(L+1)N_{x}\times M$ matrix

[TABLE]

Below we describe a classical importance sampling implementation of the particle smoother [GDW12]. The crucial difference to filtering is that after every assimilation step, the entire history is updated according to the new importance weight. Just as in filtering, it is important to perform resampling to mitigate weight degeneracy, although in smoothing, resampling is typically undertaken trajectory-wise [Sär13]. A single cycle of a classical particle smoother at time $k$ is therefore structured as follows:

Prediction:333We mention that the prediction step can utilize alternative proposal densities [AMGC02]. However, we will not discuss such alternative proposals in this paper as they add another layer of computational complexity while generally not being able to beat the curse of dimensionality of a filtering/smoothing problem [SBBA08, SBM15, MHS17]. Given $M$ samples ${\bf x}_{0:k-1|k-1}^{(i)}\in\mathbb{R}^{kN_{x}}$ from the smoothing distribution at time $k-1$, one first produces a forecast sample ${\bf x}_{k|k-1}^{(i)}\in\mathbb{R}^{N_{x}}$ at time $k$ using (1) with ${\bf x}_{k-1}={\bf x}_{k-1|k-1}^{(i)}$ for all $i=1,\ldots,M$. This leads to the concatenated trajectory samples ${\bf x}^{(i)}_{0:k|k-1}\in\mathbb{R}^{(k+1)N_{x}}$ and the prior density for the subsequent Bayesian inference step is approximated via the corresponding empirical measure

[TABLE]

Weight computation: Importance weights are then computed according to an incoming observation ${\bf y}_{k}$:

[TABLE]

Resampling: Lastly, one performs resampling by sampling with replacement the entire state trajectories ${\bf x}^{i}_{0:k|k-1}$, $i=1,\ldots,M$, according to the weighted empirical measure

[TABLE]

instead of just the states at the current time (as is the case in filtering) [Kit96]. As a result one obtains $M$ equally weighted trajectories ${\bf x}_{0:k|k}^{(i)}$.

We remark that the resampling with replacement step can be thought of as applying an $M\times M$ transformation matrix ${\bf D}_{k|k}$ to the matrices ${\bf X}_{l|k-1}$ at each time instance $l=0,\ldots,k$, that is,

[TABLE]

The entries of ${\bf D}_{k|k}$ are drawn at random from zero and one such that ${\bf D}_{k|k}^{\rm T}\mathbb{1}=\mathbb{1}$ and

[TABLE]

where $\mathbb{1}\in\mathbb{R}^{M}$ denotes a vector of ones,

[TABLE]

denotes the vector of importance weights (12), and $\mathbb{E}$ stands for taking the expectation value. If the $(i,j)$-th entry of ${\bf D}_{k|k}$ is one, then this implies that the $i$-th ensemble member ${\bf x}_{l|k-1}^{(i)}$ at time $k$ gets replaced by the $j$-th member, that is ${\bf x}_{l|k}^{(i)}={\bf x}_{l|k-1}^{(j)}$ for all $l=0,\ldots,k$.

It is crucial to point out that a direct implementation of this classical particle smoother is known to degenerate as soon as $k$ is too large (even with resampling) [Sär13]. More specifically, given $M$ samples ${\bf x}_{0|0}^{(i)}$ from the initial distribution $p({\bf x}_{0})$, collected in the ensemble matrix ${\bf X}_{0|0}$, one finds that the above algorithm leads to

[TABLE]

and ${\bf X}_{0|k}$ will contain many identical ensemble members, leading to a largely diminished effective sample size for the smoothing distribution at initial time. Upon using the compact notation (10) with $L=k$ one can rewrite the complete trajectory resampling step in matrix form as

[TABLE]

from the columns of which one can read off the smoothed sample trajectories ${\bf x}_{0:k|k}^{(i)}$ if desired. Again, many of these trajectories will contain identical state samples ${\bf x}_{l|k}^{(i)}$ as $l\ll k$.

For this reason smoothing is often carried out for a relatively small fixed-lag $L$, that is, (18) gets replaced by

[TABLE]

with $L=0$ leading back to filtering. Consequently we will focus on estimating the joint smoothing distributions $p({\bf x}_{k-L:k}|{\bf y}_{1:k})$ or the marginal distributions $p({\bf x}_{l}|{\bf y}_{1:k})$, $l=k-L,\ldots,k$, in this paper. Furthermore, recent transformation approaches to smoothing, such as the ESRS and the NETS, can be put into the above framework (19) with different choices of the matrix ${\bf D}_{k|k}$. Most importantly, the entries of ${\bf D}_{k|k}$ are no longer drawn at random and can take real values. We will also find that an extension of the ETPF to smoothing requires an important modification to the above framework.

3 Linear Ensemble Transform Smoothers

We now introduce the concept of LETSs, using the framework developed in the previous section. LETSs provide a generalisation of the associated family of linear ensemble transform filters (LETFs) [RC15]. An LETS is any smoother in which the ensemble of samples ${\bf X}_{k-L:k|k-1}$ of the prior $p({\bf x}_{k-L:k}|{\bf y}_{1:k-1})$ are moved towards an observation via a linear transformation in the form of (19), where the transformation matrix must be of size $M\times M$. Note that $L=0$ corresponds to an LETF, and the transformation (19) describes a standard filtering update in commonly employed filters such as the stochastic EnKF [Eve06] and the ensemble square root filter (ESRF) [TAB*+*03] [RC15].

The remainder of this section is devoted to three specific examples of LETSs. We first demonstrate that the ESRS and the NETS both fit into the framework provided by (19). We then propose the smoother variant of the ETPF, showing that it is also an LETS, but that ${\bf D}_{k|k}$ in (19) requires an important extension to ensure consistency.

3.1 Ensemble Square Root Smoother

Smoother extensions of the class of ensemble square root filters [TAB*+*03] have been investigated, for example [CBV*+*10b, Fai09]. Here we briefly summarize the extension of the ESRF in [RC13] in the context of LETSs. The resulting smoother is a type of deterministic ensemble Kalman smoother, reformulated such that it can be seen as an LETS. We assume for simplicity that the forward operator in (2) is linear, that is, $h(x)={\bf H}x$. The ESRS is designed so that the empirical first two moments match the Kalman mean and covariance matrix under a linear model and observation setting. Then the coefficients of the linear transformation of the ESRS are given by

[TABLE]

with

[TABLE]

and

[TABLE]

where the matrix of ensemble deviations and the ensemble mean are defined by

[TABLE]

and

[TABLE]

respectively. The ESRS is obtained from (19) with ${\bf D}_{k|k}={\bf D}^{\rm ESRS}_{k|k}$. Note that the transformation ${\bf D}^{\rm{ESRS}}_{k|k}$ is equal to the filter transformation of an ESRF at time $k$. In other words in an implementation of the ESRS one only needs to compute the filter transformations and multiply them with past states up to lag $L$. For more details on how to derive the transformation matrix (20) see, for example, [RC15].

The ESRS is easy to implement and applicable even under relatively small samples sizes. However in general it is not consistent, that is, it does not generate samples of the true posterior distribution other than the case of linear Gaussian systems even in the limit of $M\to\infty$.

3.2 Nonlinear Ensemble Transform Smoother

Along the lines of the ESRS, the NETS [KTAN17] is also constructed to match the first two moments correctly for finite $M$. The difference to the ESRS is that the NETS is matched to the empirical mean and covariance of the classical particle filter for fixed $M$. This ansatz yields the following linear transformation

[TABLE]

where ${\bf w}_{k|k}$ is the vector of importance weights (16) at time $k$ and

[TABLE]

is a diagonal matrix with the importance weights as diagonal elements. The NETS is obtained from (19) with ${\bf D}_{k|k}={\bf D}^{\rm NETS}_{k|k}$. Again the transform matrix ${\bf D}^{\rm NETS}_{k|k}$ is equivalent to that of the corresponding NETF. See [KTAN17] for more details.

The ${\bf\Omega}_{k|k}$ in (25) represents a suitably chosen orthogonal rotation matrix subject to ${\bf\Omega}_{k|k}\mathbb{1}=\mathbb{1}$. The choice of this matrix is discussed in further detail in Section 4.3. Despite matching the first two moments of the particle smoother for a finite $M$ the filter does not converge towards the true posterior in the limit $M\to\infty$.

3.3 Ensemble Transform Particle Smoother

The ETPF relies on a transformation which involves solving an optimal transport problem designed to maximize the correlation between ${\bf X}_{k|k}$ and ${\bf X}_{k|k-1}$, whilst obtaining a set of samples distributed in accordance with the importance weights ${\bf w}_{k|k}$ [RC15]. The resulting transformation matrix, denoted by ${\bf D}^{\rm ETPF}_{k|k}$, is given by

[TABLE]

with $d_{ij}$ corresponding to the $(i,j)$-th element of the transformation matrix. Following the formulation of the ESRS and the NETS, one could be tempted to define the ETPS extension of the ETPF by setting ${\bf D}_{k|k}={\bf D}^{\rm ETPF}_{k|k}$ in (19). However such a naive approach does not take into account the geometry of the underlying trajectory samples ${\bf x}_{k-L:k|k-1}^{(i)}$ and can lead to statistically inconsistent updates even in the limit $M\to\infty$ as demonstrated by the following simple experiment.

Example: We consider the following one dimensional model with $f=0$ and ${\bf Q}=1$ in (1), i.e,

[TABLE]

The observations are according to (2) with ${\bf R}=1$ and $h(x)=x$, i.e,

[TABLE]

The purpose of this example is to demonstrate the inconsistency of the naive approach and it is thus sufficient to only consider one filtering step (at time $t=1$). Let us assume that $y_{1}=0$. Note that $x_{0}\sim\mathcal{N}(0,1)$ and the prior distribution at time 1 is $p(x_{1})=\mathcal{N}(0,1)$ as it is independent of the state at $t=0$. The filtering distribution at $t=1$ is

[TABLE]

whereas the marginal smoothing distribution at $t=0$ is

[TABLE]

Note that (31) is equal to the initial distribution $p(x_{0})$ because $x_{1}$ and $x_{0}$ are independent of each other and thus the observation $y_{1}$ does not affect $x_{0}$. However, this is not what one obtains from the naive smoothing extension of the ETPF to this problem. Instead one finds that the variance of $p(x_{0}|y_{1})$ is systematically reduced to approximately $0.5$ for $M\to\infty$. This is shown in Figure 1 which shows the variance of (31) for $60$ individual independent runs for different ensemble sizes $M\in\{10,100,200,1000\}$ for the naive implementation of the ETPS, i.e., the update is done via the ETPF update given (27). For comparison, the variance obtain via the ETPS proposed below in 32 is also displayed, which clearly converges to the correct variance as $M\rightarrow\infty$.

Instead the desired transformation matrix ${\bf D}_{k-L:k|k}^{\rm ETPS}$ should be obtained as the solution of the following linear transport problem:

[TABLE]

with the minimization performed over the set of admissible $M\times M$ matrices ${\bf D}$ defined by

[TABLE]

Here $d_{ij}$ stands for the $(i,j)$-th entry of the matrix ${\bf D}$. The ETPS is now obtained from (19) with ${\bf D}_{k|k}={\bf D}^{\rm ETPS}_{k-L:k|k}$.

The transformation of the ETPS can be interpreted as a deterministic replacement of the resampling step of sequential importance resampling. The ETPS analysis step is consistent in the sample limit $M\to\infty$ as the transformation ${\bf D}^{\rm ETPS}_{k-L:k|k}$ converges to a transfer plan $\Psi$ inducing a deterministic coupling between the random variable associated with the prior trajectories ${\bf x}_{k-L:k|k-1}^{(i)}$ and the random variable associated with the posterior trajectories ${\bf x}_{k-L:k|k}^{(i)}$. This follows from a general approximation result of [McC95].

It is further important to note that the associated coupling is maximizing the correlation between the random variable associated with the prior and the one associated with the posterior. Lastly, note that the ETPS transformation matrix in (32) is computed using both the current weights ${\bf w}_{k|k}$ and information regarding the ensemble of trajectories ${\bf x}_{k-L:k|k-1}^{(i)}$ over the fixed-lag window. Hence the transform matrix ${\bf D}^{\rm ETPS}_{k-L:k|k}$ is different from that of the corresponding ETPF at time $k$ which relies on the ensemble of states ${\bf x}_{k|k-1}^{(i)}$ only.

We demonstrate in Section 5, that it is also preferable to replace ${\bf D}_{k|k}^{\rm NETS}$ in the NETS by a sample trajectory based ${\bf D}_{k-L:k|k}^{\rm NETS}$ through an optimal choice of the rotation matrix $\mathbf{\Omega}_{k|k}$ in (25). See also Section 4.3 below. Hence we close this section by noting that

[TABLE]

provides the appropriate general LETS framework with (19) for ${\bf D}_{k-L:k|k}={\bf D}_{k|k}$ as a special case.

4 Numerical tools for moderate sample sizes

Smoothing is naturally a higher dimensional problem than the associated filtering problem and thus it is important to address this additional complexity in the context of moderate sample sizes. In other words, smoothers typically require additional techniques to make them feasible for practical estimation problems. We will discuss three approaches in this section that can help to significantly improve the accuracy and computational feasibility of the proposed family of smoothers especially in the small to moderate sample size setting. We also discuss how the LETS framework is beneficial in such developments, since such techniques do not need to be developed for each algorithm individually, but can be applied to any smoother that can be seen as an LETS. Further techniques for spatially extended systems will be discussed in Section 6.

4.1 Temporal Localisation

Localisation has become a powerful tool for beating the curse of dimensionality for sequential filtering algorithms such as the ensemble Kalman filter. These techniques have recently been extended to particle filters [CR15, Pot15, RK17] and to smoothing algorithms. For example a weight localization procedure for use in Variational Particle Smoothers is proposed in [MHP18], where the weights are computed for each state component independently (or some block form of it). [Boc16] investigate how a localization operator can be modified over the smoothing window using the underlying dynamics of the system.

We note that the fixed-lag smoother can already been viewed as a temporal localisation of the global update step (18). Still, the optimal transport problem over an $(L+1)N_{x}$ dimensional space, as defined by (32) and (33), becomes high-dimensional for time-lags $L\gg 1$ even if the dimension of state space, $N_{x}$, is moderate. Hence the ETPS might perform poorly under those circumstances even if the associated ETPF does well on a given problem. This degradation of performance can be counteracted by a further localization step in time. More precisely, we replace the single $(L+1)N_{x}$ dimensional transport problem by a sequence of $N_{x}$-dimensional transport problems

[TABLE]

$l=k-L,\ldots,k$, with the minimization performed over the $M\times M$ matrices ${\bf D}$ subject to (33). We then replace (14) by

[TABLE]

for $l=k-L,\ldots,k$. We note that the localized update (36) is now consistent with the ETPF filtering update since ${\bf D}_{k|k}^{\rm ETPF}={\bf D}_{k|k}^{\rm ETPS}$. However, while the marginal distributions $p({\bf x}_{l}|{\bf y}_{1:k})$ are correctly approximated by (35), the correlation between samples ${\bf x}_{l|k}^{(i)}$ at different time-levels $l,l^{\prime}\in\{k-L,\ldots,k\}$ can be misrepresented, which can lead to problems when applied to a model of the form (4) of order $n$, $n>1$, where the joint distribution of state variables over a time interval $[k-n,k-1]$ matters in the prediction step. We explore this issue in Section 5 in more detail.

Having to solve $L+1$ optimal transport problems (35) can become computationally demanding. Our numerical experiments for the Lorenz 63 system from Section 5 suggest that this cost can be reduced by using

[TABLE]

provided the sample paths are strongly correlated in time, as is the case for many deterministic systems (i.e. with ${\bf Q}=0$). This issue is investgated further in Section 5. It should be mentioned that the update formulation of LETSs in (19) involves modifying the ensemble at each time $l$ in the smoothing window. When applied to deterministic systems, the resulting ensembles will not necessarily satisfy the underlying system of equations. This issue is well known in the area of ensemble Kalman filtering/smoothing, and methods to ensure system balances are maintained is an active area of research. One could also adopt an alternative approach whereby the goal is to estimate the smoothing update at the beginning of the window only, similar to the iEnKS [BS14] and others.

4.2 Hybrid Smoothers

Ensemble smoothers with an underlying Gaussian assumption are often favoured over consistent smoothers, as they are typically more robust in high dimensional applications. Yet depending on the problem, a Gaussian approximation can lead to a very crude estimate of the true posterior. One way to improve accuracy while maintaining some level of robustness is to combine different smoothers. This strategy has proven to be rather successful in the filtering context [FK13, PWR19]. Our hybrid formulations are based on the idea of a split likelihood [CRR16]:

[TABLE]

where $\alpha\in[0,1]$ is a tuning parameter that determines the contribution of each smoother. The benefit of the LETS framework is that it allows one to develop a hybrid smoother using any pair of smoothers that can be written in the form of (19). Specifically, this is done by a two-step transformation to form a mixture that benefits from desirable properties of each individual smoother. For example, a hybrid of the ETPS and ESRS would be developed as

[TABLE]

Here ${\bf D}^{\rm{ETPS}}_{k-L:k|k}(\alpha)$ denotes the transformation matrix corresponding to the tempered likelihood $p({\bf y}_{k}|{\bf x}_{k})^{\alpha}$ and ${\bf D}^{\rm{ESRS}}_{k|k}(1-\alpha)$ to $p({\bf y}_{k}|{\bf x}_{k})^{(1-\alpha)}$, respectively. Note that $p({\bf y}_{k}|{\bf x}_{k})^{\alpha}$, $\alpha\in(0,1]$, for given likelihood (6) corresponds to replacing the measurement error covariance matrix ${\bf R}$ by ${\bf R}/\alpha$ in the ESRS and ETPS formulations. The hybrid example in (39) is particularly valuable as two LETSs with different strengths are combined yielding a filter superior to the individual LETSs (for optimal $M$ and $\alpha$) [CRR16]. It has also been found numerically that lower RMSEs are achieved when the ETPF is used first [CRR16]. The choice of $\alpha$ determines the influence of each smoother, with either [math] or $1$ indicating purely one of the smoothers. The optimal $\alpha$ can be tuned using an appropriate performance metric, for example, the RMSE.

4.3 Adaptive spread correction and rotation

In practice, LETSs are often employed in a setting far from the asymptotic limit ($M\rightarrow\infty$) which often results in an underestimation of the spread. This underestimation is typically counteracted by the use of artificial inflation of the sample spread (sometimes also referred to as particle rejuvenation). These procedures are usually tied to a tuning factor that needs to be adjusted.

An alternative approach is to correct the sample spread adaptively by matching it to the second moment of the classical particle smoother in order to improve the overall accuracy of LETSs. This ansatz is justified by the fact that the classical particle smoother converges weakly to the true posterior and thus can be a good estimate of the true measure (if $M$ is not too small). We describe this approach below. Any smoother belonging to the subfamily of LETSs (34) with transformation matrices ${\bf D}_{k-L:k|k}\in\mathcal{D}_{1}$,

[TABLE]

match the first moment of the particle filter over a single assimilation step. Here the vector of weights ${\bf w}_{k|k}$ is defined by (16).

For example, the ETPS as well as the smoother associated with the simple transformation ${{\bf D}}_{k-L:k|k}={\bf w}_{k|k}{\mathbb{1}}^{\rm T}$ are in $\mathcal{D}_{1}$. Provided ${\bf D}_{k-L:k|k}\in\mathcal{D}_{1}$, it follows that this leads to equivalence of the empirical means, that is,

[TABLE]

However smoothers (34) with transformation ${\bf D}_{k-L:k|k}\in\mathcal{D}_{1}$ do not necessarily match the second moment, that is,

[TABLE]

is not equal to

[TABLE]

For (42) and (43) to be equal the following equation has to be satisfied

[TABLE]

where ${\bf W}_{k|k}$ is defined in (26).

Remark It is important to note that for small ensemble sizes the spread given by (42) associated with importance sampling is expected to be underestimated as well. Yet ultimately the ETPS is underestimating the spread even stronger and in the context of filtering it has been confirmed numerically that an adjustment to the particle filter significantly improves the results even for a small number of particles [AdWR16]. Thus an adjustment of the ETPS to the particle smoother spread can be a beneficial modification even when far from the ensemble limit.

The transformation can be corrected via an additive term ${\bf{\Delta}}_{k|k}$ that vanishes in the particle limit $M\to\infty$ similar to the adjustment in the filtering case [AdWR16]. That is, for any ${{\bf D}_{k-L:k|k}}\in\mathcal{D}_{1}$ the associated second order correct transformation is computed as

[TABLE]

where ${\bf\Delta}_{k-L:k|k}\in\mathbb{R}^{M\times M}$ is chosen such that ${\bf\Delta}_{k-L:k|k}^{\rm T}\mathbb{1}={\bf 0}$, ${\bf\Delta}_{k-L:k|k}\mathbb{1}={\bf 0}$, and is a solution of the algebraic Riccati equation

[TABLE]

Such Riccati equations are in general non-trivial to solve, however for the purposes of adjusting the transformation it is often sufficient to approximate the solution by numerically solving the corresponding dynamic Riccati equation [AdWR16]. The above developments demonstrate the usefulness of the LETS framework, as one only needs the form of the transformation matrix (as well as the particle smoother weights) in order to construct a second order correction. Further, it is important to note that

[TABLE]

where ${\bf\Omega}_{k-L:k|k}$ is an $M\times M$ orthogonal matrix with the property that ${\bf\Omega}_{k-L:k|k}\mathbb{1}=\mathbb{1}$ is also a second order accurate transformation since ${\bf\Delta}_{k-L:k|k}{\bf\Omega}_{k-L:k|k}$ is also a solution to (LABEL:riccati). This fact is particularly relevant as the second order corrected version of ${{\bf D}}_{k-L:k|k}={\bf w}_{k|k}{\mathbb{1}}^{\rm T}$ leads to the NETS (25) and the correction term takes the form

[TABLE]

In [AdWR16] it was shown how the performance of the NETF ($L=0$) can be improved significantly by a suitably chosen ${\bf\Omega}_{k-L:k|k}$. This approach can be extended to the smoothing problem whereby ${\bf\Omega}_{k-L:k|k}$ is chosen so that ${\bf D}_{k-L:k|k}^{\bf\Omega}$ is the minimizer of

[TABLE]

over the set of suitable rotation matrices ${\bf\Omega}$ for ${\bf\Delta}_{k-L:k|k}$ given by (48).

The desired rotation matrix can be determined via the singular value decomposition of an $M\times M$ matrix, that is,

[TABLE]

with the $(L+1)N_{x}\times M$ matrix of ensemble deviations defined by

[TABLE]

Then the optimal rotation in (47) for ${\bf D}_{k-L:k|k}={\bf w}_{k|k}{\mathbb{1}}^{\rm T}$ is given by

[TABLE]

The reader is referred to [AdWR16] for further details.

In Section 5, we demonstrate how an optimal rotation of this form can significantly improve the accuracy of the NETS with random rotation introduced in [KTAN17].

Two fixed lag smoothers belonging to the family of LETSs, namely the ETPS and the NETS with optimal rotation, are summarized in form of pseudocodes below. More specifically, the ETPS given in Algorithm 1 combines the ETPS of Section 3.3 with the additional feature of a second order correction as described above. In order to compute a solution to the optimal transport problem given in step (6) of Algorithm 1, one can employ available solvers such as the Earth Mover’s Distances algorithm given in [PW09]. The associated computational complexity is of order $\mathcal{O}(M^{3}\ln(M))$. For the theoretical aspects associated with optimal transport problems please see volumes such as [Vil03, Vil09] and references therein. In the context of filtering the optimal transport problem has been treated in [RC15]. As it is computationally expensive to solve the optimal transport problem it is possible to utilize the so called Sinkhorn approximation (see [Cut13]) of the actual solution ${\bf D}^{\rm{ETPS}}_{k-L:k|k}$ with computational complexity of order $\mathcal{O}(M^{2}C(\lambda))$. The Sinkhorn approximation involves adding a regualrisation term to the underlying transport problem, i.e,

[TABLE]

where $\lambda>0$ is a regularisation parameter that must be tuned and $d_{ij}^{0}$ are the entries of ${\bf D}^{0}_{k-L:k|k}={\bf w}_{k|k}\mathbb{1}^{T}.$ Note that this form of approximation does not ensure that the smoother is consistent in the ensemble limit anymore. Further it is important to point out that the an update with a transformation computed via the Sinkhorn approximation underestimates the spread even more than with the original ETPS transformation. Thus it is prudent to pair the Sinkhorn approximation with the the adaptive spread correction. Although the computational complexity of solving (6) is not significantly affected by the smoother lag $L$ as the transformation matrix has dimension $M\times M$ regardless of the value of $L\geq 1$, the computational burden of solving the optimal transport problem is significant. Thus an accuracy complexity trade-off in form of the Sinkhorn approximation is justifiable. For a more detailed discussion in the context filtering we refer to [AdWR16].

The NETS has already been introduced in [KTAN17] and the novelty in the disclosed pseudocode (see Algorithm 2) is the use of an optimal rotation matrix ${\bf\Omega}_{k-L:k|k}$ in (47) for ${\bf D}_{k-L:k|k}={\bf w}_{k|k}\mathbb{1}^{\rm T}$ and ${\bf\Delta}_{k-L:k|k}$ given by (48).

5 Numerical examples

In order to complement the above derivations we conducted three sets of numerical experiments each targeting a different investigational purpose.

5.1 Classical: Lorenz 63

The Lorenz 63 system [Lor63] is a deterministic ordinary differential equation in $\mathbb{R}^{3}$ with ${\bf x}_{k}=(x_{k},y_{k},z_{k})^{T}$ given by

[TABLE]

Forward Euler time-stepping with step-size $\Delta t$ leads to a system of the form (1) with $f({\bf x}_{k})={\bf x}_{k}+\Delta t\psi({\bf x}_{k})$ and ${\bf Q}={\bf 0}$. It is a classical toy example in the context of data assimilation as it is chaotic and solutions diverge rapidly without any filtering updates. Here we mainly want to compare some of the LETSs equipped with the various implementation variants proposed in this paper. In detail a standard implementation of the ESRS, NETS with different rotations (random [TA15] and optimal as proposed in Algorithm 2 based on the ideas in [AdWR16]) and the ETPS (with and without second order correction see right panel of Figure 2). Further a hybrid smoother combining the ESRS with the ETPS is applied to the system. For the computations we consider the case where only the first component is observed, at intervals of $\Delta t_{obs}=0.12$ with observation error variance $R=8$. The discretised system is evaluated using $\Delta t=0.01$. Additionally, the initial conditions are uncertain, with an ensemble generated from $N({\bf x}_{0},0.5{\bf I})$. Due to the deterministic model dynamics, particle rejuvenation is carried out as follows:

[TABLE]

where $\beta$ is a tuning parameter taken to be 0.2 in this experiment and ${\bf\xi}_{j}$ are i.i.d Gaussian variables with mean ${\bf\mu}={\bf 0}$ and covariance matrix ${\bf I}$. This choice of parameters has been checked to be robust to small variations and is a reasonable choice compared to values used in other studies.

A sequence of $K=10,000$ observations is generated and assimilation carried out using 50 identical experiments for each algorithm, due to the randomness associated with rejuvenation. Smoother performance for a given fixed lag $L$ is assessed using the time averaged RMSE:

[TABLE]

where $K$ is the length of assimilation period (i.e., RMSE is evaluated only at time points where observations are available); ${\bf x}^{\text{ref}}_{k-L}$ refers to the true state at time $k-L$ and ${\bf\hat{x}}_{k-L|k}$ is a summary statistic of the smoothed ensemble at time $k-L$ given data up to time $k$. We consider the traditional RMSE computed on the ensemble mean, denoted $\text{RMSE}_{\mu}$ with ${\bf\hat{x}}_{k-L|k}={\bf m}_{k-L|k}$. Since the densities in the nonlinear examples in this section may be skewed or display multi-modality, we also assess performance by a RMSE as per (61), but computed on the sample mode rather than the ensemble mean ${\bf m}_{k-L|k}$, denoted $\text{RMSE}_{mo}$. The sample mode is estimated from the maximum of a kernel density estimate fitted to the marginal ensembles at each time. Finally, we examine the time and space-averaged Continuous Ranked Probability Score (CRPS) [Her00] for each fixed lag $L$ to examine reliability and resolution of the posterior ensembles.

The main interest here is the performance of various smoothers for small to moderate ensemble sizes. Figure 2 shows the time averaged RMSEs and CRPS for the fixed lag $L=6$, where RMSE is calculated at observation time points only. The optimal fixed lag (i.e. with lowest $\text{RMSE}_{\mu}$) was 5 or 6 across all ensemble sizes $M=15$ to 35. In this experiment, 1 lag unit is equal to 1 assimilation interval. For the ETPS, larger values of lamba were shown to provide improved RMSE values (results are shown for $\lambda=40$, and $\lambda=100$ for the ETPS and ETPS with no second order correction respectively). Note that the choice of $\lambda$ also affects run time, i.e., large $\lambda$ values increase run time but also improve the approximation of the true optimal transport problem. In a real world application, an ideal lambda would need to be determined based on a desired balance between runtime and accuracy. As the number of ensemble members is not greater than $M=40$ we do not choose $\lambda>40$. The only exception is the ETPS run without spread correction in the right panel of Figure 2 where the regularisation is chosen to be $\lambda=100$ to indicate that even for larger $\lambda$, one is not able to obtain good RMSEs without the spread correction.

Firstly, this experiment demonstrates the potential of the rotation techniques discussed in Section 4.3. The NETS with optimal rotation provides improved results compared to the random rotation proposed in [KTAN17], to the extent that performance is similar to the 2nd order corrected ETPS. Consistent with filtering results presented in [AdWR16], the ESRS also provides the best smoother performance for low ensemble sizes, although the ETPS and NETS provide lower RMSEs for larger ensemble sizes both in terms of the ensemble mean and mode. Such improvements can be attributed to the fact that the ETPS and NETS are better equipped to handle nonlinear relationships between observed and unobserved variables, due to the use of importance sampling rather than the linear correlation in the Kalman update. However, there is a delicate balance between bias and variance here. Larger ensemble sizes are required for the NETS and ETPS to become comparable to the ESRS in terms of CRPS. This can be attributed to the slight underdispersiveness in the ETPS and NETS ensembles in this experiment. Similar to the filtering case, the ETPS without second order correction gives the worst performance of all the smoothers. We also consider a hybrid smoother formulation as described in Section 4.2, where a smoother update is first computed using the second order accurate ETPS which is then used as the prior for an update using the ESRS. Figure 3 shows $\text{RMSE}_{\mu}$ values for the ETPS-ESRS hybrid smoother at optimal lag $(L=6)$ when the ETPS update is computed prior to the ESRS update. Improved $\text{RMSE}_{\mu}$ values compared to both the standard ETPS $(\alpha=1)$ and ESRS $(\alpha=0)$ are achieved in the hybrid framework.

5.1.1 Optimal lag

Due to sampling errors, it is often found that the optimal smoother performance is obtained for finite lags. This phenomenon is often observed for the Lorenz systems when the ensemble size is small (see [KTAN17] for an example with the Lorenz 96 system [Lor96]) and also Figure 4 where the $\text{RMSE}_{\mu}$ increases for $L>6$. This may also be related to the autocorrelation properties of the underlying system; if the autocorrelation is small or negligible after a finite time, then it is possible that smoothing beyond this time lag will have minimal benefits. These autocorrelation characteristics may be affected by the underlying chaotic properties of the system and how perturbations propagate in time. Note the correspondence between the optimal smoother lag for the Lorenz 63 $(L=5$ to $6)$ and the autocorrelation lag $(\tau)$ beyond which the autocorrelation function (ACF) of the $x$ variable is small $(<0.2)$ (see Figure 5). It is worth noting that the ACF provides only a measure of linear correlation and therefore cannot be used alone to determine the optimal smoother lag. The potential of the localisation schemes can also be seen in Figure 4 (right). The extreme localisation scheme of (37 (see temp loc const.) provides significantly improved $\text{RMSE}_{\mu}$ at longer time lags. We note that in practice, one would not continue to smooth past the optimal lag time and that this example is merely for demonstration purposes.

5.2 Long Memory: Mackey-Glass model

An important application of smoothing is in non-Markovian systems. These systems typically have strong temporal correlations in the state variables, meaning that fixed lag smoothing over a long window has the potential to improve state estimation compared to filtering. To investigate this further, we examine the Mackey-Glass model (originally proposed by [MG77], which is a single-variable non-linear delay differential equation used to model various physiological processes such as blood cell production. It is particularly useful for such systems because the evolution of the state is known to depend on the state at previous times, not just the current time. The system has infinite degrees of freedom since it is a delay differential equation, however its strange attractor is of finite dimension. Here, we consider the following parameterization of the model which produces periodic and chaotic solutions

[TABLE]

with $\phi=0.2,\gamma=0.1,\kappa=10$ and $\nu=17$. We use a 4th order Runge-Kutta scheme with $\Delta t=0.1$ to discretise the system, resulting in a finite difference system in the form of (4) with $n=\frac{\nu}{\Delta t}=170$. Observations of the state variable are available at intervals of $\Delta t_{obs}=8$, with observation error variance $R=0.05$. Here we consider the case where observations are temporally sparse since the system consists of only a single state variable. The initial condition samples $x_{0:n}^{(i)}$ are generated by perturbing $x_{0}$ with random noise from $\mathcal{N}(0,0.1)$ and then propagating it through the model until time index $n$. As with the Lorenz 63 system, particle rejuvenation using (60) with $\beta=0.005$444This choice is robust with respect to small changes. is carried out since the model dynamics are deterministic. A sequence of $K=12,000$ observations is generated and assimilation carried out using 50 identical experiments for each smoothing algorithm. The ETPS was solved using the Sinkhorn approximation (53) with $\lambda=40$ (see discussion in L63 section). In this experiment we consider only the NETS with optimal rotation proposed in Section 4.3.

Figure 6 shows the time averaged RMSEs and CRPS of the state variable against fixed lag $L$ for various smoothers for $M=50$. In this experiment, 1 lag unit is equal to 1 model time step, since the system is higher order Markov. RMSE is also evaluated at every time step rather than only when observations are available, so that we can evaluate performance at the hidden states also. Once again, the ETPS algorithm provides a significant improvement over the ESRS, and a slight improvement over the NETS with the optimal determinstic rotation scheme, in terms of all performance metrics. There is also an overall trend in reduction of all metrics with increasing fixed lag size, unlike in the Lorenz 63 system. This can be attributed to the presence of strong temporal correlations in the system even at long lag times (see Figure 5). Figure 6 (right) shows that the temporal localisation schemes discussed in Section 4.1 leads to only a minor loss of accuracy, with the computationally cheaper scheme of (37) still providing competitive results compared to the pathwise ETPS and localisation scheme of (35). Finally, note that RMSE does not decrease monotonically with lag because the presence of a higher order memory term in the model equation means that smoother solutions for lag $L$ are not solely a deterministic transformation of smoother solutions at lag $L-1$.

6 Spatially extended systems

In many applications the state space of the considered systems can be very high dimensional which can be challenging if the ensemble size is considerably smaller, i.e., $N_{x}\gg M$. Localisation techniques have been successfully employed in these settings [HKS07, HM01, NSBG14] to address this impacts of small ensemble sizes with an extended state space. In Section 4.1 we already highlighted one form of temporal localisation for the ETPS which has been numerically investigated in Sections 5.2 and 5.1. Here we will further explore localisation and look into options that allow to localise jointly in space and time for the various smoothers mentioned in this manuscript. Moreover an atmospheric toy model often used to investigate the feasibility of newly developed methods for extended state spaces is considered.

6.1 Temporal and Spatial Localisation

In the following we restrict our attention to problems where the state vector is defined on a spatial domain that is a set of discrete points in the Euclidean space $\mathbb{R}^{n}$, denoted by $\{{\bf p}_{j}\}_{j=1:N_{x}}$. That is, the $s$-th component of the vector ${\bf x}_{k}$, denoted ${\bf x}_{k}[s]$, gives the value of some variable at the position ${\bf p}_{s}\in\mathbb{R}^{n}$. For example, this may be a scalar PDE in one space dimension (i.e., ${\bf p}\in\mathbb{R}$) or on a two dimensional grid. Similarly, the observation vector is defined on a spatial domain in the same Euclidean space, but potentially on a different set of discrete points $\{{\bf q}_{j}\}_{j=1:N_{y}}$. We also assume that the observation locations are fixed in time and that observation errors are independent in space and time (i.e. ${\bf R}$ is diagonal). Furthermore, we use the notation ${\bf z}[i]$ to indicate the $i$-the component of the vector ${\bf z}$ while ${\bf z}(r)$ indicates that the vector ${\bf z}$ is a function of the variable $r$.

The key idea behind localisation is that spurious temporal correlations [NSBG14] stemming from sampling errors due to small ensemble sizes relative to the dimension of the state space are artificially reduced [Eve06]. Typically localisation is implemented by either adjusting the empirical covariance matrix (which is referred to as ${\bf B}$-localisation in the context of NWP [HM01, HWS01]) or by localising in observation space, known as ${\bf R}$-Localisation. Direct adjustment of the covariance matrix is not extendable to particle filters thus we will focus on localisation in observation space.

The general motivation behind ${\bf R}$-Localisation is that most physical systems exhibit a spatial and temporal decay in correlation. It is thus pertinent to determine whether there is a decay in correlation, and if so, how variables are jointly correlated across time and space. This is clearly highly dependent on the characteristics of the dynamical system, so it is crucial to introduce a distance measure $d_{l}({\bf p}_{s},{\bf q}_{u})$ between the points ${\bf p}_{s}$ and ${\bf q}_{u}$ depending on time $l$ tailored to the underlying physical model.

Distance in space as well as spatial correlations between state and observation components are often relatively easily to determine for many systems of interest. However temporal correlations are much harder to identify as they depend on the dynamical flow of the system, making localisation for smoothers more difficult. Yet there have been several suggestions in recent literature on how to approach this problem. For example, in [Boc16] the temporal localisation function is governed by a Liouville equation for the iterative ensemble Kalman smoother [BS14], while a weight localisation technique [Pot15] based on the distance of the current state to the observation is presented for a variational particle smoother in [MHP18]. A scheme that allows to adaptively incorporate how the covariance is affected by the flow is discussed in [BH11] in form of ${\bf B}$-localisation.

Here we suggest two choices of distance measures for ${\bf R}$-Localisation that could be used in the context of smoothing. As the optimal-lag is often relatively short, one ansatz is to assume that the spatio-temporal dependence does not change significantly within the smoothing interval. With this idea in mind, one option is to keep the metric stationary in time i.e.,

[TABLE]

where $m(\cdot,\cdot)$ refers to an appropriate distance metric, for example, the Euclidean distance, which also incorporates the underlying spatial boundary conditions of the system. This ansatz is similar to the localisation function chosen in [MHP18]. An alternative approach is to take the spatio-temporal dynamics into account in the localisation function. One way of doing so is to consider the autocorrelation structure of the underlying system. Given a time series of model simulations $\{\mathcal{X}_{t}\}_{t=1:T}$, recall that the sample autocorrelation matrix at autocorrelation lag $\tau$ is defined by

[TABLE]

where ${\bf\Sigma_{\tau}}$ is the sample autocovariance matrix given by

[TABLE]

with $\hat{\mu}=\frac{1}{T}\sum_{t=1}^{T}\mathcal{X}_{t}$, and ${\bf V}$ is a diagonal matrix with the $i$-th entry on the diagonal given by $\sqrt{{\bf\Sigma}_{0}[i,i]}$.

Our proposed localisation scheme is to take the distance between any given observation point ${\bf q}_{u}$ and the spatial grid point ${\bf p}_{s^{*}}$ at time $k$, which is the grid point that is maximally correlated in magnitude with ${\bf p}_{s}$ at time $l$, i.e.

[TABLE]

[TABLE]

The smoothing lag variable is directly related to the autocorrelation lag, i.e. $l\equiv k-\tau$, where $\tau\in\{0,1,\cdots,L\}$. Furthermore, notice that (66) is independent of $k$.

The next step of ${\bf R}$-Localisation is to adjust ${\bf R}^{\rm{loc}}_{k}(g_{sl})^{-1}$ according to the influence of the observation ${\bf y}_{k}$ on the individual components of the state ${\bf x}_{k-L:k}[g_{sl}]$ for all $s=1,\cdots,N_{x}$ and $l=k-L,\cdots,k$, where $g_{sl}$ is the index of the vector ${\bf x}_{k-L:k}$ corresponding to the $s$-th position coordinate at the $l$-th time in the smoothing window $k-L$ to $k$. More specifically

[TABLE]

where $C_{k}(g_{sl})$ is an $N_{y}\times N_{y}$ diagonal matrix with zeroes on the off-diagonal elements (since measurement errors are assumed independent) with the $u$-th diagonal element equal to

[TABLE]

for all $u=1,\cdots N_{y}$. $\rho$ is a localisation function (e.g., Gaspari-Cohn [GC99]) and $r_{loc}$ is the localisation radius which must be specified a priori. For the NETS and the ETPS the localised inverse error covariance ${\bf R}_{k}^{\rm{loc}}(g_{sl})^{-1}$ enters through the unnormalized weights

[TABLE]

for all $s=1,\cdots,N_{x}$ and $l=k-L,\cdots,k$. This form of localisation is an extension of localisation in the ETPF [RC13, CR15]. There have been several other suggestions on how to localise particle filters [Pot15] as well as variants of the particle smoother [MHP18]. Note that localisation does not only help to mitigate particle degeneracy, which generally occurs as soon as the smoothing interval is too large, but it can also help to reduce run time. This is because the update for R-localisation is done individually for each component (or for smaller blocks) and can thus be parallelised. Furthermore, note that the ensemble update via ${\bf\tilde{D}}_{l|k}(g_{sl})$ is now done individually for each component $g_{sl}$ of the state vector ${\bf x}_{k-L:k|k}$ which can mean reduced computational complexity for certain smoothers. For example, in the case of the ETPS, the aforementioned localisation scheme is combined with the localisation scheme in Section 4.1, which results in solving the following one dimensional optimal transport problem

[TABLE]

subject to

[TABLE]

This can be solved via a sorting algorithm which considerably decreases the run time [RC13]. Similarly, for the ESRS, such a localisation scheme means that the smoother can be parallelised since the update is carried out for each ${\bf x}_{k-L:k|k}[g_{sl}]$ independently. Consequently localisation is often deliberately employed when dealing with high-dimensional state spaces and is thus a crucial feature to make smoothing feasible.

6.2 Lorenz 96

The last numerical example we consider is the Lorenz 96 system (L96) [Lor96] which is given by

[TABLE]

where ${\bf x}[{-1}]={\bf x}[{N_{x}-1}]$, ${\bf x}[0]={\bf x}[{N_{x}}]$ and ${\bf x}[N_{x}+1]={\bf x}[1]$ for each spatial component $s$ and here $N_{x}=40$. It is another classical toy model which is often utilised to test algorithms in a higher dimensional context. The following setting is used for the Lorenz 96 runs. We observe every 2nd grid point in space, with temporal observation intervals of $\Delta t_{obs}=0.11$, observation error variance ${\bf R}=8$ and forcing constant $F=8$. A forward Euler scheme is used to discretize the model, with time step $\Delta t=0.005$. Additionally, the initial conditions are uncertain, with an ensemble generated from $\mathcal{N}({\bf x}_{0},0.5{\bf I})$. Specifically we would like to investigate the spatial as well as the time-wise localisation scheme of the proposed smoothers in this section. The proposed distance functions (63) and (66) are compared by means of different ESRS runs where the ensemble size is set to $M=30$ and $500000$ data assimilation steps are computed for each lag ($L=0$ to $L=8$). Further two different localisation radii: ${r}_{loc}=1$ and ${r}_{loc}=8$ are considered to show the effect of a larger radius compared to a very small one. In practice this value is often chosen according to the underlying dynamics of the system. The distance metric is chosen as $m({\bf p}_{s},{\bf q}_{u}):=min\{|s-u|,|s-u-N_{x}|,|s-u+N_{x}|\}$, which ensures that the periodic boundary conditions are taken into account. The resulting RMSE values are displayed in the panels of Figure 9. Finally, the Gaspari-Cohn function [GC99] is used for $\rho$.

The autocorrelation localisation scheme using the distance function given in (66) performs equally well as the stationary distance defined in (63) for lags [math] to $2$ for both radii due to the fact that $s^{*}=s$ (see Figure 7). Yet as soon as $s^{*}$ and $s$ differ, the autocorrelation localisation scheme produces lower RMSEs than the stationary localisation (see Figure 9). Naturally the choice of radius also effects the performance of the two schemes and the benefits of the localisation scheme based on autocorrelation coefficients is more prominent for radius $1$ as can be seen in the left panel of Figure 9. Thus in particular in a set up with small spatial localisation radii and longer smoother lags, the autocorrelation scheme would be preferred to the stationary one. Yet the stationary scheme is a valid choice for small lags and large spatial localisation radii and does not require the computation of the autocorrelation of the underlying dynamical model.

Finally a hybrid of an ESRS and the second order corrected ETPS is tested on the L96 system. For the ETPS-ESRS hybrid runs we set $L=10$ and ${r}_{loc}=8$, $\lambda=40$ and consider the stationary localisation scheme based on the distance given in (63). A total of 50,000 assimilation steps are performed with a burn-in period of $1000$ time instances. Figure 8 displays the results for $\alpha\in\{0,0.1,\dots,0.9,1\}$ for $M=30$ and $40$ and the best RMSE values are obtained for $\alpha=0.4$ in the smoothing case. This suggests that a hybrid smoother with localisation could be a useful tool to increase the accuracy of the underlying state estimation problem in a high dimensional setting.

7 Conclusions

A general formulation (34) for smoothers induced by the family of LETFs is derived. Smoothing extensions for some members of this class were already available e.g., ESRS [CBV*+*10a] or NETS [KTAN17], yet this general extension to a complete class of smoothers allows one to universally investigate the underlying properties of this family rather than having to study them individually. Furthermore, we show how the LETS framework allows one to develop the extension of the ETPF [AdWR16] to the smoother case. Important features for high dimensional state and parameter estimation such as localisation, hybrid implementations and second order corrections are introduced for the LETSs and are studied numerically.

Application to the Lorenz 63, Mackey-Glass model and Lorenz 96 systems confirmed various improvements that can be achieved via the proposed techniques. Firstly, the proposed second order corrected ETPS has superior performance in terms of RMSEs and CRPS for moderate ensemble sizes. Similar to previous studies, the ESRS has superior performance for very low ensemble sizes (relative to the state dimension), where it is known that Gaussian assumptions can lead to better tracking capability than attempting to capture the full posterior via importance sampling. This can be seen in the 40-variable Lorenz96 ($M=30$ & $M=40$) and Lorenz63 ($M\leq 20$) experiments. Secondly, the proposed numerical techniques are shown to lead to improvements on existing smoothers. The optimal rotation proposed here significantly improves RMSE scores of the NETS, such that it is almost equivalent to the RMSE of the second order corrected ETPS (see Figure 2) . In line with earlier results obtained in the context of filtering, the second order corrected ETPS has a significantly higher accuracy than the ETPS without any correction (see Figure 2, right panel). Additionally, the optimal combination of the ESRS and ETPS through the hybrid framework has better performance in a mean squared error sense than the individual smoothers (see Figures 4 and 8). This suggests that an optimal choice of hybrid smoother could be used to provide more accurate and robust posterior estimates in high dimensional nonlinear applications.

It is crucial to mention however that additional computations associated with smoothing lead to an increase in run time compared to filtering. It is therefore necessary to find the right trade off between accuracy and computational cost. In this regard, we have proposed a number of temporal and spatio-temporal localisation schemes and evaluated them numerically for the Lorenz 96 model (see Figure 9). The results are promising and suggest that smoothing (in particular a hybrid smoother) is a feasible option through appropriate localisation in a high dimensional setting.

Acknowledgement

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1294 “Data Assimilation”, Project (A02) “Long-time stability and accuracy of ensemble transform filter algorithms”.