The Ensemble Kalman Filter: A Signal Processing Perspective

TL;DR

This paper reviews the ensemble Kalman filter (EnKF) from a signal processing perspective, explaining its derivation, challenges, and relations to other filters, aiming to facilitate its adoption and inspire new research directions.

Contribution

It provides a comprehensive, signal processing-oriented overview of EnKF, including derivations, extensions, and a survey of literature, tailored for signal processing researchers.

Findings

EnKF handles high-dimensional, nonlinear, non-Gaussian problems effectively.

The paper clarifies the relation between EnKF, sigma-point KF, and particle filters.

Simulation examples demonstrate EnKF's practical performance on benchmark problems.

Abstract

The ensemble Kalman filter (EnKF) is a Monte Carlo based implementation of the Kalman filter (KF) for extremely high-dimensional, possibly nonlinear and non-Gaussian state estimation problems. Its ability to handle state dimensions in the order of millions has made the EnKF a popular algorithm in different geoscientific disciplines. Despite a similarly vital need for scalable algorithms in signal processing, e.g., to make sense of the ever increasing amount of sensor data, the EnKF is hardly discussed in our field. This self-contained review paper is aimed at signal processing researchers and provides all the knowledge to get started with the EnKF. The algorithm is derived in a KF framework, without the often encountered geoscientific terminology. Algorithmic challenges and required extensions of the EnKF are provided, as well as relations to sigma-point KF and particle filters. The relevant EnKF literature is summarized in an extensive survey and unique simulation examples, including popular benchmark problems, complement the theory with practical insights. The signal processing perspective highlights new directions of research and facilitates the exchange of potentially beneficial ideas, both for the EnKF and high-dimensional nonlinear and non-Gaussian filtering in general.

Tables

Table 1. Table 1: Averaged errors ε ¯ ¯ 𝜀 \bar{\varepsilon} for different EnKF.

$N$	$c$	$\rho$	$\overline{\epsilon }$
$1000$	$1$	no	$0.29$
$40$	$1$	no	$0.44$
$40$	$1.05$	no	$0.33$
$40$	$1$	yes	$0.29$
$40$	$1.02$	yes	$0.28$
$20$	$>1$	no	$>1$
$20$	$1.01$	yes	$0.3$
$10$	$1.05$	yes	$0.34$

Equations

$$x_{k+1}$$

$$y_k$$

$$x_{k+1}$$

$$y_k$$

$${\hat{x}}_{k+1|k}$$

$$P_{k+1|k}$$

$${\hat{y}}_{k|k-1}$$

$$S_k$$

$${\hat{x}}_{k|k}$$

$$=(I-K_{k}H){\hat{x}}_{k|k-1}+K_k{y}_k,$$

$$P_{k|k}$$

$$K_k$$

$${\bar{x}}_{k|k}$$

$${\bar{P}}_{k|k}$$

$${\bar{x}}_{k|k}$$

$${\bar{P}}_{k|k}$$

$${\tilde{X}}_{k|k}=X_{k|k}-{\bar{x}}_{k|k}{1}^T=X_{k|k}(I_N-\frac{1}{N}11^T)$$

$$X_{k+1|k}=F{X}_{k|k}+G{V}_k.$$

$$X_{k+1|k}=f(X_{k|k},V_k),$$

$$X_{k|k}=(I-{\bar{K}}_{k}H)X_{k|k-1}+{\bar{K}}_k{y}_k{1}^T.$$

$$X_{k|k}=(I-{\bar{K}}_{k}H)X_{k|k-1}+{\bar{K}}_k{y}_k{1}^T-{\bar{K}}_k{E}_k.$$

$$Y_{k|k-1}=H{X}_{k|k-1}+E_k$$

$$X_{k|k}=X_{k|k-1}+{\bar{K}}_k(y_k{1}^T-Y_{k|k-1}).$$

$$Y_{k|k-1}=h(X_{k|k-1},E_k),$$

$${\tilde{X}}_{k|k-1}$$

$${\tilde{Y}}_{k|k-1}$$

$${\bar{M}}_k$$

$${\bar{S}}_k$$

$${\bar{K}}_k{\tilde{Y}}_{k|k-1}{\tilde{Y}}_{k|k-1}^T={\tilde{X}}_{k|k-1}{\tilde{Y}}_{k|k-1}^T.$$

$$X_{k|k}=X_{k|k-1}+({\tilde{X}}_{k|k-1}{\tilde{Y}}_{k|k-1}^T)({\tilde{Y}}_{k|k-1}{\tilde{Y}}_{k|k-1}^T{)}^{-1}(y_k{1}^T-Y_{k|k-1})$$

$$x_{k|k}^{(i)}=x_{k|k-1}^{(i)}+\frac{{\sum }_{j=1}^N(x_{k|k-1}^{(j)}{)}^2}{{\sum }_{j=1}^N(x_{k|k-1}^{(j)}+e_k^{(j)}{)}^2}(y_k-x_{k|k-1}^{(i)}-e_{k|k-1}^{(i)}),$$

$${\tilde{Z}}_{k|k-1}=H{\tilde{X}}_{k|k-1}$$

$${\bar{M}}_k$$

$${\bar{S}}_k$$

$${\bar{S}}_k=\left[\begin{array}{cc}\frac{1}{\sqrt{N-1}}{\tilde{Z}}_{k|k-1}& R^{\frac{1}{2}}\end{array}\right]\left[\begin{array}{c}\frac{1}{\sqrt{N-1}}{\tilde{Z}}_{k|k-1}^T\\ R^{\frac{\mathrm{T}}{2}}\end{array}\right].$$

$$P_{k|k}$$

$$=P_{k|k-1}^{\frac{1}{2}}(I-P_{k|k-1}^{\frac{\mathrm{T}}{2}}{H}^T{S}_k^{-1}H{P}_{k|k-1}^{\frac{1}{2}})P_{k|k-1}^{\frac{\mathrm{T}}{2}}$$

Full text

This self-contained review paper is aimed at signal processing researchers and provides all the knowledge to get started with the EnKF. The algorithm is derived in a KF framework, without the often encountered geoscientific terminology. Algorithmic challenges and required extensions of the EnKF are provided, as well as relations to sigma-point KF and particle filters. The relevant EnKF literature is summarized in an extensive survey and unique simulation examples, including popular benchmark problems, complement the theory with practical insights. The signal processing perspective highlights new directions of research and facilitates the exchange of potentially beneficial ideas, both for the EnKF and high-dimensional nonlinear and non-Gaussian filtering in general.

1 Introduction

Numerical weather prediction [1] is an extremely high-dimensional geoscientific state estimation problem. The state $x$ comprises physical quantities (temperature, wind speed, air pressure, etc.) at many spatially distributed grid points, which often yields a state dimension $n$ in the order of millions. Consequently, the Kalman filter (KF) [2, 3] or its nonlinear extensions [4, 5] that require the storage and processing of $n\times n$ covariance matrices cannot be applied directly. It is well-known that the application of particle filters [6, 7] is not feasible either. In contrast, the ensemble Kalman filter [8, 9] (EnKF) was specifically developed as algorithm for high-dimensional $n$. The EnKF

• •

is a random-sampling implementation of the KF;

• •

reduces the computational complexity of the KF by propagating an ensemble of $N<n$ state realizations;

• •

can be applied to nonlinear state-space models without the need to compute Jacobian matrices;

• •

can be applied to continuous-time as well as discrete-time state transition functions;

• •

can be applied to non-Gaussian noise densities;

• •

is simple to implement;

• •

does not converge to the Bayesian filtering solution for $N\rightarrow\infty$ in general;

• •

often requires extra measures to work in practice.

Also in the field of stochastic signal processing (SP) and Bayesian state estimation, high-dimensional problems become more and more relevant. Examples include SLAM [10] where $x$ contains an increasing number of landmark positions, or extended target tracking [11, 12] where $x$ can contain many parameters to describe the shape of the target. Furthermore, scalable SP algorithms are required to make sense of the ever increasing amount of data from sensors in everyday devices. EnKF approaches hardly appear in the relevant SP journals, though. In contrast, vivid theoretical development is documented in geoscientific journals under the umbrella term data assimilation (DA) [1]. Hence, a relevant SP problem is being addressed with only little participation from the SP community. Conversely, much of the DA literature makes little reference to relevant SP contributions. It is our intention to bridge this interesting gap.

There is further overlap that motivates for a closer investigation of the EnKF. First, the basic EnKF [9] can be applied to nonlinear and non-Gaussian state-space models because it is entirely sampling based. In fact, the state evolution in geoscientific applications is typically governed by large nonlinear black box prediction models derived from partial differential equations. Furthermore, satellite measurements in weather applications are nonlinearly related to the states [1]. Hence, the EnKF has long been investigated as nonlinear filter. Second, the EnKF literature contains so called localization methods [13, 14] to systematically approach high-dimensional problems by only acting on a part of the state vector in each measurement update. These ideas can be directly transferred to sigma point filters [5]. Third, the EnKF offers several interesting opportunities to apply SP techniques, e.g., via the application of bootstrap or regularization methods in the EnKF gain computation.

The contributions of this paper aim at making the EnKF more accessible to SP researchers. We provide a concise derivation of the EnKF based on the KF. A literature review highlights important EnKF papers with their respective contributions, and facilitates an easier access to the extensive and rapidly developing DA literature on the EnKF. Moreover, we put the EnKF in context with popular SP algorithms such as sigma point filters [4, 5] and the particle filter [6, 7]. Our presentation forms a solid basis for further developments and the transfer of beneficial ideas and techniques between the fields of SP and DA.

The structure of the paper is as follows. After an extensive literature review in Sec. 2, the EnKF is developed from the KF in Sec. 3. Algorithmic properties and challenges of the EnKF and the available approaches to face them are discussed in Sec. 4 and 5, respectively. Relations to other filtering algorithms are discussed in Sec. 6. The theoretical considerations are followed by numerical simulations in Sec. 7 and some concluding remarks in Sec. 8.

2 Filtering and EnKF literature

The following literature review provides important landmarks for the EnKF novice.

State-space models and the filtering problem have been investigated since the 1960s. Early results include the Kalman filter (KF) [2] as algorithm for linear systems, and the Bayesian filtering equations [15] as theoretical solution for nonlinear and non-Gaussian systems. Because the latter approach cannot be implemented in general, approximate filtering algorithms are required. With a leap in computing capacity, the 1990s saw major developments. The sampling-based sigma point Kalman filters [4, 5] started to appear. Furthermore, particle filters [6, 7] were developed to approximately implement the Bayesian filtering equations through sequential importance sampling.

The first EnKF [8] was proposed in a geoscientific journal in 1994 and introduced the idea of propagating ensembles to mimic the KF. A systematic error that resulted in an underestimated uncertainty was later corrected by processing “perturbed measurements”. This randomization is well motivated in [9] but also used in [13].

Interestingly, [8] remains the most cited EnKF paper111With over 3000 citations between 1994 and 2016., followed by the overview article [16] and the monograph [17] by the same author. Other insightful overviews from a geoscientific perspective are [18, 19]. Many practical aspects of operational EnKF for weather prediction and re-analysis are described in [19, 20, 21]. Whereas the aforementioned papers were mostly published in geoscientific outlets, a special issue of the IEEE Control Systems Magazine appeared with review articles [22, 23, 24] and an EnKF case study [25]. Still, the above material was written by EnKF researchers with a geoscientific focus and in the application specific terminology. Furthermore, references to the recent SP literature and other nonlinear KF variants [5] are scarce.

Some attention has been devoted to the EnKF also beyond the geosciences. Convergence properties for $N\rightarrow\infty$ have been established using different theoretical analyses of the EnKF [26, 27, 28]. Statistical perspectives are provided in the thesis [29] and the review [30]. A recommended SP view that connects the EnKF with Bayesian filtering and particle methods, including convergence results for nonlinear systems, is [31]. Examples of the EnKF as tool for tomographic imaging and target tracking are described in [32] and [33], respectively. Brief introductory papers that connect the EnKF with more established SP algorithms include [34] and [35]. The latter also served as basis for this article.

The majority of EnKF advances is still documented in geoscientific publications. Notable contributions include deterministic EnKF that avoid the randomization of [9] and propagate an ensemble of deviations from the ensemble mean [36, 37, 38, 16]. Their common basis as square root EnKF and the relation to square root KF [3] is discussed in [39]. The computational advantages in high-dimensional EnKF with small ensembles ($N\ll n$) come at the price of adverse effects, including the risk of filter divergence. The often encountered underestimation of uncertainty can be counteracted with ensemble inflation [40]. A scheme with two EnKF in parallel that provide each other with gain matrices to reduce unwanted “inbreeding” has been suggested in [13]. The benefit of such a double EnKF is, however, debated [41, 38]. The low-rank approximation of covariance matrices can yield spurious correlations between supposedly uncorrelated state components and measurements. Localization techniques such as local measurement updates [13, 16, 42] or covariance tapering [14, 43] let the measurement only affect a part of the state vector. In other words, localization effectively reduces the dimension of each measurement update. Inflation and localization are essential components of operational EnKF [19]. A list of further advances in the geoscientific literature is provided in the appendix of [17].

An interesting development for SP researchers is the reconsideration of particle filters (PF) for high-dimensional geoscientific problems, with seemingly little reference to SP literature. An early example is [44]. The well-known challenges, mostly related to the problem of importance sampling in high dimensions, are reviewed in [45, 46]. Several recent approaches [47, 48, 49] were successfully tested on a popular EnKF benchmark problem [50] that is also investigated in the simulation examples of this paper. Combinations of the EnKF with the deterministic sampling of sigma point filters [5] are given in [51] and [52]. However, the benefit of the unscented transformation [53, 5] in [52] is debated in [54]. Ideas to combine the EnKF with Gaussian mixture approaches are given in [55, 56, 57].

3 A Signal Processing Introduction to the Ensemble Kalman Filter

The underlying framework of our filter presentation are discrete-time state-space models [3, 15]. The Kalman filter and many EnKF variants are built upon the linear model

[TABLE]

with the $n$-dimensional state $x$ and the $m$-dimensional measurement $y$. The initial state $x_{0}$, the process noise $v_{k}$, and the measurement noise $e_{k}$ are described by $\operatorname{E}(x_{0})=\hat{x}_{0}$, $\operatorname{E}(v_{k})=0$, $\operatorname{E}(e_{k})=0$, $\operatorname{cov}(x_{0})=P_{0}$, $\operatorname{cov}(v_{k})=Q$, and $\operatorname{cov}(e_{k})=R$. In the Gaussian case, these moments completely characterize the distributions of $x_{0}$, $v_{k}$, and $e_{k}$.

Nonlinear relations in the state evolution and measurement equations can be described by a more general model

[TABLE]

More general noise and initial state distributions can, for example, be characterized by probability density functions $p(x_{0})$, $p(v_{k})$, and $p(e_{k})$.

Both (1) and (2) can be time-varying but the time indices for functions and matrices are omitted for convenience.

3.1 A brief Kalman filter review

The KF is an optimal linear filter [3] for (1) that propagates state estimates $\hat{x}_{k|k}$ and covariance matrices $P_{k|k}$.

The KF time update or prediction is given by

[TABLE]

The above parameters can be used to predict the output of (1) and its uncertainty via

[TABLE]

The measurement update adjusts the prediction results according to

[TABLE]

with a gain matrix $K_{k}$. Here, (5b) resembles a deterministic observer and only requires all eigenvalues of $(I-K_{k}H)$ inside the unit circle to obtain a stable filter. The optimal $K_{k}$ in the minimum variance sense is given by

[TABLE]

where $M_{k}$ is the cross-covariance between the state and output predictions. Alternatives to the covariance update (5c) exist, but the shown Joseph form [3] will simplify the derivation of the EnKF. Furthermore, it is valid for all gain matrices $K_{k}$ beyond (6) and numerically well-behaved. It should be noted that it is numerically advisable to obtain $K_{k}$ by solving $K_{k}S_{k}=M_{k}$ rather than explicitly computing $S_{k}^{-1}$ [58].

3.2 The ensemble idea

The central idea of the EnKF is to propagate an ensemble of $N<n$ (often $N\ll n$) state realizations $\smash{\{x^{(i)}_{k}\}_{i=1}^{N}}$ instead of the $n$-dimensional estimate $\hat{x}_{k|k}$ and the $n\times n$ covariance $P_{k|k}$ of the KF. The ensemble is processed such that

[TABLE]

Reduced computational complexity is achieved because the explicit computation of $\bar{P}_{k|k}$ is avoided in the EnKF recursion. Of course, this reduction comes at the price of a low rank approximation in (7b) that entails some negative effects and requires extra measures.

For our development it is convenient to treat the ensemble as an $n\times N$ matrix $X_{k|k}$ with columns $x^{(i)}_{k}$. This allows for the compact notation of the ensemble mean and covariance

[TABLE]

where $\mathbb{1}=[1,\ldots,1]^{T}$ is an $N$-dimensional vector and

[TABLE]

is an ensemble of deviations from $\bar{x}_{k|k}$, sometimes called ensemble anomalies [17]. The matrix multiplication in (9) provides a compact way to write the anomalies, but is not the most efficient way to compute them.

3.3 The EnKF time update

The EnKF time update is referred to as forecast in the geoscientific literature. In analogy to (3), a prediction ensemble $X_{k+1|k}$ is computed that carries the information in $\hat{x}_{k+1|k}$ and $P_{k+1|k}$. An ensemble of $N$ independent process noise realizations $\smash{\{v^{(i)}_{k}\}_{i=1}^{N}}$ with zero mean and covariance $Q$, stored as matrix $V_{k}$, is used in

[TABLE]

An extension to nonlinear state transitions (2a) is given by

[TABLE]

where we generalized $f$ to act on the columns of its input matrices. Apparently, the EnKF time update amounts to a one-step-ahead simulation of $X_{k|k}$. Consequently, also continuous-time dynamics can be considered by, for example, numerically solving partial differential equations to obtain $X_{k+1|k}$. Also non-Gaussian initial state and process noise distributions with arbitrary densities $p(x_{0})$ and $p(v_{k})$ can be employed as long as they allow sampling. Perhaps because of this flexibility, the time update is often omitted in the EnKF literature [9, 13].

3.4 The EnKF measurement update

The EnKF measurement update is referred to as analysis in the geoscientific literature. A prediction or forecast ensemble $X_{k|k-1}$ is processed to obtain the filtering ensemble $X_{k|k}$ that encodes the KF mean and covariance. We assume that a gain $\bar{K}_{k}=K_{k}$ is given and postpone its computation to the next section.

With $\bar{K}_{k}$ available, the KF update (5b) can be applied to each ensemble member as follows [8]

[TABLE]

The resulting ensemble average (8a) approximates the correct $\hat{x}_{k|k}$. However, with $y_{k}\mathbb{1}^{T}$ known, the sample covariance (8b) gives only the first term of (5c) and therefore fails to resemble $P_{k|k}$. A solution [9] is to account for the missing term $\bar{K}_{k}R\bar{K}_{k}^{T}$ by adding artificial zero-mean measurement noise realizations $\smash{\{e^{(i)}_{k}\}_{i=1}^{N}}$ with covariance $R$, stored as matrix $E_{k}$, according to

[TABLE]

Then, $X_{k|k}$ correctly resembles $\hat{x}_{k|k}$ and $P_{k|k}$. The model (1) is implicit in (13) because the matrix $H$ appears. If we, in analogy to (4), define a predicted output ensemble

[TABLE]

that encodes $\hat{y}_{k|k-1}$ and $S_{k}$, we can reformulate (13) to an update that resembles (5a):

[TABLE]

In contrast to (13), the update (15) is entirely sampling based. As a consequence, we can extend the algorithm to nonlinear measurement models (2b) by replacing (14) with

[TABLE]

where we generalized $h$ to accept matrix inputs similar to (11).

In the EnKF literature, the prevailing view of inserting artificial noise is that perturbed measurements $y_{k}\mathbb{1}^{T}-E_{k}$ are processed. This might appear unusual from an SP perspective since it suggests that information is distorted before processing. The introduction of output ensembles $Y_{k|k-1}$, in contrast, yields a direct connection to (4) and highlights the similarities between (15) and (5a).

An interesting point [57] is that the measurement enters linearly in (13) and (15) and merely shifts the ensemble locations. This highlights the EnKF roots in the linear KF in which $P_{k|k}$ also remains unchanged by $y_{k}$.

3.5 The EnKF gain

The optimal gain (6) in the KF is computed from the covariance matrices of the predicted state and output. In the EnKF, the required $M_{k}$ and $S_{k}$ are not available but must be approximated from the prediction ensembles (10) or (11), and (14) or (16).

A straightforward way to compute the EnKF gain $\bar{K}_{k}$ is to first compute the deviations or anomalies

[TABLE]

The computations (17) are entirely sampling-based, which is useful for the nonlinear case but introduces extra sampling errors. An obvious improvement for additive measurement noise $e_{k}$ with covariance $R$ is given in Sec. 5.2, together with the square root EnKF that avoid the insertion of $E_{k}$ altogether.

Similar to the KF, the gain $\bar{K}_{k}$ should be obtained from the solution of a linear system of equations

[TABLE]

4 Some Properties and Challenges of the EnKF

After a brief review of convergence results and the computational complexity of the EnKF, we discuss adverse effects that can occur in EnKF with finite ensemble size $N$.

4.1 Asymptotic convergence results

In linear systems the EnKF mean and covariance (7) converge to the KF results (5) as $N\rightarrow\infty$. This result has been established from different theoretical perspectives [26, 27, 31, 28].

For nonlinear systems the convergence is not as tangible. An investigation of the EnKF as particle system is given in [31], with the outcome that the EnKF does not give the Bayesian filtering solution except for the linear Gaussian case. An illustration of this property is given in the example of Sec. 7.2.

4.2 Computational complexity

For the complexity analysis we assume that we are only interested in the filtering results and that $n>N>m$, that is, the number of measurements is less than the ensemble size and state dimension.

The KF propagates the $n$-dimensional mean vector $\hat{x}_{k|k}$ and the $n\times n$ covariance matrix $P_{k|k}$ with $n(n+1)/2$ unique entries. These storage requirements of $\mathcal{O}(n^{2}/2)$ dominate for large $n>m$. The EnKF requires the storage of only $nN$ values. The space required to store the Kalman gain and other intermediate results is similar for the KF and EnKF. A reduction via sequential processing of measurements, as explained in Sec. 5.1, is possible for both.

For large $n$ the computational bottleneck of the KF is the covariance time update (3b). Without considering any potential structure in $F$, slightly less than $\mathcal{O}(n^{3})$ floating point operations (flops) are required. Contemporary matrix multiplication routines [58] achieve a reduction to roughly $\mathcal{O}(n^{2.4})$. The EnKF time update requires the propagation of $N$ realizations. If each propagation costs $\mathcal{O}(n^{2})$ flops, then time update is achieved in $\mathcal{O}(n^{2}N)$ flops.

The computation of the KF gain requires $\mathcal{O}(n^{2}m)$ flops for the computation of $M_{k}$ and $S_{k}$. The solution of (6) for $K_{k}$ amounts to $\mathcal{O}(m^{3})$. The actual measurement update (5) adds further $\mathcal{O}(n^{2}m)$ flops. For large $n$, the total cost is $\mathcal{O}(n^{2}m)$. In contrast, the EnKF parameters $\bar{M}_{k}$ and $\bar{S}_{k}$ can be computed in $\mathcal{O}(nmN)$ flops which, again, dominates the total cost of the measurement update for large $n$. So, the EnKF reduces the flop count by a factor $\frac{N}{n}$.

4.3 Sampling and coupling effects for finite ensemble size

A serious issue in the EnKF is a commonly noted tendency to underestimate the state uncertainty when using $N<n$ ensemble members [13, 19, 18]. In other words, the EnKF becomes over-confident and is likely to diverge [3] for too small $N$. A number of causes and related effects can be noted.

First, an ensemble $X_{k|k-1}$ with too few members might not cover the relevant regions of the state-space well enough after the time update (10). The underestimated spread persists in the measurement update (13) or (15) and also $X_{k|k}$ shows too little spread.

Second, the ensemble can only transport limited information and provide a sampling covariance $\bar{P}_{k|k}$, (7b) or (8b), of at most rank $N-1$. Consequently, identically zero entries of $P_{k|k}$ are difficult to reproduce and unwanted spurious correlations show up in $\bar{P}_{k|k}$. An example would be an unreasonably large correlation between the temperature at two distant locations on the globe. Of course, these correlations also affect $\bar{M}_{k}$ and $\bar{S}_{k}$, and thus the EnKF gain $\bar{K}_{k}$ in (18). As a result, state components that are actually uncorrelated to $y_{k}$ are erroneously updated in (13) or (15). Again, this leads to a reduction in ensemble spread.

Third, the ensemble members are nonlinearly coupled because the gain (18) is computed from the ensemble. This “inbreeding” [13] increases with each measurement update. An interesting side effect is that the ensemble is not independent and Gaussian, even for linear Gaussian problems. To illustrate this, we combine (18) and (15) to obtain

[TABLE]

and consider a linear model (1) with $n=1$, $H=1$, and a zero-mean $X_{k|k-1}$. Then, one member of $X_{k|k}$ is given by

[TABLE]

which clearly shows the nonlinear dependencies that impede Gaussianity of $x_{k|k}^{(i)}$. Although similar conclusions hold for the general case, concise effects on the ensemble spread are difficult to analyze. Some special cases ($n=1$ and $n=m$, $H=I$, $R\propto I$) are investigated in [26] and shown to produce an underestimated $\bar{P}_{k|k}$.

Finally, the random sampling in the measurement update by inserting measurement noise in (14) or (16) adds to the EnKF error budget. The inherent sampling errors can be reduced by using the square root EnKF of Sec. 5.2.

Experiments suggest that there is a threshold for $N$ above which the EnKF works. A good example is given in [42]. Sec. 5 discusses methods such as inflation and localization that can reduce this minimum $N$.

5 Important Extensions to the EnKF

The previous section highlighted some of the challenges of the EnKF. Here, we summarize the important extensions that are often essential to achieve a working EnKF with only few ensemble members.

5.1 Sequential updates

For the KF it is algebraically equivalent to carry out $m$ measurement updates (5) with the scalar components of $y_{k}$ instead of a batch update with the $m$-dimensional $y_{k}$, if the measurement noise covariance $R$ is diagonal [3]. Although often treated as a side note only, this technique is very useful. It yields a more flexible algorithm with regard to the availability of measurements at each time step $k$ and reduces the computational complexity. After all, the Kalman gain (6) merely requires a scalar division for each component of $y_{k}$. An extension to block-diagonal $R$ is imminent.

Motivated by the large number of measurements in geoscientific problems, sequential updates have also been suggested for the EnKF [14]. Because of the randomness inherent to the EnKF, there is no algebraic equivalence between sequential and batch updates. Hence, the order in which measurements are processed has an effect on the filtering results.

Furthermore, an unusual alternative interpretation of sequential updates can be found in the EnKF literature. Namely, measurement updates are carried out “grid point by grid point” [13, 16, 42], that is, an iteration is carried out over state rather than measurement components. We will return to this aspect in Sec. 5.4.

5.2 Model knowledge in the EnKF and square-root filters

The sampling based derivation of the EnKF in equations (10) through (18) facilitates a compact presentation. However, the randomization through $E_{k}$ in (14) or (16) adds Monte Carlo sampling errors to the EnKF budget. This section discusses how these errors can be reduced for linear systems (1). Similar results for nonlinear systems with additive noise follow easily. The interpretation of ensembles as (rectangular) matrix square roots is a common theme in the following approaches. In (8b), for instance, $\tfrac{1}{\sqrt{N-1}}\widetilde{X}_{k|k}$ can be seen as an $n\times N$ square root of $\bar{P}_{k|k}$.

A first thing to note is that the cross covariance $M_{k}$ in the KF and its ensemble equivalent $\bar{M}_{k}$ should not be influenced by additive measurement noise $e_{k}$. Therefore, it is reasonable to replace $\widetilde{Y}_{k|k-1}$ with

[TABLE]

The Kalman gain $\bar{K}_{k}$ is then computed as in the KF (6). Alternatively, a matrix square-root $R^{\frac{1}{2}}$ with $R^{\frac{1}{2}}R^{\frac{\mathrm{T}}{2}}=R$ can be used to factorize

[TABLE]

A QR decomposition [58] of the right matrix then yields a triangular $m\times m$ square root of $\bar{S}_{k}$, and the computation of $\bar{K}_{k}$ simplifies to forward and backward substitution. Such ideas have their origin in sigma point KF variants [59].

The KF permits offline computation of the covariance matrices $P_{k|k}$ for all $k$ because they do not depend on the measurements. In an EnKF for a linear system (1) we can mimic this behavior by propagating zero-mean ensembles $\widetilde{X}_{k|k}$ that only carry the information of $P_{k|k}$. This is the central idea of different square root EnKF [39] which were suggested in [36, 37, 38]. The name stems from a relation to square root KF [3] which propagate $n\times n$ matrix square roots $\smash{P^{\frac{1}{2}}_{k|k}}$ with $\smash{P^{\frac{1}{2}}_{k|k}P^{\frac{\mathrm{T}}{2}}_{k|k}=P_{k|k}}$. Most importantly, the artificial measurement noise and the inherent sampling error can be avoided.

The following derivation [39] rewrites an alternative expression for (5c) using a square root $\smash{P^{\frac{1}{2}}_{k|k-1}}$ and its ensemble approximation $\tfrac{1}{N-1}\widetilde{X}_{k|k-1}$:

[TABLE]

which requires the left hand side to be positive definite. This property is easily established for the positive definite $\bar{S}_{k}$ of (21c) after realizing that the left hand side of (23b) is a Schur complement [58] of a positive definite matrix.

Finally, the $N\times N$ matrix $\Pi_{k}^{\frac{1}{2}}$ can be used to create a deviation ensemble

[TABLE]

that correctly encodes $P_{k|k}$ without using any random perturbations. Other variants update the deviation ensemble via a multiplication from the left [36], which is more costly for large $n$. Some more conditions on $\Pi_{k}^{\frac{1}{2}}$ must be met for $\widetilde{X}_{k|k}$ to remain zero-mean [60, 61].

The actual filtering is achieved by updating a single estimate according to

[TABLE]

where $\bar{K}_{k}$ is computed from the deviation ensembles.

There are indications that in nonlinear and non-Gaussian systems the sampling based EnKF variants should be preferable over their square root counterparts: A low-dimensional example is studied in [62]; the impression is confirmed for a high-dimensional problem in [63].

5.3 Ensemble inflation

Ensemble or covariance inflation is a measure to counteract the tendency of the EnKF to underestimate the state uncertainty for small $N$, and an important ingredient in operational EnKF [18]. The spread of the prediction ensemble $X_{k|k-1}$ is increased according to

[TABLE]

with a factor $c>1$. In the EnKF context, this heuristic has been proposed in [40]. Related concepts are dithering in the PF [7] and the “fudge factor” to increase $P_{k|k-1}$ in the KF [64]. Extensions to adaptive inflation, where $c$ is adjusted online, are discussed in [23].

5.4 Localization

Localization is a technique to address the issue of spurious correlations in the EnKF, and a crucial feature of operational EnKF [18, 19]. The underlying idea applies equally well to the EnKF and the KF, and can be used to systematically update only a part of the state vector with each measurement.

In order to explain the concept, we regard the KF measurement update for a linear system (1) with a low-dimensional222We assume that the components can be processed sequentially. measurement $y_{k}$.

Let $x=x_{k|k-1}$ and $P=P_{k|k-1}$ for notational convenience. It is possible to permute the state components such that

[TABLE]

Only the part $x_{1}$ appears in the measurement equation (1b) $y_{k}=H_{1}x_{1}+e_{k}$. While $x_{2}$ is correlated to $x_{1}$, there is zero correlation between $x_{1}$ and $x_{3}$. As a consequence, many submatrices of $P$ vanish in the computation of

[TABLE]

A KF measurement update (5) with the above $K_{k}$ does not affect the $x_{3}$ estimate or covariance. Hence, there is a lower-dimensional measurement update that only alters the statistics of $x_{1}$ and $x_{2}$.

Localization in the EnKF enforces the above structure using two prevailing techniques, local updates [13, 16, 42] and covariance tapering [14, 43]. Both rely on prior knowledge of the covariance structure. For example, the state components are often connected to geographic locations in geoscientific applications. From the underlying physics it is reasonable to assume zero correlation between distant states. Unfortunately, this viewpoint is not transferable to high-dimensional problems in general.

Local updates were introduced for the sampling based EnKF in [13] and for different square root EnKF in [16, 42]. Nonlinear measurement functions (2b) are linearized in the latter two. All of the above references update the state vector “grid point by grid point”, which appears unusual from a KF perspective [3]. In an iteration, local state vectors of small dimension ($<N$) are chosen and updated with a subset of supposedly relevant measurements. These “full rank” updates avoid some of the problems associated with small $N$ and large $n$. However, discontinuities between state components are introduced [65]. Some heuristics to combine the local ensembles and further implementation details can be found in [42, 66].

Under the assumption of the structure in (27), a local analysis would amount to an EnKF update of the $x_{1}$- and $x_{2}$-components only, to avoid errors in $x_{3}$.

Covariance tapering was introduced in [13]. It contradicts the EnKF idea in the sense that the ensemble covariance $\bar{P}_{k|k-1}$ of $X_{k|k-1}$ is processed. However, it will become clear that not all entries of $\bar{P}_{k|k-1}$ must be computed. Prior knowledge of a covariance structure as in (27) is used to create an $n\times n$ matrix $\rho$ with entries in $[0,1]$, and a tapered covariance $(\rho\circ\bar{P}_{k|k-1})$ is computed. Here, $\circ$ denotes the element-wise Hadamard or Schur product [58]. A typical $\rho$ has ones on the diagonal and decays smoothly to zero for unwanted off-diagonal elements. The standard choice uses a compactly supported correlation function from [67] and is discussed in [14, 43, 65]. Subsequently, the Kalman gain is computed as in the KF (6) using

[TABLE]

where we assumed a linear measurement relation (1b).

There are some technicalities associated with the tapering operation. Only positive semi-definite $\rho$ guarantee that $(\rho\circ\bar{P}_{k|k-1})$ is a valid covariance [26]. Full rank $\rho$ yield an increased rank in $(\rho\circ\bar{P}_{k|k-1})$ [14]. However, low rank $\rho$ do not necessarily decrease the rank of $(\rho\circ\bar{P}_{k|k-1})$. A closely related problem to finding valid (positive semi-definite or definite) $\rho$ is the creation of covariance functions and kernels in Gaussian processes [68]. Here, a methodology to create more complicated kernels from simpler ones could be used to create $\rho$.

Unfortunately, the Hadamard product cannot be formulated as an operation on the ensembles in general. Still, the computational requirements can be limited by only working with the non-zero elements of $(\rho\circ\bar{P}_{k|k-1})$. Furthermore, it is common to avoid the computation of $\bar{P}_{k|k-1}$ using

[TABLE]

instead of (29a) and to skip the tapering in $S_{k}$ altogether [43]. After all, for low-dimensional $y_{k}$ (small $m$) $\bar{M}_{k}$ has the strongest influence on the gain $\bar{K}_{k}$. Also the matrix $\rho_{M}$ is constructed from prior knowledge about the correlation. In the geoscientific context, where the state components and measurements are associated with geographic locations, this is easy. In general, however, it might not be possible to devise an appropriate $\rho_{M}$. Other variants [14, 26, 65] with tapering for $\bar{S}_{k}$ exist and have in common that they are only identical to (29) for $H=I$.

Some relations between local updates and covariance tapering are discussed in [65]. For the structure in (27) we can suggest a rank-1 taper $\rho$ that establishes a correspondence between the two concepts. Let $r_{1}$ and $r_{2}$ be vectors of the same dimensions as $x_{1}$ and $x_{2}$, respectively, that contain all ones. Let $r_{3}$ be a zero vector of the same dimension as $x_{3}$ and $r^{T}=[r_{1}^{T},r_{2}^{T},r_{3}^{T}]$. Then $\rho=rr^{T}$ removes all entries from $\bar{P}_{k|k-1}$ that would disappear in (28) anyhow. Furthermore, the Hadamard product for the rank-1 $\rho$ can be written as an operation on the ensemble $\widetilde{X}_{k|k-1}$ using

[TABLE]

The multiplication with $\operatorname{diag}(r)$ merely removes the rows corresponding to $x_{3}$, which establishes an equivalence between local updates and covariance tapering. By picking a smoothly decaying $r$ we can furthermore avoid the discontinuities associated with local updates.

5.5 The EnKF gain and least squares

A parallel to least squares problems can be disclosed by closer inspection of the equation (18) that is used to compute the EnKF gain $\bar{K}_{k}$. Perhaps more apparent in the transpose of (18), in

[TABLE]

that are to be solved for each row of $\bar{K}_{k}$ and $\widetilde{X}_{k|k-1}$. Hence, the EnKF iteration can be carried out without explicitly computing any sample covariance matrices if instead efficient solutions to the problem (32b) are employed. Furthermore, the problem (32b) could be modified using regularization [69] to enforce sparsity in $\bar{K}_{k}$. This would be an alternative approach to the localization methods discussed earlier. Related ideas to improve the Kalman gain using bootstrap methods [69] for computing $\bar{M}_{k}$ and $\bar{S}_{k}$ in (17) are discussed in [70, 71].

6 Relations to other algorithms

The EnKF for nonlinear systems (2) differs from other sampling based nonlinear filters such as sigma point KF [5] or particle filters (PF) [7]. One reason for this is that the EnKF approximates the KF algorithm (with the side effect that it can be applied to (2)) rather than trying to solve the nonlinear filtering problem directly.

The biggest difference between the EnKF and sigma point filters [5] such as the unscented KF [4, 53] or divided difference KF [59] is the measurement update. Whereas the EnKF updates its ensembles, the latter carry out the KF measurement update (5) using approximately computed mean values and covariance matrices. That is, the samples or sigma points are condensed into a filtering estimate $\hat{x}_{k|k}$ and its covariance $P_{k|k}$, which entails a loss of information and can be seen as an inherent Gaussian assumption on the filtering density $p(x_{k}|y_{1:k})$. In contrast, the EnKF can preserve more information and deviations from Gaussianity in the ensemble. Similarities appear in the gain computations of the EnKF and sigma point KF. In both, the Kalman gain appears as a function of the sampling covariance matrices, although with the deterministic sigma points and weights in the latter. With their origin in the KF, both sigma point filters and the EnKF can be expected to share difficulties with multimodal posterior distributions.

Similar to the EnKF, the PF propagates $N$ state realizations that are called particles. For the bootstrap particle filter [6], the prediction step corresponds to the EnKF time update (11). Apart from that, however, the differences dominate. First, the PF is designed as an approximate solution of the Bayesian filtering equations [15] using sequential importance sampling [7]. For $N\rightarrow\infty$, the PF solution recovers the true filtering density. Second, the samples in basic PF variants are generated from a proposal distribution only once every time instance and then left untouched. The measurement update amounts to updating the particle weights, which leads to a degeneracy problem for large $n$. In the EnKF, in contrast, the ensemble members are influenced by the time and the measurement update. Third, the PF relies on a crucial resampling step that is not present in the EnKF.

In summary, the EnKF appears as a distinct algorithm besides sigma point KF and PF. Its properties and potential for nonlinear problems remain to be fully investigated. Existing results that the EnKF does not converge to the Bayesian filtering recursion [31] remain to be interpreted in a constructive manner.

7 Instructive Simulation Examples

Four examples are discussed in greater detail, among them one popular benchmark problem of the SP and DA literature each.

7.1 A scalar linear Gaussian model

The first example illustrates the tendency of the EnKF to underestimate the state uncertainty. A related example is studied in [38]. We compare the EnKF variance $\bar{P}_{k|k}$ to the $P_{k|k}$ of the KF via Monte Carlo simulations on the simple scalar state-space model

[TABLE]

A trajectory of (33) is simulated and a KF is used to compute the optimal variances $P_{k|k}$. Because the model is time-invariant, the $P_{k|k}$ quickly converge to a constant value. For $k>3$ $P_{k|k}=0.0092$ is obtained.

Next, 10000 Monte Carlo experiments with a sampling based EnKF with $N=5$ are performed. The distribution of obtained $\bar{P}_{k|k}$ for $k=10$ is illustrated in Fig. 1. The vertical lines indicate the $P_{k|k}$ of the KF and the median and mean of the $\bar{P}_{k|k}$ outcomes.

The average $\bar{P}_{k|k}$ over the Monte Carlo realizations is close to the desired $P_{k|k}$. However, there is a large spread among the $\bar{P}_{k|k}$ and the distribution is skewed toward zero with its median below $P_{k|k}$. Although $N>n$, there is a tendency to underestimate $P_{k|k}$.

In order to clarify the reason for this behavior and whether it has to do with the coupling between the EnKF $\bar{K}_{k}$ and the ensemble members, we repeat the experiment with an EnKF that uses the gain of the stationary KF for all $k$. The resulting outcomes are illustrated in Fig. 2.

Now, the average $\bar{P}_{k|k}$ is correct. However, the median shows that there is still more probability mass below $P_{k|k}$. The tendency to underestimate $P_{k|k}$ and the remaining spread must be due to random sampling errors. For larger $N$ the effect vanishes, and the median and mean of $\bar{P}_{k|k}$ appear similar for $N\geq 10$.

7.2 The particle filter benchmark

In the second example we show that the EnKF does not converge to the Bayesian filtering solution in nonlinear systems as $N\rightarrow\infty$ [31]. A well-known benchmark problem from the PF literature [6] is used. The model is specified by

[TABLE]

with independent $v_{k}\sim\mathcal{N}(0,10)$, $e_{k}\sim\mathcal{N}(0,1)$, and $x_{0}\sim\mathcal{N}(0,1)$. Because the model is scalar, the Bayesian filtering densities $p(x_{k}\,|\,y_{1:k})$ can be computed numerically using point mass filters (PMF) [72]. A sampling based EnKF with $N=500$ is tested and kernel density estimates are used to obtain approximations of $p(x_{k}\,|\,y_{1:k})$ from the ensembles. For comparison, we include a closely related sigma point KF variant that uses Monte Carlo integration with $N=500$ samples [5]. The only difference to the EnKF is that this Monte Carlo KF (MCKF) carries out the KF measurement update (5) to propagate a mean and a variance. We illustrate the results as Gaussian densities.

Fig. 3 shows the prediction results for $k=150$. The PMF reference solution is bimodal with one mode close to the true state. The reason for this lies in the squared $x_{k}$ in (34b).

The EnKF prediction resembles the PMF well except for the random variations in the kernel density estimate. The MCKF cannot represent the multimodality but the Gaussian bell covers the relevant regions.

The filtering results for $k=150$ are shown in Fig. 4.

The PMF reference solution has much narrower peaks after including $y_{k}$. The EnKF provides a skewed density that does not resemble $p(x_{k}\,|\,y_{1:k})$ even though the EnKF prediction approximated $p(x_{k}\,|\,y_{1:k-1})$ well. This is the main take-away result and confirms [31]. Again, the MCKF exhibits a large variance.

Further filtering results for the PMF and EnKF are shown in Fig. 5.

It can be seen that the EnKF solutions sometimes resemble the PMF very well, but not always. Similar statements can be made for the prediction results. Dots in Fig. 5 illustrate the mean values as state estimates. Especially for the PMF, it can be seen that the mean (though optimal in a minimum variance sense [3]) is debatable for multimodal densities. Often, all estimates are quite close. Fig. 6 provides the estimation error densities obtained from $100$ Monte Carlo experiments with $151$ time steps each.

The PMF mean estimates exhibit a larger peak around [math]. The estimation errors for the EnKF and MCKF appear similar. This is surprising because the latter employs a Gaussian approximation at each time step. Both error densities have heavier tails than the PMF density. All estimation errors appear unbiased.

7.3 Batch smoothing using the EnKF

We here show how to use the EnKF as smoothing algorithm by estimating batches of states. This allows us to compare its performance for $N<n$ in problems of arbitrary dimension.

First, we formulate an “augmented state” that comprises an entire trajectory of $L+1$ steps,

[TABLE]

with dimension $n=(L+1)n_{x}$. Second, we note that the measurements $y_{k}$, $k=1,\dotsc,L$, have uncorrelated measurement noise and known relations to the components of $\xi$. For linear systems (1), the predicted mean and covariance of $\xi$ can be easily derived, and smoothed estimates of all $x_{k}$, $k=0,\dotsc,L$, can be obtained by sequentially processing all $y_{k}$ in KF measurement updates for $\xi$.

Also other smoothing variants and the Rauch-Tung-Striebel (RTS) algorithm can be derived from state augmentation approaches [3]. Due to its sequential nature, however, the RTS smoother does not provide joint covariance matrices of $x_{k}$ and $x_{k+i}$ for $i\neq 0$. Except for this and the higher computational complexity of working with $\xi$, the batch and RTS smoothers are equivalent for (1).

An EnKF approach to batch smoothing mimics the above. A prediction ensemble for $\xi$ is obtained by simulating $N$ trajectories for random process noise and initial state realizations. This can also be carried out for nonlinear models (2). Then, sequential EnKF measurement updates are performed for all $y_{k}$.

For our experiments we use a tracking problem with a constant velocity model [64] and position measurements. The low-dimensional state is given by

[TABLE]

With $n_{x}=4$ and $L=49$ we obtain $n=200$ as dimension of $\xi$. The RTS solution is compared to EnKF of ensemble size $N=\{10,20,50\}$. Monte Carlo errors are reduced using (21) in the gain computations.

A realization of a true trajectory and its measurements is provided in Fig. 7 together with the RTS estimate and an ensemble of $N=50$ trajectories.

The latter are the initial ensemble of an EnKF. The ensemble is well gathered around the initial position but fans out wildly. Fig. 8 shows the ensemble after an update with $y_{L}$ only.

The measurement at the end of the trajectory provides an anchor point and quickly reduces the spread of the ensemble. Fig. 9 shows the result after processing all measurements in sequential order from first to last.

The true trajectory and the RTS estimate are mostly covered well by the ensemble. The EnKF with $N=50$ appears consistent in this respect. Position errors for the RTS and the EnKF are provided in Fig. 10.

The EnKF performs slightly worse than the RTS but still gives good results for $N=50$, without extra inflation or localization.

The next experiment explores the EnKF for $N=10$. Fig. 11 shows the ensemble after processing all measurements.

The ensemble is compactly gathered but does not cover the true trajectory well. The EnKF is overconfident.

A last experiment explores how well an EnKF with $N=20$ captures the uncertainty of the state estimate. Furthermore, we discuss effects of the order in which the measurements are processed. Specifically, we compare the ensemble covariance of the positions $\mathsf{x}_{k}$ to the exact $\operatorname{cov}(\mathsf{x}_{k},\mathsf{x}_{i})$, $i,k=0,\dotsc,L$, obtained by KF updates for the augmented state $\xi$.

The exact covariance after processing all measurements is illustrated in Fig. 12.

Row $k$ in the matrix defines the covariance function between $\mathsf{x}_{k}$ and the remaining $\mathsf{x}$ positions. The banded structure indicates that subsequent positions are more related than, say, $\mathsf{x}_{0}$ and $\mathsf{x}_{L}$.

Fig. 13 shows the corresponding EnKF covariance after processing the measurements from $y_{1}$ to $y_{L}$.

The off-diagonal elements do not decay uniformly as in Fig. 12 and spurious positive and negative correlations appear. Furthermore, the correct temporal order of measurements entails an unwanted structure. Later $\mathsf{x}_{k}$ are rated more uncertain according to the lighter areas in the lower right corner of Fig. 13.

A covariance after processing the measurements in random order is shown in Fig. 14.

The spurious correlations persist but the diagonal elements appear more homogeneous.

From the above experiments we conclude that the EnKF can provide good estimates for ensembles with $N<n$. However, there is a minimum $N$ required to obtain consistent results without further measures such as localization or inflation. We have shown adverse effects such as ensembles with too little spread and spurious correlations.

7.4 The 40-dimensional Lorenz model

Our final example is a benchmark problem from the EnKF literature. We investigate the 40-dimensional Lorenz-96 model333Also known as the Lorenz-96, L95, L96, or L40 model. from [50] that is used in, e.g., [36, 38, 42, 66, 60, 47, 49]. The state $\mathsf{x}$ mimics an atmospheric quantity at equally spaced locations along a circle. Its evolution is specified by the nonlinear differential equation

[TABLE]

where $j=1,\dotsc,40$ indexes the components of $\mathsf{x}$, with the convention that $\mathsf{x}(0)=\mathsf{x}(40)$ etc. Instead of the commonly used forcing term $\mathsf{F}(j)=8$, we assume time-dependent $\mathsf{F}_{k}(j)\sim\mathcal{N}(8,1)$ that are constant for time intervals $\mathsf{T}=0.05$ only and act as process noise. A Runge-Kutta method (RK4) is used to discretize (37) to obtain the nonlinear state difference equation (2a) with $x_{k}=\mathsf{x}_{k}$ and $v_{k}=\mathsf{F_{k}}$. The step size $\mathsf{T}$ corresponds to about six hours if $\mathsf{x}$ were an atmospheric quantity on a latitude circle of the earth [50]. Although the model (37) is said to be chaotic, the effects are only mild for short integration times $\mathsf{T}$. In our experiments all $n=40$ states are measured with additive Gaussian noise $e_{k}\sim\mathcal{N}(0,I)$. The initial state is Gaussian with $x_{0}\sim\mathcal{N}(0,P_{0})$, where $P_{0}$ is drawn from a Wishart distribution with seed matrix $I_{n}$ and $n$ degrees of freedom. Fig. 15 illustrates how the state evolves over several time steps.

There is a tendency for peaks to move “westwards” as $k$ increases.

We note that there are also alternative approaches for estimating $\mathsf{x}$, for example, by first linearizing and then discretizing (37). However, we adopt the RK4 discretization of the EnKF literature that yields a state transition that is easy to evaluate but difficult to linearize. Because of this, the EKF [3] cannot be applied easily and we obtain a challenging benchmark problem.

We use sampling based EnKF to estimate long state sequences of $L=10^{4}$ time steps. Following [42, 38], the performance is assessed by the error

[TABLE]

where $\hat{x}_{k|k}$ is the ensemble mean. We use the average $\varepsilon_{k}$ for $k=100,\dotsc,L$, denoted by $\bar{\varepsilon}$, as quantitative performance measure for different EnKF. Useful EnKF must yield $\bar{\varepsilon}<1$, which is the error when simply taking $\hat{x}_{k|k}=y_{k}$.

First, we compute a reference solution using an EnKF with $N=1000$. Without any localization or inflation $\bar{\varepsilon}=0.29$ is achieved. Fig. 16 shows the sample covariance $\bar{P}_{k|k-1}$ of a prediction ensemble $X_{k|k-1}$, our best guess of the true covariance.

The banded structure reveals that the problem is suitable for localization. Hence, we construct a matrix $\rho$ for covariance tapering from a compactly supported correlation function [67] that is also used in [14, 43, 38, 26] and appears to be the standard choice. The chosen $\rho$ is a Toeplitz matrix because the components of $x_{k}$ are at equidistant locations, and shown in Fig. 17.

Next, EnKF with different ensemble sizes $N$, ensemble inflation factors $c$, with or without tapering, are compared. The obtained errors $\bar{\varepsilon}$ are summarized in Table 1.

For $N=n=40$, we obtain a worse $\bar{\varepsilon}$ than for $N=1000$. While inflation without tapering does reduce the error slightly, the covariance tapering even yields a better result that the EnKF with $N=1000$. Further improvements are obtained by combining inflation and tapering.

Fig. 18 shows the estimation error $x_{k}-\hat{x}_{k|k}$ for $k=10^{4}$, $N=40$, $c=1.02$, and tapering with $\rho$. In the background, the ensemble deviations $\widetilde{X}_{k|k}$ are illustrated. The estimation error is mostly contained in the intervals spanned by the ensemble, hence the EnKF is consistent. Tests on EnKF with $N=20$ reveal convergence problems, even with inflation the initial estimation error persists. With the help of tapering, however, a competitive error can be achieved. Even further reduction to $N=10$ is possible with tapering and inflation. The required inflation factor $c$ must be increased to counteract the lack of ensemble spread. Similar to Fig. 18, Fig. 19 illustrates the estimation error and deviation ensemble for $k=10^{4}$, $N=10$, $c=1.05$, and tapering with $\rho$. Although the obtained error is larger than for $N=40$, the ensemble deviations represent the estimation uncertainty well.

A number of lessons have been learned from related experiments. As alternative to the $\rho$ in Fig. 17, a simpler taper that contains only ones and zeros to enforce the banded structure was used. Although this $\rho$ was indefinite, a reduction in $\bar{\varepsilon}$ was achieved without any numerical issues. Hence, the specific structure of $\rho$ appears secondary. The smooth $\rho$ of Fig. 17 remains preferable in terms of $\bar{\varepsilon}$, though. Sequential processing of the measurements did not degrade the performance. Experiments without process noise give the lower errors $\bar{\varepsilon}$ from, e.g., [42, 38].

8 Concluding Remarks

With this paper, we have given a comprehensive and easy to understand introduction to the EnKF for signal processing researchers. The origin of the EnKF in the KF and its simple implementation have been demonstrated. The unique literature review provides quick access to the most relevant papers in the plethora of geoscientific EnKF publications. Furthermore, we have discussed the challenges related to small ensembles for high-dimensional states, $N<n$, and the available solutions such as localization or inflation. Finally, we have tested the EnKF on signal processing and EnKF benchmark problems.

With its scalability and simple implementation, even for nonlinear and non-Gaussian problems, the EnKF stands out as viable candidate for many state estimation problems. Furthermore, localization ideas and advanced concepts for estimating covariance matrices and the EnKF gain from the limited information in the ensembles provide new research directions for the EnKF and high-dimensional filters in general, hopefully with an increased participation from the signal processing community.

Acknowledgments

This work was supported by the project Scalable Kalman Filters granted by the Swedish Research Council (VR).