Revising the stochastic iterative ensemble smoother

TL;DR

This paper simplifies and improves the stochastic iterative ensemble smoother (EnRML) for large nonlinear inverse problems, addressing computational and noise issues while maintaining its accuracy.

Contribution

It provides a streamlined formulation of EnRML, clarifies its linearization and subspace properties, and enhances its computational efficiency without altering the results.

Findings

EnRML is made more computationally efficient and less noisy.

The relation between ensemble linearizations and average sensitivity is clarified.

Benchmark results show improved performance on the Lorenz-96 model.

Abstract

Ensemble randomized maximum likelihood (EnRML) is an iterative (stochastic) ensemble smoother, used for large and nonlinear inverse problems, such as history matching and data assimilation. Its current formulation is overly complicated and has issues with computational costs, noise, and covariance localization, even causing some practitioners to omit crucial prior information. This paper resolves these difficulties and streamlines the algorithm, without changing its output. These simplifications are achieved through the careful treatment of the linearizations and subspaces. For example, it is shown (a) how ensemble linearizations relate to average sensitivity, and (b) that the ensemble does not lose rank during updates. The paper also draws significantly on the theory of the (deterministic) iterative ensemble Kalman smoother (IEnKS). Comparative benchmarks are obtained with the Lorenz-96 model with these two smoothers and the ensemble smoother using multiple data assimilation (ES-MDA).

Equations

$$\mathbf{y}$$

$$p(\mathbf{x}|\mathbf{y})\propto p(\mathbf{y}|\mathbf{x})p(\mathbf{x}),$$

$$p(\mathbf{x})$$

$$=|2\pi {\mathbf{C}}_{\mathbf{x}}{|}^{-\frac{1}{2}}{e}^{-\frac{1}{2}\Vert \mathbf{x}-{\mathbf{\mu }}_{\mathbf{x}}{\Vert }_{{\mathbf{C}}_{\mathbf{x}}}^2}.$$

$$p(\mathbf{\delta })$$

$$p(\mathbf{y}|\mathbf{x})$$

$$J_{\mathbf{x},n}(\mathbf{x})$$

$$\nabla J_{\times }$$

$${\mathbf{C}}_{\times }^{-1}$$

$$\begin{array}{rl}{\mathbf{x}}_{n,i+1}& ={\mathbf{x}}_{\times }-{\mathbf{C}}_{\times }\nabla J_{\times }\\ & ={\mathbf{x}}_{\times }+{\mathbf{\Delta }}_{\times }^{\text{prior}}+{\mathbf{\Delta }}_{\times }^{\text{lklhd}},\end{array}$$

$${\mathbf{\Delta }}_{\times }^{\text{prior}}$$

$${\mathbf{\Delta }}_{\times }^{\text{lklhd}}$$

$${\mathbf{\Delta }}_{\times }^{\text{prior}}$$

$${\mathbf{\Delta }}_{\times }^{\text{lklhd}}$$

$${\mathbf{K}}_{\times }$$

$${\mathbf{C}}_{\mathbf{y}}$$

$$\mathbf{E}$$

$$\mathbf{D}$$

$$\mathbf{\Pi }\mathbf{\Pi }=\mathbf{\Pi }={\mathbf{\Pi }}^{\mathsf{T}}.$$

$${\mathbf{\Pi }}_{\mathbf{A}}\mathbf{A}$$

$${\mathbf{\Pi }}_{\mathbf{A}}=\mathbf{A}{\mathbf{A}}^{+}=({\mathbf{A}}^{\mathsf{T}}{)}^{+}({\mathbf{A}}^{\mathsf{T}}).$$

$$\mathbf{X}$$

$${\text{accentset}\mathbf{C}}_{\mathbf{x}}$$

$${\text{accentset}\mathbf{M}}_i$$

$$\underset{N\to \infty }{\lim }\text{accentset}\mathbf{M}$$

$${\text{accentset}\mathbf{K}}_i$$

$$\displaystyle={\mathbf{{X}}}{\mathbf{{Y}}}_{i}^{{\mathsf{T}}}\big{(}{\mathbf{{Y}}}_{i}{\mathbf{{Y}}}_{i}^{{\mathsf{T}}}+({N{-}1}){\mathbf{{C}}}_{\bm{\delta}}\big{)}^{-1}\,,$$

$${\mathbf{Y}}_i$$

$${\text{accentset}\mathbf{K}}_i$$

$${\text{accentset}\mathbf{C}}_{\mathbf{w},i}$$

$$\text{accentset}{\mathbf{\Delta }}_{\times }^{\text{lklhd}}$$

$$\mathbf{x}(\mathbf{w})=\text{accentset}\mathbf{x}+\mathbf{X}\mathbf{w},$$

$$\text{accentset}{\mathbf{\Delta }}_{\times }^{\text{prior}}$$

Full text

1 Introduction

Ensemble (Kalman) smoothers are approximate methods used for data assimilation (state estimation in geoscience), history matching (parameter estimation for petroleum reservoirs), and other inverse problems constrained by partial differential equations. Iterative forms of these smoothers, derived from optimization perspectives, have proven useful in improving the estimation accuracy when the forward operator is nonlinear. Ensemble randomized maximum likelihood (EnRML) is one such method.

This paper rectifies several conceptual and computational complications with EnRML, detailed in section 1.1. As emphasized in section 1.2, these improvements are largely inspired by the theory of the iterative ensemble Kalman smoother (IEnKS). Readers unfamiliar with EnRML may jump to the beginning of the derivation, starting in section 2, which defines the inverse problem and the idea of the randomized maximum likelihood method. Section 3 derives the new formulation of EnRML, which is summarized by LABEL:algo:GN_EnRML of section 3.7. Section 4 shows benchmark experiments obtained with various iterative ensemble smoothers. Appendix A provides proofs of some of the mathematical results used in the text.

1.1 Ensemble randomized maximum likelihood (EnRML): obstacles

The Gauss-Newton variant of EnRML was given by Gu and Oliver [2007]; Chen and Oliver [2012], with an important precursor from Reynolds et al. [2006]. This version explicitly requires the ensemble-estimated “model sensitivity” matrix, herein denoted ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$. As detailed in section 3, this is problematic because ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ is noisy and requires the computation of the pseudo-inverse of the “anomalies”, ${\mathbf{{X}}}_{i}^{+}$, for each iteration, $i$.

A Levenberg-Marquardt variant was proposed in the landmark paper of Chen and Oliver [2013b]. Its main originality is a partial resolution to the above issue by modifying the Hessian (beyond the standard trust-region step regularization): the prior ensemble covariance matrix is replaced by the posterior covariance (of iteration $i$): ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}\leftarrow{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{x}},i}$. Now the Kalman gain form of the likelihood increment is “vastly simplified”, because the linearization ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ only appears in the product {\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{x}},i}{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}^{\mathchoice{\raisebox{-1.50694pt}{\displaystyle{\mathsf{T}}}}{\raisebox{-1.50694pt}{\textstyle{\mathsf{T}}}}{\raisebox{-1.50694pt}{\scriptstyle{\mathsf{T}}}}{\raisebox{-1.50694pt}{\scriptscriptstyle{\mathsf{T}}}}}, which does not require ${\mathbf{{X}}}_{i}^{+}$. For the prior increment, on the other hand, the modification breaks its Kalman gain form. Meanwhile, the precision matrix form, i.e. their equation (10), is already invalid because it requires the inverse of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{x}},i}$. Still, in their equation (15), the prior increment is formulated with an inversion in ensemble space, and also unburdened of the explicit computation of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$. Intermediate explanations are lacking, but could be construed to involve approximate inversions. Another issue is that the pseudo-inverse of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}$ is now required (via ${\mathbf{{X}}}$), and covariance localization is further complicated.

An approximate version was therefore also proposed in which the prior mismatch term is omitted from the update formula altogether. This is not principled, and severely aggravates the chance of overfitting and poor prediction skill. Therefore, unless the prior mismatch term is relatively insignificant, overfitting must be prevented by limiting the number of steps or by clever stopping criteria. Nevertheless, this version has received significant attention in history matching.

This paper revises EnRML; without any of the above tricks, we formulate the algorithm such that there is no explicit computation of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$, and show how the product ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$ may be computed without any pseudo-inversions of the matrix of anomalies. Consequently, the algorithm is simplified, computationally and conceptually, and there is no longer any reason to omit the prior increment. Moreover, the Levenberg-Marquardt variant is a trivial modification of the Gauss-Newton variant. The above is achieved by improvements to the derivation, notably by

(a) improving the understanding of the sensitivity (i.e. linearizations) involved,

(b) explicitly and rigorously treating issues of rank deficiency and subspaces, and

(c) avoiding premature insertion of singular value decompositions (SVD).

1.2 Iterative ensemble Kalman smoother (IEnKS)

The contributions of this paper (listed by the previous paragraph) are original, but draw heavily on the theory of the IEnKS of Sakov et al. [2012]; Bocquet and Sakov [2012, 2014]. Relevant precursors include [Zupanski, 2005], as well as the iterative, extended Kalman filter [e.g. Jazwinski, 1970].

It is informally known that EnRML can be seen as a stochastic flavour of the IEnKS [Sakov et al., 2012]. Indeed, while the IEnKS update takes the form of a deterministic, “square-root” transformation, based on a single objective function, EnRML uses stochastic “perturbed observations”, associated with an ensemble of randomized objective functions.

Another notable difference is that the IEnKS was developed in the atmospheric literature, while EnRML was developed in the literature on subsurface flow. Thus, typically, the IEnKS is applied to (sequential) state estimation problems such as filtering for chaotic dynamical systems, while EnRML is applied to (batch) parameter estimation problems, such as nonlinear inversion for physical constants and boundary conditions. For these problems, EnRML is sometimes referred to as the iterative ensemble smoother (IES). As shown by Gu and Oliver [2007], however, EnRML is easily reformulated for the sequential problem. Vice-versa, the IEnKS may be formulated for the batch problem.

The improvements to the EnRML algorithm herein render it very similar to the IEnKS, also in computational cost. It thus fully establishes that EnRML is the stochastic ”counterpart” to the IEnKS. In spite of the similarities, the theoretical insights and comparative experiments of this paper should make it interesting also for readers already familiar with the IEnKS.

2 RML

Randomized maximum likelihood (RML) [Kitanidis, 1995; Oliver, 1996; Oliver et al., 2008] is an approximate solution approach to a class of inverse problems. The form of RML described here is a simplification, common for large inverse problems, without the use of a correction step (such as Metropolis-Hastings). This restricts the class of problems for which it is unbiased, but makes it more tractable [Oliver, 2017]. Similar methods were proposed and studied by Bardsley et al. [2014]; Liu et al. [2017]; Morzfeld et al. [2018].

2.1 The inverse problem

Consider the problem of estimating the unknown, high-dimensional state (or parameter) vector ${\bm{x}}\in\mathbb{R}^{M}$, given the observation ${\bm{y}}\in\mathbb{R}^{P}$. It is assumed that the (generic and typically nonlinear) forward observation process may be approximated by a computational model, $\mathcal{M}$, so that

[TABLE]

where the error, ${\bm{\delta}}$, is random and gives rise to a likelihood, $\mathop{}\!p({\bm{y}}|{\bm{x}})$.

In the Bayesian paradigm, prior information is quantified as a probability density function (pdf) called the prior, denoted $\mathop{}\!p({\bm{x}})$, and the truth, ${\bm{x}}$, is considered a draw thereof. The inverse problem then consists of computing and representing the posterior which, in principle, is given by pointwise multiplication:

[TABLE]

quantifying the updated estimation of ${\bm{x}}$. Due to the noted high dimensionality and nonlinearity, this can be challenging, necessitating approximate solutions.

The prior is assumed Gaussian, with mean ${\bm{\mu}}_{\bm{x}}$ and covariance ${\mathbf{{C}}}_{\bm{x}}$, i.e.

[TABLE]

For now, the prior covariance matrix, ${\mathbf{{C}}}_{\bm{x}}\in\mathbb{R}^{M\times M}$, is assumed invertible such that the corresponding norm, $\|{\bm{x}}\|_{{\mathbf{{C}}}_{\bm{x}}}^{2}={\bm{x}}^{{\mathsf{T}}}{\mathbf{{C}}}_{\bm{x}}^{-1}{\bm{x}}$, is defined. Note that vectors are taken to have column orientation, and that ${\bm{x}}^{{\mathsf{T}}}$ denotes the transpose.

The observation error, ${\bm{\delta}}$, is assumed drawn from:

[TABLE]

whose covariance, ${\mathbf{{C}}}_{\bm{\delta}}\in\mathbb{R}^{P\times P}$, will always be assumed invertible. Then, assuming ${\bm{\delta}}$ and ${\bm{x}}$ are independent and recalling equation 1,

[TABLE]

2.2 Randomize, then optimize

The Monte-Carlo approach offers a convenient representation of distributions as samples. Here, the prior is represented by the “prior ensemble”, $\{{\bm{x}}_{n}\}_{n=1}^{N}$, whose members (sample points) are assumed independently drawn from it. RML is an efficient method to approximately “condition” (i.e. implement 2 on) the prior ensemble, using optimization. Firstly, an ensemble of perturbed observations, $\{{\bm{y}}_{n}\}_{n=1}^{N}$, is generated as ${\bm{y}}_{n}={\bm{y}}+{\bm{\delta}}_{n}$, where ${\bm{\delta}}_{n}$ is independently drawn according to equation 4.

Then, the $n$-th “randomized log-posterior”, $J_{{\bm{x}},n}$, is defined by Bayes’ rule 2, except with the prior mean and the observation replaced by the $n$-th members of the prior and observation ensembles:

[TABLE]

The two terms are referred to as the model mismatch (log-prior) and data mismatch (log-likelihood), respectively.

Finally, these log-posteriors are minimized. Using the Gauss-Newton iterative scheme (for example) requires 7a its gradient and 7b its Hessian approximated by first-order model expansions, both evaluated at the current iterate, labelled ${\bm{x}}_{n,i}$ for each member $n$ and iteration $i$. To simplify the notation, define ${\bm{x}}_{\times}={\bm{x}}_{n,i}$. Objects evaluated at ${\bm{x}}_{\times}$ are similarly denoted; for instance, ${\mathbf{{M}}}_{\times}=\mathcal{M}^{\prime}({\bm{x}}_{\times})\in\mathbb{R}^{P\times M}$ denotes the Jacobian of $\mathcal{M}$ evaluated at ${\bm{x}}_{\times}$, and

[TABLE]

Application of the Gauss-Newton scheme yields:

[TABLE]

where the prior (or model) and likelihood (or data) increments are respectively given by:

[TABLE]

which can be called the “precision matrix” form.

Alternatively, by corollaries of the well known Woodbury matrix identity, the increments can be written in the “Kalman gain” form:

[TABLE]

where ${\mathbf{{I}}}_{M}\in\mathbb{R}^{M\times M}$ is the identity matrix, and ${\mathbf{{K}}}_{\times}\in\mathbb{R}^{M\times P}$ is the gain matrix:

[TABLE]

with

[TABLE]

As the subscript suggests, ${\mathbf{{C}}}_{{\bm{y}}}$ may be identified (in the linear case) as the prior covariance of the observation, ${\bm{y}}$, of equation 1; it is also the covariance of the innovation, ${\bm{y}}-\mathcal{M}({\bm{\mu}}_{\bm{x}})$. Note that if $P\ll M$, then the inversion of ${\mathbf{{C}}}_{{\bm{y}}}\in\mathbb{R}^{P\times P}$ for the Kalman gain form 10 is significantly cheaper than the inversion of ${\mathbf{{C}}}_{\times}\in\mathbb{R}^{M\times M}$ for the precision matrix form 9.

3 EnRML

Ensemble-RML (EnRML) is an approximation of RML in which the ensemble is used in its own update, by estimating ${\mathbf{{C}}}_{\bm{x}}$ and ${\mathbf{{M}}}_{\times}$. This section derives EnRML, and gradually introduces the new improvements.

Computationally, compared to RML, EnRML offers the simultaneous benefits of working with low-rank representations of covariances, and not requiring a tangent-linear (or adjoint) model. Both advantages will be further exploited in the new formulation of EnRML.

Concerning their sampling properties, a few points can be made. Firstly (due to the ensemble covariance), EnRML is biased for finite $N$, even for a linear-Gaussian problem, for which RML will sample the posterior correctly. This bias arises for the same reasons as in the ensemble Kalman filter [EnKF, van Leeuwen, 1999; Sacher and Bartello, 2008]. Secondly (due to the ensemble linearization), EnRML effectively smoothes the likelihood. It is therefore less prone to getting trapped in local maxima of the posterior [Chen and Oliver, 2012]. Sakov et al. [2018] explain this by drawing an analogy to the secant method, as compared to the Newton method. Hence, it may reasonably be expected that EnRML yields constructive results if the probability mass of the exact posterior is concentrated around its global maximum. Although this regularity condition is rather vague, it would require that the model be “not too nonlinear” in this neighbourhood. Conversely, EnRML is wholly inept at reflecting multimodality introduced through the likelihood, and so RML may be better suited when local modes feature prominently, as is quite common in problems of subsurface flow [Oliver and Chen, 2011]. However, while RML has the ability to sample multiple modes, it is difficult to predict to what extent their relative proportions will be accurate (without the costly use of a correction step such as Metropolis-Hastings). Further comparison of the sampling properties of RML and EnRML was done by Evensen [2018].

3.1 Ensemble preliminaries

For convenience, define the concatenations:

[TABLE]

which are known as the “ensemble matrix” and the “perturbation matrix”, respectively.

Projections sometimes appear through the use of linear regression. We therefore recall [Trefethen and Bau, 1997] that a (square) matrix ${\mathbf{{\Pi}}}$ is an orthogonal projector if

[TABLE]

For any matrix ${\mathbf{{A}}}$, let ${\mathbf{{\Pi}}}_{{\mathbf{{A}}}}$ denote the projector whose image is the column space of ${\mathbf{{A}}}$, implying that

[TABLE]

Equivalently, ${\mathbf{{\Pi}}}_{{\mathbf{{A}}}}^{\perp}{\mathbf{{A}}}={\mathbf{{0}}}$, where ${\mathbf{{\Pi}}}^{\perp}_{{\mathbf{{A}}}}={\mathbf{{I}}}-{\mathbf{{\Pi}}}_{{\mathbf{{A}}}}$ is called the complementary projector. The (Moore-Penrose) pseudo-inverse, ${\mathbf{{A}}}^{+}$, may be used to express the projector:

[TABLE]

Here, the second equality follows from the first by equation 15 and $({\mathbf{{A}}}^{+})^{{\mathsf{T}}}=({\mathbf{{A}}}^{{\mathsf{T}}})^{+}$. The formulae simplify further in terms of the SVD of ${\mathbf{{A}}}$.

Now, denote $\mathds{1}\in\mathbb{R}^{N}$ the (column) vector of ones. The matrix of anomalies, ${\mathbf{{X}}}\in\mathbb{R}^{M\times N}$, is defined and computed by subtracting the ensemble mean, ${\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}={\mathbf{{E}}}\mathds{1}/N$, from each column of ${\mathbf{{E}}}$. It should be appreciated that this amounts to the projection:

[TABLE]

where ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}={\mathbf{{I}}}_{N}-{{\mathbf{{\Pi}}}_{\mathds{1}}}$, with ${{\mathbf{{\Pi}}}_{\mathds{1}}}=\mathds{1}\mathds{1}^{{\mathsf{T}}}/N$.

Definition 1** (The ensemble subspace).**

*The flat (i.e. affine subspace) given by: $\{{\bm{x}}\in\mathbb{R}^{M}\;:\;[{\bm{x}}-{\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}]\in\operatorname{col}({\mathbf{{X}}})\}$. *

Similarly to section 2, iteration index ($i>0$) subscripting on ${\mathbf{{E}}}$, ${\mathbf{{X}}}$, and other objects, is used to indicate that they are conditional (i.e. posterior). The iterations are initialized with the prior ensemble: ${\bm{x}}_{n,0}={\bm{x}}_{n}$.

3.2 The constituent estimates

The ensemble estimates of ${\mathbf{{C}}}_{\bm{x}}$ and ${\mathbf{{M}}}_{\times}$ are the building blocks of the EnRML algorithm. The canonical estimators are used, namely the sample covariance 19a, and the least-squares linear regression coefficients 19b. They are denoted with the overhead bar:

[TABLE]

The anomalies at iteration $i$ are again given by ${\mathbf{{X}}}_{i}={\mathbf{{E}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$, usually computed by subtraction of ${\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}_{i}$. The matrix $\mathcal{M}({\mathbf{{E}}}_{i})$ is defined by the column-wise application of $\mathcal{M}$ to the ensemble members. Conventionally, $\mathcal{M}({\mathbf{{E}}}_{i})$ would also be centred in equation 19b, i.e. multiplied on the right by ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$. However, this operation (and notational burden) can be neglected, because ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}{\mathbf{{X}}}_{i}^{+}={\mathbf{{X}}}_{i}^{+}$, which follows from ${\mathbf{{\Pi}}}({\mathbf{{A}}}{\mathbf{{\Pi}}})^{+}=({\mathbf{{A}}}{\mathbf{{\Pi}}})^{+}$ [valid for any matrix ${\mathbf{{A}}}$ and projector ${\mathbf{{\Pi}}}$, as shown by Maciejewski and Klein, 1985].

Note that the linearization (previously ${\mathbf{{M}}}_{\times}$, now ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$) no longer depends on the ensemble index, $n$. Indeed, it has been called “average sensitivity” since the work of Zafari and Reynolds [2005]; Reynolds et al. [2006]; Gu and Oliver [2007]. However, this intuition has not been rigorously justified.111The formula 19b for ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ is sometimes arrived at via a truncated Taylor expansion of $\mathcal{M}$ around ${\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}_{i}$. This is already an approximation, and still requires further, indeterminate approximations to obtain any other interpretation than $\mathcal{M}^{\prime}({\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}_{i})$: the Jacobian evaluated at the ensemble mean.

This is accomplished by the following theorem.

Theorem 1** (Regression coefficients versus derivatives).**

Suppose the ensemble is drawn from a Gaussian. Then

[TABLE]

with “almost sure” convergence, and expectation $(\mathop{}\!\mathbb{E})$ in ${\bm{x}}$, which has the same distribution as the ensemble members. Regularity conditions and proof in appendix A.

A corollary of Theorem 1 is that ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}\approx\frac{1}{N}\sum_{n=1}^{N}\mathcal{M}^{\prime}({\bm{x}}_{n})$, justifying the “average sensitivity/derivative/gradient” description. The theorem applies for the ensemble of any Gaussian, and hence also holds for ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$. On the other hand, the generality of Theorem 1 is restricted by the Gaussianity assumption. Thus, for generality and precision, ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ should simply be labelled “the least-squares (linear) fit” of $\mathcal{M}$, based on ${\mathbf{{E}}}_{i}$.

Note that the computation 19b of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ seemingly requires calculating a new pseudo-inverse, ${\mathbf{{X}}}_{i}^{+}$, at each iteration, $i$; this is addressed in section 3.6.

The prior covariance estimate (previously ${\mathbf{{C}}}_{\bm{x}}$, now ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}$) is not assumed invertible, in contrast to section 2. It is then not possible to employ the precision matrix forms 9 because ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}^{-1}$ is not defined. Using the ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}^{+}$ in its stead is flawed and damaging because it is zero in the directions orthogonal to the ensemble subspace, so that its use would imply that the prior is assumed infinitely uncertain (i.e. flat) as opposed to infinitely certain (like a delta function) in those directions. Instead, one should employ ensemble subspace formulae, or equivalently (as shown in the following, using corollaries of the Woodbury identity), the Kalman gain form.

3.3 Estimating the Kalman gain

The ensemble estimates 19 are now substituted into the Kalman gain form of the update, equation 10 to 12. The ensemble estimate of the gain matrix, denoted ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{K}}}}_{i}$, thus becomes:

[TABLE]

where ${\mathbf{{Y}}}_{i}\in\mathbb{R}^{P\times N}$ has been defined as the prior (i.e. unconditioned) anomalies, under the action of the $i$-th iterate linearization:

[TABLE]

A Woodbury corollary can be used to express ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{K}}}}_{i}$ as:

[TABLE]

with

[TABLE]

The reason for labelling this matrix with the subscript ${\bm{w}}$ is revealed later. For now, note that, in the common case of $N\ll P$, the inversion in equation 24 is significantly cheaper than the inversion in equation 21. Another computational benefit is that ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i}$ is non-dimensional, improving the conditioning of the optimization problem [Lorenc, 1997].

In conclusion, the likelihood increment 10b is now estimated as:

[TABLE]

This is efficient because ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ does not explicitly appear in ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{K}}}}_{i}$ (neither in formula 21 nor 23), even though it is implicitly present through ${\mathbf{{Y}}}_{i}$ 22, where it multiplies ${\mathbf{{X}}}$. This absence

(a) is reassuring, as the product ${\mathbf{{Y}}}_{i}$ constitutes a less noisy estimate than just ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ alone [Chen and Oliver, 2012; Emerick and Reynolds, 2013b, figures 2 and 27, resp.];

(b) constitutes a computational advantage, as will be shown in section 3.6;

(c) enables leaving the type of linearization made for $\mathcal{M}$ unspecified, as is usually the case in EnKF literature.

3.4 Estimating the prior increment

In contrast to the likelihood increment 10b, the Kalman gain form of the prior increment 10a explicitly contains the sensitivity matrix, ${\mathbf{{M}}}_{\times}$. This issue was resolved by Bocquet and Sakov [2012] in their refinement of Sakov et al. [2012] by employing the change of variables:

[TABLE]

where ${\bm{w}}\in\mathbb{R}^{N}$ is called the ensemble “controls” [Bannister, 2017], also known as the ensemble “weights” [Ott et al., 2004], or “coefficients” [Bocquet and Sakov, 2013].

Denote ${\bm{w}}_{\times}$ an ensemble coefficient vector such that ${\bm{x}}({\bm{w}}_{\times})={\bm{x}}_{\times}$, and note that ${\bm{x}}({\bm{e}}_{n})={\bm{x}}_{n}$, where ${\bm{e}}_{n}$ is the $n$-th column of the identity matrix. Thus, $[{\bm{x}}_{n}-{\bm{x}}_{\times}]={\mathbf{{X}}}[{\bm{e}}_{n}-{\bm{w}}_{\times}]$, and the prior increment 10a with the ensemble estimates becomes:

[TABLE]

where there is no explicit ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$, which only appears implicitly through ${\mathbf{{Y}}}_{i}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$, as defined in equation 22 Alternatively, applying the subspace formula 23 and using ${\mathbf{{I}}}_{N}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i}({\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i})^{-1}$ yields:

[TABLE]

3.5 Justifying the change of variables

Lemma 1** (Closure).**

Suppose ${\mathbf{{E}}}_{i}$ is generated by EnRML. Then, each member (column) of ${\mathbf{{E}}}_{i}$ is in the (prior) ensemble subspace. Moreover, $\operatorname{col}({\mathbf{{X}}}_{i})\subseteq\operatorname{col}({\mathbf{{X}}})$.

Lemma 1 may be proven by noting that ${\mathbf{{X}}}$ is the leftmost factor in ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{K}}}}_{i}$, and using induction on equations 10a and 10b. Alternatively, it can be deduced [Raanes et al., 2019] as a consequence of the implicit assumption on the prior that ${\bm{x}}\sim\mathop{}\!\mathcal{N}({\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}},{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}})$. A stronger result, namely $\operatorname{col}({\mathbf{{X}}}_{i})=\operatorname{col}({\mathbf{{X}}})$, is conjectured in appendix A, but Lemma 1 is sufficient for the present purposes: it implies that there exists ${\bm{w}}_{\times}\in\mathbb{R}^{N}$ such that ${\bm{x}}({\bm{w}}_{\times})={\bm{x}}_{\times}$ for any ensemble member and any iteration. Thus, the lemma justifies the change of variables 26.

Moreover, using the ensemble coefficient vector (${\bm{w}}$) is theoretically advantageous as it inherently embodies the restriction to the ensemble subspace. A practical advantage is that ${\bm{w}}$ is relatively low-dimensional compared to ${\bm{x}}$, which lowers storage and accessing expenses.

3.6 Simplifying the regression

Recall the definition of equation 22: ${\mathbf{{Y}}}_{i}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$. Avoiding the explicit computation of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ used in this product between the iteration-$i$ estimate ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ and the initial (prior) ${\mathbf{{X}}}$ was the motivation behind the modification ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}\leftarrow{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{x}},i}$ by Chen and Oliver [2013b]. Here, instead, by simplifying the expression of the regression, it is shown how to compute ${\mathbf{{Y}}}_{i}$ without first computing ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$.

3.6.1 The transform matrix

Inserting the regression ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ 19b into the definition 22,

[TABLE]

where ${\mathbf{{T}}}_{i}^{+}={\mathbf{{X}}}_{i}^{+}{\mathbf{{X}}}$ has been defined, apparently requiring the pseudo-inversion of ${\mathbf{{X}}}_{i}$ for each $i$. But, as shown in section A.2,

[TABLE]

which only requires the one-time pseudo-inversion of the prior anomalies, ${\mathbf{{X}}}$. Then, since the pseudo-inversion of ${\mathbf{{T}}}_{i}\in\mathbb{R}^{N\times N}$ for ${\mathbf{{Y}}}_{i}$ 29 is a relatively small calculation, this saves computational time.

The symbol ${\mathbf{{T}}}$ has been chosen in reference to deterministic, square-root EnKFs. Indeed, multiplying equation 30 on the left by ${\mathbf{{X}}}$ and recalling equations 17 and 1 produces ${\mathbf{{X}}}_{i}={\mathbf{{X}}}{\mathbf{{T}}}_{i}$. Therefore, the “transform matrix”, ${\mathbf{{T}}}_{i}$, describes the conditioning of the anomalies (and covariance).

Conversely, equation 29 can be seen as the “de-conditioning” of the posterior observation anomalies. This interpretation of ${\mathbf{{Y}}}_{i}$ should be contrasted to its definition 22, which presents it as the prior state anomalies “propagated” by the linearization of iteration $i$. The two approaches are known to be “mainly equivalent” in the deterministic case [Sakov et al., 2012]. To our knowledge, however, it has not been exploited for EnRML before now, possibly because the proofs (section A.2) are a little more complicated in this stochastic case.

3.6.2 From the ensemble coefficients

The ensemble matrix of iteration $i$ can be written:

[TABLE]

where the columns of ${\mathbf{{W}}}_{i}\in\mathbb{R}^{N\times N}$ are the ensemble coefficient vectors 26. Multiplying equation 31 on the right by ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$ to get the anomalies produces:

[TABLE]

This seems to indicate that ${\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$ is the transform matrix, ${\mathbf{{T}}}_{i}$, discussed in the previous subsection. However, they are not fully equal: inserting ${\mathbf{{X}}}_{i}$ from 32 into 30 yields:

[TABLE]

i.e. they are distinguished by ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}={\mathbf{{X}}}^{+}{\mathbf{{X}}}$: the projection onto the row space of ${\mathbf{{X}}}$.

Section A.3 shows that, in most conditions, this pesky projection matrix vanishes when ${\mathbf{{T}}}_{i}$ is used in equation 29:

[TABLE]

In other words, the projection ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}$ can be omitted unless $\mathcal{M}{}$ is nonlinear and the ensemble is larger than the unknown state’s dimensionality.

A well known result of Reynolds et al. [2006] is that the first step of the EnRML algorithm (with ${\mathbf{{W}}}_{0}={\mathbf{{I}}}_{N}$) is equivalent to the EnKF. However, this is only strictly true if there is no appearance of ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}$ in EnRML. The following section explains why EnRML should indeed always be defined without this projection.

3.6.3 Linearization chaining

Consider applying the change of variables 26 to ${\bm{w}}$ at the very beginning of the derivation of EnRML. Since ${\mathbf{{X}}}\mathds{1}=0$, there is a redundant degree of freedom in ${\bm{w}}$, meaning that there is a choice to be made in deriving its density from the original one, given by $J_{{\bm{x}},n}({\bm{x}})$ in equation 6. The simplest choice [Bocquet et al., 2015] results in the log-posterior:

[TABLE]

Application of the Gauss-Newton scheme with the gradients and Hessian of $J_{{\bm{w}},n}$, followed by a reversion to ${\bm{x}}$, produces the same EnRML algorithm as above.

The derivation summarized in the previous paragraph is arguably simpler than that of the last few pages. Notably,

(a) it does not require the Woodbury identity to derive the subspace formulae;

(b) there is never an explicit ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ to deal with;

(c) the statistical linearization of least-squares regression from ${\mathbf{{W}}}_{i}$ to $\mathcal{M}({\mathbf{{E}}}_{i})$ directly yields equation 34, except that there are no preconditions.

While the case of a large ensemble (${N{-}1}>M$) is not typical in geoscience, the fact that this derivation does not produce a projection matrix (which requires a pseudo-inversion) under any conditions begs the questions: Why are they different? Which version is better?

The answers lie in understanding the linearization of the map ${\bm{w}}\mapsto\mathcal{M}({\bm{{\accentset{\rule{2.79996pt}{0.56001pt}}{x}}}}+{\mathbf{{X}}}{\bm{w}})$, and noting that, similarly to analytical (infinitesimal) derivatives, the chain rule applies for least-squares regression. In effect, the product ${\mathbf{{Y}}}_{i}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$, which implicitly contains the projection matrix ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}$, can be seen as an application of the chain rule for the composite function $\mathcal{M}({\bm{x}}({\bm{w}}))$. By contrast, equation 34 – but without the precondition – is obtained by direct regression of the composite function. Typically, the two versions yield identical results (i.e. the chain rule). However, since the intermediate space, $\operatorname{col}({\mathbf{{X}}})$, is of lower dimensions than the initial domain ($M<{N{-}1}$), composite linearization results in a loss of information, manifested by the projection matrix. Therefore, the definition ${\mathbf{{Y}}}_{i}=\mathcal{M}({\mathbf{{E}}}_{i})\,({\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp})^{+}$ is henceforth preferred to ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$.

Numerical experiments, as in section 4 but not shown, indicate no statistically significant advantage for either version. This corroborates similar findings by Sakov et al. [2012] for the deterministic flavour. Nevertheless, there is a practical advantage: avoiding the computation of ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}$.

3.6.4 Inverting the transform

In square-root ensemble filters, the transform matrix should have $\mathds{1}$ as an eigenvector [Sakov and Oke, 2008; Livings et al., 2008]. By construction, this also holds true for ${\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$, with eigenvalue 0. Now, consider adding ${\mathbf{{0}}}={\mathbf{{X}}}{{\mathbf{{\Pi}}}_{\mathds{1}}}$ to equation 32, yielding another valid transformation:

[TABLE]

The matrix ${\mathbf{{\Omega}}}_{i}$, in contrast to ${\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$ and ${\mathbf{{T}}}_{i}$, has eigenvalue 1 for $\mathds{1}$ and is thus invertible. This is used to prove equation 34 in section A.3, where ${\mathbf{{Y}}}_{i}$ is expressed in terms of ${\mathbf{{\Omega}}}_{i}^{-1}$.

Numerically, the use of ${\mathbf{{\Omega}}}_{i}$ in the computation 34 of ${\mathbf{{Y}}}_{i}$ was found to yield stable convergence of the new EnRML algorithm in the trivial example of $\mathcal{M}({\bm{x}})=\alpha{\bm{x}}$. By contrast, the use of $({\mathbf{{W}}}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp})^{+}$ exhibited geometrically growing (in $i$) errors when $\alpha>1$. Other formulae for the inversion are derived in section A.4; the one found to be the most stable is $({\mathbf{{W}}}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp})^{+}={\mathbf{{W}}}^{-1}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$; it is therefore preferred in LABEL:algo:GN_EnRML.

Irrespective of the inverse transform formula used, it is important to retain all non-zero singular values. This absence of a truncation threshold is a tuning simplification compared with the old EnRML algorithm, where ${\mathbf{{X}}}$ and/or ${\mathbf{{X}}}_{i}$ was scaled, decomposed, and truncated. If, by extreme chance or poor numerical subroutines, the matrix ${\mathbf{{W}}}_{i}$ is not invertible (this never occurred in any of the experiments except by our explicit intervention; cf. the conjecture in appendix A), its pseudo-inversion should be used; however, this must also be accounted for in the prior increment by multiplying the formula on line 8 on the left by the projection onto ${\mathbf{{W}}}_{i}$.

3.7 Algorithm

To summarize, LABEL:algo:GN_EnRML provides pseudo-code for the new EnRML formulation. The increments $\accentset{\rule{3.5pt}{0.62997pt}}{{\mathbf{{\Delta}}}}^{\textnormal{{lklhd}}}$ 25 and $\accentset{\rule{3.5pt}{0.62997pt}}{{\mathbf{{\Delta}}}}^{\textnormal{{prior}}}$ 28 can be recognized by multiplying line 10 on the left by ${\mathbf{{X}}}$. For aesthetics, the sign of the gradients has been reversed. Note that there is no need for an explicit iteration index. Nor is there an ensemble index, $n$, since all $N$ columns are stacked into the matrix ${\mathbf{{W}}}$. However, in case $M$ is large, ${\mathbf{{Y}}}$ may be computed column-by-column to avoid storing ${\mathbf{{E}}}$.

Line 6 is typically computed by solving ${\mathbf{{Y}}}^{\prime}{\mathbf{{W}}}=\mathcal{M}({\mathbf{{E}}})$ for ${\mathbf{{Y}}}^{\prime}$ and then subtracting its column mean. Alternative formulae are discussed in section 3.6.4. Line 9 may be computed using a reduced (or even truncated) SVD of ${\mathbf{{C}}}_{\bm{\delta}}^{-1/2}{\mathbf{{Y}}}$, which is relatively fast for $N$ both larger and smaller than $P$. Alternatively, the Kalman gain forms could be used.

The Levenberg-Marquardt variant is obtained by adding the trust-region parameter $\lambda>0$ to $({N{-}1})$ in the Hessian, line 9, which impacts both the step length and direction.

Localization may be implemented by local analysis [Hunt et al., 2007; Sakov and Bertino, 2011]; also see Bocquet [2016]; Chen and Oliver [2017]. Here, tapering is applied by replacing the local-domain ${\mathbf{{C}}}_{\bm{\delta}}^{-1/2}$ (implicit on lines 7 and 9) by $\rho\circ{\mathbf{{C}}}_{\bm{\delta}}^{-1/2}$, with $\circ$ being the Schur product, and $\rho$ a square matrix containing the (square-root) tapering coefficients, $\rho_{m,l}\in[0,1]$. If the number of local domains used is large, so that the number of ${\mathbf{{W}}}$ matrices used becomes large, then it may be more efficient to revert to the original state variables, and explicitly compute the sensitivities ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$ using the local parts of $\mathcal{M}({\mathbf{{E}}}_{i})$ and ${\mathbf{{X}}}_{i}$.

Inflation and model error parameterizations are not included in the algorithm, but may be applied outside of it. We refer to Sakov et al. [2018]; Evensen [2019] for model error treatment with iterative methods.

4 Benchmark experiments

The new EnRML algorithm produces results that are identical to the old formulation, at least up to round-off and truncation errors, and for $N-1\leq M$. Therefore, since there are already a large number of studies of EnRML with reservoir cases [e.g. Chen and Oliver, 2013a; Emerick and Reynolds, 2013b], adding to this does not seem necessary.

However, there do not appear to be any studies of EnRML with the Lorenz-96 system [Lorenz, 1996] in a data assimilation setting. The advantages of this case are numerous:

(a) the model is a surrogate of weather dynamics, and as such holds relevance in geoscience;

(b) the problem is (exhaustively) sampled from the system’s invariant measure, rather than being selected by the experimenter;

(c) the sequential nature of data assimilation inherently tests prediction skill, which helps avoid the pitfalls of point measure assessment, such as overfitting;

(d) its simplicity enhances reliability and reproducibility, and has made it a literature standard, thus facilitating comparative studies.

Comparison of the benchmark performance of EnRML will be made to the IEnKS, and to ensemble multiple data assimilation (ES-MDA)222Note that this is MDA in the sense of Emerick and Reynolds [2013a]; Stordal [2015]; Kirkpatrick et al. [1983], where the annealing itself yields iterations, and not in the sense of quasi-static assimilation [Pires et al., 1996; Bocquet and Sakov, 2014; Fillion et al., 2018], where it is used as an auxiliary technique. . Both the stochastic and the deterministic (square-root) flavours of ES-MDA are included, which in the case of only one iteration (not shown), result in exactly the same ensembles as EnRML and IEnKS, respectively. Not included in the benchmark comparisons is the version of EnRML where the prior increment is dropped (cf. section 1.1). This is because the chaotic, sequential nature of this case makes it practically impossible to achieve good results without propagating prior information. Similarly, as they lack a dynamic prior, this precludes “regularizing, iterative ensemble smoothers” [Iglesias, 2015], [Luo et al., 2015],333Their Lorenz-96 experiment only concerns the initial conditions. [Mandel et al., 2016]444Their Lorenz-96 experiment seems to have failed completely, with most of the benchmark scores (their Figure 5) indicating divergence, which makes it pointless to compare benchmarks. Also, when reproducing their experiment, we obtain much lower scores than they report for the EnKF. One possible explanation is that we include, and tune, inflation. , even if their background is well-tuned, and their stopping condition judicious. Because they require the tangent-linear model, ${\mathbf{{M}}}_{\times}$, RML and EDA/En4DVar [Tian et al., 2008; Bonavita et al., 2012; Jardak and Talagrand, 2018] are not included. For simplicity, localization will not be used, nor covariance hybridization. Other, related methods may be found in the reviews of Bannister [2017]; Carrassi et al. [2018].

4.1 Setup

The performances of the iterative ensemble smoother methods are benchmarked with “twin experiments”, using the Lorenz-96 dynamical system, which is configured with standard settings [e.g. Ott et al., 2004; Bocquet and Sakov, 2014], detailed below. The dynamics are given by the $M=40$ coupled ordinary differential equations:

[TABLE]

for $m=1,\ldots,M$, with periodic boundary conditions. These are integrated using the fourth-order Runge-Kutta scheme, with time steps of 0.05 time units, and no model noise, to yield the truth trajectory, ${\bm{x}}(t)$. Observations of the entire state vector are taken $\Delta t_{\text{obs}}\,$ time units apart with unit noise variance, meaning ${\bm{y}}(t)={\bm{x}}(t)+{\bm{\delta}}(t)$, for each $t=k\cdot\Delta t_{\text{obs}}\,$, with $k=0,1,\ldots,20\,000$, and ${\mathbf{{C}}}_{\bm{\delta}}={\mathbf{{I}}}_{M}$.

The iterative smoothers are employed in the sequential problem of filtering, aiming to estimate ${\bm{x}}(t)$ as soon as ${\bm{y}}(t)$ comes in. In so doing, they also tackle the smoothing problem for ${\bm{x}}(t{-}\Delta t_{\text{DAW}}\,)$, where the length of the data assimilation window, $\Delta t_{\text{DAW}}\,$, is fixed at a near-optimal value [inferred from Figures 3 and 4 of Bocquet and Sakov, 2013] that is also cost efficient (i.e. short). This window is shifted by $1\cdot\Delta t_{\text{obs}}\,$ each time a new observation becomes available. A post-analysis inflation factor is tuned for optimal performance for each smoother and each ensemble size, $N$. Also, random rotations are used to generate the ensembles for the square-root variants. The number of iterations is fixed, either at $3$ or $10$. No tuning of the step length is undertaken: it is $1/3$ or $1/10$ for ES-MDA, and $1$ for EnRML and the IEnKS.

The methods are assessed by their accuracy, as measured by root-mean squared error:

[TABLE]

which is recorded immediately following each analysis of the latest observation ${\bm{y}}(t)$. The “smoothing” error [assessed with ${\bm{x}}(t{-}\Delta t_{\text{DAW}}\,)$] is also recorded. After the experiment, the instantaneous ${\text{RMSE}}(t)$ are averaged for all $t>20$. The results can be reproduced using Python-code scripts hosted online at https://github.com/nansencenter/DAPPER/tree/paper_StochIEnS. This code reproduces previously published results in the literature. For example, our benchmarks obtained with the IEnKS can be cross-referenced with the ones reported by Bocquet and Sakov [2014, Figure 7a].

4.2 Results

A table of RMSE averages is compiled for a range of $N$, and then plotted as curves for each method, in Figure 1. The upper panels report the analysis RMSE scores, while the lower panels report the smoothing RMSE scores. The smoothing scores are systematically lower, but the relative results are highly similar. Moving right among the panels increases $\Delta t_{\text{obs}}\,$, and thus the nonlinearity; naturally, all of the RMSE scores also increase. As a final “sanity check”, note that the performances of all of the ensemble methods improve with increasing $N$, which needs to be at least $15$ for tolerable performance, corresponding to the rank of the unstable subspace of the dynamics plus $1$ [Bocquet and Carrassi, 2017].

For experiments with $\Delta t_{\text{obs}}\,\leq 0.4$, using $3$ iterations is largely sufficient, since its markers are rarely significantly higher than those of $10$ iterations. On the other hand, for the highly nonlinear experiment where $\Delta t_{\text{obs}}\,=0.6$, there is a significant advantage in using $10$ iterations.

The deterministic (square-root) IEnKS and ES-MDA score noticeably lower RMSE averages than the stochastic IEnKS (i.e. EnRML) and ES-MDA, which require $N$ closer to $30$ for good performance. This is qualitatively the same result as obtained for non-iterative methods [e.g. Sakov and Oke, 2008]. Also tested (not shown) was the first-order-approximate deterministic flavour of ES-MDA [Emerick, 2018]; it performed very similarly to the square-root flavour.

Among the stochastic smoothers, the one based on Gauss-Newton (EnRML) scores noticeably lower averages than the one based on annealing (ES-MDA) – when the nonlinearity is strong ($\Delta t_{\text{obs}}\,\geq 0.4$), and for small $N$. A similar trend holds for the deterministic smoothers: the IEnKS performs better than ES-MDA for $\Delta t_{\text{obs}}\,=0.6$. The likely explanation for this result is that EnRML/IEnKS can iterate indefinitely, while ES-MDA may occasionally suffer from not “reaching” the optimum.

Furthermore, the performance of EnRML/IEnKS could possibly be improved by lowering the step lengths, to avoid causing “unphysical” states, and to avoid “bouncing around” near the optimum. The tuning of the parameter that controls the step length, (e.g. the trust-region parameter and the MDA-inflation parameter) has been the subject of several studies [Chen and Oliver, 2012; Bocquet and Sakov, 2012; Ma et al., 2017; Le et al., 2016; Rafiee and Reynolds, 2017]. However, our superficial trials with this parameter (not shown) yielded little or no improvement.

5 Summary

This paper has presented a new and simpler (on paper and computationally) formulation of the iterative, stochastic ensemble smoother known as ensemble randomized maximum likelihood (EnRML). Notably, there is no explicit computation of the sensitivity matrix ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$, while the product ${\mathbf{{Y}}}_{i}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$ is computed without any pseudo-inversions of the matrix of state anomalies. This fixes issues of noise, computational cost, and covariance localization, and there is no longer any temptation to omit the prior increment from the update. Moreover, the Levenberg-Marquardt variant is now a trivial modification of the Gauss-Newton variant.

The new EnRML formulation was obtained by improvements to the background theory and derivation. Notably, Theorem 1 established the relation of the ensemble-estimated, least-squares linear regression coefficients, ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}$, to “average sensitivity”. Section 3.6 then showed that the computation of its action on the prior anomalies, ${\mathbf{{Y}}}_{i}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}_{i}{\mathbf{{X}}}$, simplifies into a de-conditioning transformation, ${\mathbf{{Y}}}_{i}=\mathcal{M}({\mathbf{{E}}}_{i})\,{\mathbf{{T}}}_{i}^{+}$. Further computational gains resulted from expressing ${\mathbf{{T}}}_{i}$ in terms of the coefficient vectors, ${\mathbf{{W}}}_{i}$, except that it also involves the “annoying” ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}^{{\mathsf{T}}}}$. Although it usually vanishes, the appearance of this projection is likely the reason why most expositions of the EnKF do not venture to declare that its implicit linearization of $\mathcal{M}$ is that of least-squares linear regression. Section 3.6.3 showed that the projection is merely the result of using the chain rule for indirect regression to the ensemble space, and argued that it is preferable to use the direct regression of the standard EnKF.

The other focus of the derivation was rank issues, with ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}$ not assumed invertible. Using the Woodbury matrix lemma, and avoiding implicit pseudo-inversions and premature insertion of SVDs, it was shown that the rank deficiency invalidates the Hessian form of the RML update, which should be restricted to the ensemble subspace. On the other hand, the subspace form and Kalman gain form of the update remain equivalent and valid. Furthermore, Theorem 2 of appendix A proves that the ensemble does not lose rank during the updates of EnRML (or EnKF).

The paper has also drawn significantly on the theory of the deterministic counterpart to EnRML: the iterative ensemble Kalman smoother (IEnKS). Comparative benchmarks using the Lorenz-96 model with these two and the ensemble multiple data assimilation (ES-MDA) smoother were shown in section 4. In the case of small ensembles and large nonlinearity, EnRML (resp. IEnKS) achieved better accuracy than stochastic (resp. deterministic) ES-MDA. Similarly to the trend for non-iterative filters, the deterministic smoothers systematically obtained better accuracy than the stochastic smoothers.

Appendix A Proofs

A.1 Preliminary

Proof of Theorem 1.

Assume $0<|{\mathbf{{C}}}_{\bm{x}}|<\infty$, and that each element of ${\mathbf{{C}}}_{\mathcal{M}({\bm{x}}),{\bm{x}}}$ and $\mathop{}\!\mathbb{E}[\mathcal{M}^{\prime}({\bm{x}})]$ is finite. Then ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}$ is a strongly consistent estimator of ${\mathbf{{C}}}_{\bm{x}}$. Likewise, ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\mathcal{M}({\bm{x}}),{\bm{x}}}\rightarrow{\mathbf{{C}}}_{\mathcal{M}({\bm{x}}),{\bm{x}}}$ almost surely, as $N\rightarrow\infty$. Thus, since ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}={\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\mathcal{M}({\bm{x}}),{\bm{x}}}\,{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{\bm{x}}^{-1}$ for sufficiently large $N$, Slutsky’s theorem yields ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{M}}}}\rightarrow{\mathbf{{C}}}_{\mathcal{M}({\bm{x}}),{\bm{x}}}\,{\mathbf{{C}}}_{\bm{x}}^{-1}$, almost surely. The equality to $\mathop{}\!\mathbb{E}[\mathcal{M}^{\prime}({\bm{x}})]$ follows directly from “Stein’s lemma” [Liu, 1994]. ∎

Theorem 2** (EnKF rank preservation).**

The posterior ensemble’s covariance, obtained using the EnKF, has the same rank as the prior’s, almost surely (a.s.).

Proof.

The updated anomalies, both for the square-root and the stochastic EnKF, can be written ${\mathbf{{X}}}^{a}={\mathbf{{X}}}{\mathbf{{T}}}^{a}$ for some ${\mathbf{{T}}}^{a}\in\mathbb{R}^{N\times N}$.

For a deterministic EnKF, ${\mathbf{{T}}}^{a}=\sqrt{N-1}{\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}}}^{-1/2}$ for the symmetric positive definite square root of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}}}$, or an orthogonal transformation thereof [Sakov and Oke, 2008]. Hence $\operatorname{rank}({\mathbf{{X}}}^{a})=\operatorname{rank}({\mathbf{{X}}})$.

For the stochastic EnKF, equations 25 and 23 may be used to show that ${\mathbf{{T}}}^{a}=({N{-}1}){\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}}}{\mathbf{{\Upsilon}}}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$, with ${\mathbf{{\Upsilon}}}={\mathbf{{I}}}_{N}+{\mathbf{{Y}}}^{{\mathsf{T}}}{\mathbf{{C}}}_{\bm{\delta}}^{-1}{\mathbf{{D}}}/({N{-}1})$. Hence, for rank preservation, it will suffice to show that ${\mathbf{{\Upsilon}}}$ is a.s. full rank.

We begin by writing ${\mathbf{{\Upsilon}}}$ more compactly:

[TABLE]

From equations 4, 14 and 38 it can be seen that column $n$ of ${\mathbf{{Z}}}$ follows the law ${\bm{z}}_{n}\sim\mathop{}\!\mathcal{N}({\bm{0}},{\mathbf{{I}}}_{P}/({N{-}1}))$. Hence, column $n$ of ${\mathbf{{\Upsilon}}}$ follows ${\bm{\upsilon}}_{n}\sim\mathop{}\!\mathcal{N}({\bm{e}}_{n},{\mathbf{{S}}}^{{\mathsf{T}}}{\mathbf{{S}}}/({N{-}1}))$, and has sample space:

[TABLE]

Now consider, for $n=0,\ldots,N$, the hypothesis:

[TABLE]

where ${\mathbf{{\Upsilon}}}_{:n}$ denotes the first $n$ columns of ${\mathbf{{\Upsilon}}}$, and ${\mathbf{{I}}}_{n:}$ denotes the last $N-n$ columns of ${\mathbf{{I}}}_{N}$. Clearly, ${\textnormal{{H}}}_{0}$ is true. Now, suppose ${\textnormal{{H}}}_{n-1}$ is true. Then the columns of $[{\mathbf{{\Upsilon}}}_{:n-1},\ {\mathbf{{I}}}_{n-1:}]$ are all linearly independent. For column $n$, this means that ${\bm{e}}_{n}\notin\operatorname{col}([{\mathbf{{\Upsilon}}}_{:n-1},\ {\mathbf{{I}}}_{n:}])$. By contrast, from equation 39, ${\bm{e}}_{n}\in\mathcal{S}_{n}$. The existence of a point in $\mathcal{S}_{n}\setminus\operatorname{col}([{\mathbf{{\Upsilon}}}_{:n-1},\ {\mathbf{{I}}}_{n:}])$ means that

[TABLE]

Since ${\bm{\upsilon}}_{n}$ is absolutely continuous with sampling space $\mathcal{S}_{n}$, equation 40 means that the probability that ${\bm{\upsilon}}_{n}\in\operatorname{col}([{\mathbf{{\Upsilon}}}_{:n-1},\ {\mathbf{{I}}}_{n:}])$ is zero. This implies ${\textnormal{{H}}}_{n}$ a.s., establishing the induction. Identifying the final hypothesis (${\textnormal{{H}}}_{N}$) with $\operatorname{rank}({\mathbf{{\Upsilon}}})=N$ concludes the proof. ∎

A corollary of Theorem 2 and Lemma 1 is that the ensemble subspace is also unchanged by the EnKF update. Note that both the prior ensemble and the model (involved through ${\mathbf{{Y}}}$) are arbitrary in Theorem 2. However, ${\mathbf{{C}}}_{\bm{\delta}}$ is assumed invertible. The result is therefore quite different from the topic discussed by Kepert [2004]; Evensen [2004], where rank deficiency arises due to a reduced-rank ${\mathbf{{C}}}_{\bm{\delta}}$.

Conjecture 1**.**

The rank of the ensemble is preserved by the EnRML update (a.s.) and ${\mathbf{{W}}}_{i}$ is invertible.

We were not able to prove Conjecture 1, but it seems a logical extension of Theorem 2, and is supported by numerical trials. The following proofs utilize Conjecture 1, without which some projections will not vanish. Yet, even if Conjecture 1 should not hold (due to bugs, truncation, or really bad luck), LABEL:algo:GN_EnRML is still valid and optimal, as discussed in sections 3.6.3 and 3.6.4.

A.2 The transform matrix

Theorem 3**.**

$({\mathbf{{X}}}^{+}{\mathbf{{X}}}_{i})^{+}={\mathbf{{X}}}_{i}^{+}{\mathbf{{X}}}$.

Proof.

Let ${\mathbf{{T}}}={\mathbf{{X}}}^{+}{\mathbf{{X}}}_{i}$ and ${\mathbf{{S}}}={\mathbf{{X}}}_{i}^{+}{\mathbf{{X}}}$. The following shows that ${\mathbf{{S}}}$ satisfies the four properties of the Moore-Penrose characterization of the pseudo-inverse of ${\mathbf{{T}}}$:

\begin{aligned} {\mathbf{{T}}}{\mathbf{{S}}}{\mathbf{{T}}}&=({\mathbf{{X}}}^{+}{\mathbf{{X}}}_{i})({\mathbf{{X}}}_{i}^{+}{\mathbf{{X}}})({\mathbf{{X}}}^{+}{\mathbf{{X}}}_{i})\\ &={\mathbf{{X}}}^{+}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}}{\mathbf{{X}}}_{i}&[{\mathbf{{\Pi}}}_{\mathbf{{A}}}={\mathbf{{A}}}{\mathbf{{A}}}^{+}]\\ &={\mathbf{{X}}}^{+}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}}{\mathbf{{X}}}_{i}&[\text{\lx@cref{creftype~refnum}{lemm:Ak_space}}]\\ &={\mathbf{{T}}}\,.&[{\mathbf{{\Pi}}}_{\mathbf{{A}}}{\mathbf{{A}}}={\mathbf{{A}}}]\end{aligned} 2. 2.

${\mathbf{{S}}}{\mathbf{{T}}}{\mathbf{{S}}}={\mathbf{{S}}}$, as may be shown similarly to point 1. 3. 3.

${\mathbf{{T}}}{\mathbf{{S}}}={\mathbf{{X}}}^{+}{\mathbf{{X}}}$, as may be shown similarly to point 1, using Conjecture 1. The symmetry of ${\mathbf{{T}}}{\mathbf{{S}}}$ follows from that of ${\mathbf{{X}}}^{+}{\mathbf{{X}}}$. 4. 4.

The symmetry of ${\mathbf{{S}}}{\mathbf{{T}}}$ is shown as for point 3.

∎

This proof was heavily inspired by appendix A of Sakov et al. [2012]. However, our developments apply for EnRML (rather than the deterministic, square-root IEnKS). This means that ${\mathbf{{T}}}_{i}$ is not symmetric, which complicates the proof in that the focus must be on ${\mathbf{{X}}}^{+}{\mathbf{{X}}}_{i}$ rather than ${\mathbf{{X}}}_{i}^{+}$ alone. Our result also shows the equivalence of ${\mathbf{{S}}}^{+}$ and ${\mathbf{{T}}}$ in general, while the additional result of the vanishing projection matrix in the case of $N-1\leq M$ is treated separately, in section A.3.

A.3 Proof of equation 34

Lemma 2**.**

${\mathbf{{\Omega}}}_{i}$* is invertible (provided ${\mathbf{{W}}}_{i}$ is).*

Proof.

We show that ${\mathbf{{\Omega}}}_{i}{\bm{u}}\neq 0$ for any ${\bm{u}}\neq 0$, where ${\mathbf{{\Omega}}}_{i}={\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}+{{\mathbf{{\Pi}}}_{\mathds{1}}}$. For ${\bm{u}}\in\operatorname{col}(\mathds{1})$: ${\mathbf{{\Omega}}}_{i}{\bm{u}}={\bm{u}}$. For ${\bm{u}}\in\operatorname{col}(\mathds{1})^{\perp}$: ${\mathbf{{\Omega}}}_{i}{\bm{u}}={\mathbf{{W}}}_{i}{\bm{u}}\neq 0$ (Conjecture 1). ∎

Recall that equation 33 was obtained by inserting ${\mathbf{{X}}}_{i}$ in the expression 30 for ${\mathbf{{T}}}_{i}$. By contrast, the following inserts ${\mathbf{{X}}}$ from equation 35 in the expression 29 for ${\mathbf{{T}}}_{i}^{+}$, yielding ${\mathbf{{T}}}_{i}^{+}={\mathbf{{X}}}_{i}^{+}{\mathbf{{X}}}={\mathbf{{X}}}_{i}{\mathbf{{X}}}_{i}{\mathbf{{\Omega}}}_{i}^{-1}={\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}{\mathbf{{\Omega}}}_{i}^{-1}={\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}{\mathbf{{\Omega}}}_{i}^{-1}$, and hence

[TABLE]

Next, it is shown that, under certain conditions, the projection matrix ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}$ vanishes:

[TABLE]

Thereafter, equation 47 of section A.4 can be used to write ${\mathbf{{\Omega}}}_{i}^{-1}$ in terms of $({\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp})^{+}$, reducing equation 42 to 34.

The case of ${N{-}1}\leq M$

In the case of ${N{-}1}\leq M$, the null space of ${\mathbf{{X}}}$ is the range of $\mathds{1}$ [with probability 1, Muirhead, 1982, Theorem 3.1.4]. By Lemma 2, the same applies for ${\mathbf{{X}}}_{i}$, and so ${\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}$ in equation 41 reduces to ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$. $\square$

The case of linearity

Let ${\mathbf{{M}}}$ be the matrix of the observation model $\mathcal{M}$, here assumed linear: $\mathcal{M}({\mathbf{{E}}}_{i})={\mathbf{{M}}}{\mathbf{{E}}}_{i}$. By equation 41, ${\mathbf{{Y}}}_{i}={\mathbf{{M}}}{\mathbf{{E}}}_{i}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}{\mathbf{{\Omega}}}_{i}^{-1}$. But ${\mathbf{{E}}}_{i}{\mathbf{{\Pi}}}_{{\mathbf{{X}}}_{i}^{{\mathsf{T}}}}={\mathbf{{X}}}_{i}={\mathbf{{E}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}$. $\square$

A.4 Inverse transforms

Recall from equation 22 that ${\mathbf{{Y}}}_{i}\mathds{1}=0$. Therefore

[TABLE]

where ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i}$ is defined in equation 24, and the identity for ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i}$ follows from that of ${\mathbf{{\accentset{\rule{3.5pt}{0.62997pt}}{C}}}}_{{\bm{w}},i}^{-1}$. Similarly, the following identities are valid also when ${\mathbf{{W}}}_{i}$ and ${\mathbf{{W}}}_{i}^{-1}$ are swapped.

[TABLE]

Equation 44 is proven inductively (in $i$) by inserting 43 in line 10 of LABEL:algo:GN_EnRML. It enables showing 45, using ${\mathbf{{\Pi}}}_{\mathds{1}}^{\perp}={\mathbf{{I}}}_{N}-{{\mathbf{{\Pi}}}_{\mathds{1}}}$. This enables showing 46, similarly to Theorem 3. Note that this implies that ${\mathbf{{Y}}}_{i}\mathds{1}=0$ also for ${\mathbf{{Y}}}_{i}=\mathcal{M}({\mathbf{{E}}}_{i})\,({\mathbf{{W}}}_{i}{\mathbf{{\Pi}}}_{\mathds{1}}^{\perp})^{+}$, and hence that the identities of this section also hold with this definition. Equations 45 and 46 can be used to show (by multiplying with ${\mathbf{{\Omega}}}_{i}$) that

[TABLE]

Acknowledgements

The authors thank Dean Oliver, Kristian Fossum, Marc Bocquet, and Pavel Sakov for their reading and comments, and Elvar Bjarkason for his questions concerning the computation of the inverse transform matrix. This work has been funded by DIGIRES, a project sponsored by industry partners and the PETROMAKS2 programme of the Research Council of Norway.