Revising the stochastic iterative ensemble smoother

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).

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 Sect. 1.1. As emphasized in Sect. 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 Sect. 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 Algorithm 1 of Sect. 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) and 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 ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$. As detailed in Sect. 3, this is problematic because ${\stackrel{\mathrm{‾}}{\mathbf{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 matrix (beyond the standard trust-region step regularization): the prior ensemble covariance matrix is replaced by the posterior covariance (of iteration i): ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}\leftarrow {\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x},i}$. Now the Kalman gain form of the likelihood increment is “vastly simplified” because the linearization ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ only appears in the product ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x},i}{\stackrel{\mathrm{‾}}{\mathbf{M}}}_i^{\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 Eq. (10), is already invalid because it requires the inverse of ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x},i}$. Still, in their Eq. (15), the prior increment is formulated with an inversion in ensemble space and also unburdened of the explicit computation of ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$. Intermediate explanations are lacking but could be construed to involve approximate inversions. Another issue is that the pseudo-inverse of ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}$ is now required (via 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 ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ and show how the product ${\stackrel{\mathrm{‾}}{\mathbf{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 (SVDs).

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) and 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, and, 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 that of 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), and Morzfeld et al. (2018).

2.1 The inverse problem

Consider the problem of estimating the unknown, high-dimensional state (or parameter) vector x∈ℝ^M given the observation y∈ℝ^P. It is assumed that the (generic and typically nonlinear) forward observation process may be approximated by a computational model, ℳ, so that

$$\begin{array}{}\text{(1)}& \mathit{y}=\mathcal{M}\left(\mathit{x}\right)+\mathit{\delta },\end{array}$$

where the error, δ, is random and gives rise to a likelihood, p(y|x).

In the Bayesian paradigm, prior information is quantified as a probability density function (PDF) called the prior, denoted p(x), and the truth, x, is considered to be a draw thereof. The inverse problem then consists of computing and representing the posterior which, in principle, is given by pointwise multiplication:

$$\begin{array}{}\text{(2)}& p\left(\mathit{x}\right|\mathit{y})\propto p(\mathit{y}\left|\mathit{x}\right)p\left(\mathit{x}\right),\end{array}$$

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

The prior is assumed Gaussian (normal), with mean μ_x and covariance C_x, i.e.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}p\left(\mathit{x}\right)& =\mathcal{N}\left(\mathit{x}\right|{\mathit{\mu }}_{\mathit{x}},{\mathbf{C}}_{\mathit{x}})\\ & =|\mathrm{2}\mathit{\pi }{\mathbf{C}}_{\mathit{x}}{|}^{-\frac{\mathrm{1}}{\mathrm{2}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}\Vert \mathit{x}-{\mathit{\mu }}_{\mathit{x}}{\Vert }_{{\mathbf{C}}_{\mathit{x}}}^{\mathrm{2}}}.\end{array}\end{array}$$

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

The observation error, δ, is assumed drawn from

$$\begin{array}{}\text{(4)}& {\displaystyle }p\left(\mathit{\delta }\right){\displaystyle }=\mathcal{N}\left(\mathit{\delta }\right|\mathbf{0},{\mathbf{C}}_{\mathit{\delta }}),\end{array}$$

whose covariance, ${\mathbf{C}}_{\mathit{\delta }}\in {\mathbb{R}}^{P\times P}$, will always be assumed invertible. Then, assuming that δ and x are independent and recalling Eq. (1),

$$\begin{array}{}\text{(5)}& p\left(\mathit{y}\right|\mathit{x})=\mathcal{N}(\mathit{y}\left|\mathcal{M}\right(\mathit{x}),{\mathbf{C}}_{\mathit{\delta }}).\end{array}$$

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”, $\mathit{\{}{\mathit{x}}_n{\mathit{\}}}_{n=\mathrm{1}}^N$, whose members (sample points) are assumed independently drawn from it. RML is an efficient method to approximately “condition” (i.e. implement Eq. 2 on) the prior ensemble, using optimization. Firstly, an ensemble of perturbed observations, $\mathit{\{}{\mathit{y}}_n{\mathit{\}}}_{n=\mathrm{1}}^N$, is generated as ${\mathit{y}}_n=\mathit{y}+{\mathit{\delta }}_n$, where δ_n is independently drawn according to Eq. (4).

Then, the nth “randomized log posterior”, J_x,n, is defined by Bayes' rule (Eq. 2), except with the prior mean and the observation being replaced by the nth members of the prior and observation ensembles:

$$\begin{array}{}\text{(6)}& J_{\mathit{x},n}\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathit{x}-{\mathit{x}}_n{\Vert }_{{\mathbf{C}}_{\mathit{x}}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathcal{M}\left(\mathit{x}\right)-{\mathit{y}}_n{\Vert }_{{\mathbf{C}}_{\mathit{\delta }}}^{\mathrm{2}}.\end{array}$$

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 (Eq. 7a) its gradient and (Eq. 7b) its Hessian approximated by first-order model expansions, both evaluated at the current iterate, labelled x_n,i for each member n and iteration i. To simplify the notation, define ${\mathit{x}}_{\times }={\mathit{x}}_{n,i}$. Objects evaluated at x_× are similarly denoted; for instance, ${\mathbf{M}}_{\times }={\mathcal{M}}^{\prime }\left({\mathit{x}}_{\times }\right)\in {\mathbb{R}}^{P\times M}$ denotes the Jacobian matrix of ℳ evaluated at x_×, and

$$\begin{array}{}\text{(7a)}& {\displaystyle }\mathrm{\nabla }{J}_{\times }& {\displaystyle }={\mathbf{C}}_{\mathit{x}}^{-\mathrm{1}}[{\mathit{x}}_{\times }-{\mathit{x}}_n]+{\mathbf{M}}_{\times }^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}}\left[\mathcal{M}\right({\mathit{x}}_{\times })-{\mathit{y}}_n],\text{(7b)}& {\displaystyle }{\mathbf{C}}_{\times }^{-\mathrm{1}}& {\displaystyle }={\mathbf{C}}_{\mathit{x}}^{-\mathrm{1}}+{\mathbf{M}}_{\times }^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}}{\mathbf{M}}_{\times }.\end{array}$$

Application of the Gauss–Newton scheme yields

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{x}}_{n,i+\mathrm{1}}& ={\mathit{x}}_{\times }-{\mathbf{C}}_{\times }\mathrm{\nabla }{J}_{\times }\\ & ={\mathit{x}}_{\times }+{\mathbf{\Delta }}_{\times }^{\mathrm{prior}}+{\mathbf{\Delta }}_{\times }^{\mathrm{lklhd}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(9a)}& {\displaystyle }{\mathbf{\Delta }}_{\times }^{\mathrm{prior}}& {\displaystyle }={\mathbf{C}}_{\times }{\mathbf{C}}_{\mathit{x}}^{-\mathrm{1}}[{\mathit{x}}_n-{\mathit{x}}_{\times }],\text{(9b)}& {\displaystyle }{\mathbf{\Delta }}_{\times }^{\mathrm{lklhd}}& {\displaystyle }={\mathbf{C}}_{\times }{\mathbf{M}}_{\times }^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}}[{\mathit{y}}_n-\mathcal{M}({\mathit{x}}_{\times }\left)\right],\end{array}$$

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:

$$\begin{array}{}\text{(10a)}& {\displaystyle }{\mathbf{\Delta }}_{\times }^{\mathrm{prior}}& {\displaystyle }=({\mathbf{I}}_M-{\mathbf{K}}_{\times }{\mathbf{M}}_{\times })[{\mathit{x}}_n-{\mathit{x}}_{\times }],\text{(10b)}& {\displaystyle }{\mathbf{\Delta }}_{\times }^{\mathrm{lklhd}}& {\displaystyle }={\mathbf{K}}_{\times }[{\mathit{y}}_n-\mathcal{M}({\mathit{x}}_{\times }\left)\right],\end{array}$$

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,

$$\begin{array}{}\text{(11)}& {\mathbf{K}}_{\times }={\mathbf{C}}_{\mathit{x}}{\mathbf{M}}_{\times }^{\mathsf{T}}{\mathbf{C}}_{\mathit{y}}^{-\mathrm{1}},\end{array}$$

with

$$\begin{array}{}\text{(12)}& {\mathbf{C}}_{\mathit{y}}={\mathbf{M}}_{\times }{\mathbf{C}}_{\mathit{x}}{\mathbf{M}}_{\times }^{\mathsf{T}}+{\mathbf{C}}_{\mathit{\delta }}.\end{array}$$

As the subscript suggests, C_y may be identified (in the linear case) as the prior covariance of the observation, y, of Eq. (1); it is also the covariance of the innovation, y−ℳ(μ_x). Note that if P≪M, then the inversion of ${\mathbf{C}}_{\mathit{y}}\in {\mathbb{R}}^{P\times P}$ for the Kalman gain form, Eq. (10a) and (10b), is significantly cheaper than the inversion of ${\mathbf{C}}_{\times }\in {\mathbb{R}}^{M\times M}$ for the precision matrix form, Eq. (9a) and (9b).

3 EnRML

EnRML is an approximation of RML in which the ensemble is used in its own update by estimating C_x and M_×. 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, 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 the extent to which 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:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathbf{E}=\left[\begin{array}{c}{\mathit{x}}_{\mathrm{1}},\mathrm{\dots }{\mathit{x}}_n,\mathrm{\dots }{\mathit{x}}_N\end{array}\right]\in {\mathbb{R}}^{M\times N},\text{(14)}& {\displaystyle }\mathbf{D}=\left[\begin{array}{c}{\mathit{\delta }}_{\mathrm{1}},\mathrm{\dots }{\mathit{\delta }}_n,\mathrm{\dots }{\mathit{\delta }}_N\end{array}\right]\in {\mathbb{R}}^{P\times N},\end{array}$$

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 Π is an orthogonal projector if

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathbf{\Pi }\mathbf{\Pi }=\mathbf{\Pi }={\mathbf{\Pi }}^{\mathsf{T}}.\end{array}$$

For any matrix A, let Π_A denote the projector whose image is the column space of A, implying that

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathbf{\Pi }}_{\mathbf{A}}\mathbf{A}{\displaystyle }=\mathbf{A}.\end{array}$$

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

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathbf{\Pi }}_{\mathbf{A}}={\mathbf{AA}}^{+}=\left({\mathbf{A}}^{\mathsf{T}}{)}^{+}\right({\mathbf{A}}^{\mathsf{T}}).\end{array}$$

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

Now, denote 1∈ℝ^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, $\stackrel{\mathrm{‾}}{\mathit{x}}=\mathbf{E}\mathbf{1}/N$, from each column of E. It should be appreciated that this amounts to the projection

$$\begin{array}{}\text{(18)}& {\displaystyle }\mathbf{X}{\displaystyle }=\mathbf{E}-\stackrel{\mathrm{‾}}{\mathit{x}}{\mathbf{1}}^{\mathsf{T}}=\mathbf{E}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂},\end{array}$$

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

Definition 1: the ensemble subspace. The flat (i.e. affine subspace) given by $\mathit{\{}\mathit{x}\in {\mathbb{R}}^M:[\mathit{x}-\stackrel{\mathrm{‾}}{\mathit{x}}]\in \mathrm{col}(\mathbf{X}\left)\mathit{\right\}}$.

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

3.2 The constituent estimates

The ensemble estimates of C_x and M_× are the building blocks of the EnRML algorithm. The canonical estimators are used, namely the sample covariance (Eq. 19a) and the least-squares linear-regression coefficients (Eq. 19b). They are denoted with the overhead bar:

$$\begin{array}{}\text{(19a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}& {\displaystyle }=\frac{\mathrm{1}}{N-\mathrm{1}}{\mathbf{XX}}^{\mathsf{T}},\text{(19b)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{M}}}_i& {\displaystyle }=\mathcal{M}\left({\mathbf{E}}_i\right){\mathbf{X}}_i^{+}.\end{array}$$

The anomalies at iteration i are again given by ${\mathbf{X}}_i={\mathbf{E}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$, usually computed by subtraction of ${\stackrel{\mathrm{‾}}{\mathit{x}}}_i$. The matrix ℳ(E_i) is defined by the column-wise application of ℳ to the ensemble members. Conventionally, ℳ(E_i) would also be centred in Eq. (19b), i.e. multiplied on the right by ${\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$. However, this operation (and notational burden) can be neglected because ${\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{\mathbf{X}}_i^{+}={\mathbf{X}}_i^{+}$, which follows from $\mathbf{\Pi }(\mathbf{A}\mathbf{\Pi }{)}^{+}=(\mathbf{A}\mathbf{\Pi }{)}^{+}$ (valid for any matrix A and projector Π, as shown by Maciejewski and Klein, 1985).

Note that the linearization (previously M_×, now ${\stackrel{\mathrm{‾}}{\mathbf{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), and Gu and Oliver (2007). However, this intuition has not been rigorously justified.¹ This is accomplished by the following theorem.

Theorem 1regression coefficients versus derivatives

Suppose the ensemble is drawn from a Gaussian. Then

$$\begin{array}{}\text{(20)}& {\displaystyle }\underset{N\to \mathrm{\infty }}{lim}\stackrel{\mathrm{‾}}{\mathbf{M}}{\displaystyle }=\mathbb{E}\left[{\mathcal{M}}^{\prime }\right(\mathit{x}\left)\right],\end{array}$$

with “almost sure” convergence, and expectation (𝔼) in x, which has the same distribution as the ensemble members. Regularity conditions and proof may be found in Appendix A.

A corollary of Theorem 1 is that $\stackrel{\mathrm{‾}}{\mathbf{M}}\approx \frac{\mathrm{1}}{N}{\sum }_{n=\mathrm{1}}^N{\mathcal{M}}^{\prime }\left({\mathit{x}}_n\right)$, justifying the “average sensitivity, derivative, or gradient” description. The theorem applies for the ensemble of any Gaussian, and hence also holds for ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$. On the other hand, the generality of Theorem 1 is restricted by the Gaussianity assumption. Thus, for generality and precision, ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ should simply be labelled the “least-squares (linear) fit” of ℳ, based on E_i.

Note that the computation (Eq. 19b) of ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ seemingly requires calculating a new pseudo-inverse, ${\mathbf{X}}_i^{+}$, at each iteration, i; this is addressed in Sect. 3.6.

The prior covariance estimate (previously C_x, now ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}$) is not assumed invertible, in contrast to Sect. 2. It is then not possible to employ the precision matrix forms, Eq. (9a) and (9b), because ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}^{-\mathrm{1}}$ is not defined. Using the ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}^{+}$ in its stead is flawed and damaging because it is zero in the directions orthogonal to the ensemble subspace, so 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, Eq. (19a) and (19b), are now substituted into the Kalman gain form of the update (Eq. 10a and 10b to Eq. 12). The ensemble estimate of the gain matrix, denoted ${\stackrel{\mathrm{‾}}{\mathbf{K}}}_i$, thus becomes

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{\mathbf{K}}}_i& ={\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}{\stackrel{\mathrm{‾}}{\mathbf{M}}}_i^{\mathsf{T}}({\stackrel{\mathrm{‾}}{\mathbf{M}}}_i{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}{\stackrel{\mathrm{‾}}{\mathbf{M}}}_i^{\mathsf{T}}+{\mathbf{C}}_{\mathit{\delta }}{)}^{-\mathrm{1}}\\ & ={\mathbf{XY}}_i^{\mathsf{T}}({\mathbf{Y}}_i{\mathbf{Y}}_i^{\mathsf{T}}+(N-\mathrm{1}){\mathbf{C}}_{\mathit{\delta }}{)}^{-\mathrm{1}},\end{array}\end{array}$$

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 ith iterate linearization:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathbf{Y}}_i{\displaystyle }={\stackrel{\mathrm{‾}}{\mathbf{M}}}_i\mathbf{X}.\end{array}$$

A Woodbury corollary can be used to express ${\stackrel{\mathrm{‾}}{\mathbf{K}}}_i$ as

$$\begin{array}{}\text{(23)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{K}}}_i{\displaystyle }=\mathbf{X}{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}{\mathbf{Y}}_i^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}},\end{array}$$

with

$$\begin{array}{}\text{(24)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}{\displaystyle }=({\mathbf{Y}}_i^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}}{\mathbf{Y}}_i+(N-\mathrm{1}){\mathbf{I}}_N{)}^{-\mathrm{1}}.\end{array}$$

The reason for labelling this matrix with the subscript w is revealed later. For now, note that, in the common case of N≪P, the inversion in Eq. (24) is significantly cheaper than the inversion in Eq. (21). Another computational benefit is that ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}$ is non-dimensional, improving the conditioning of the optimization problem (Lorenc, 1997).

In conclusion, the likelihood increment (Eq. 10b) is now estimated as

$$\begin{array}{}\text{(25)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{\Delta }}}_{\times }^{\mathrm{lklhd}}{\displaystyle }={\stackrel{\mathrm{‾}}{\mathbf{K}}}_i[{\mathit{y}}_n-\mathcal{M}({\mathit{x}}_{\times }\left)\right].\end{array}$$

This is efficient because ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ does not explicitly appear in ${\stackrel{\mathrm{‾}}{\mathbf{K}}}_i$ (neither in Eq. 21 nor Eq. 23) even though it is implicitly present through Y_i (Eq. 22), where it multiplies X. This absence (a) is reassuring, as the product Y_i constitutes a less noisy estimate than just ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ alone (Chen and Oliver, 2012; Emerick and Reynolds, 2013b, their Figs. 2 and 27, respectively), (b) constitutes a computational advantage, as will be shown in Sect. 3.6, and (c) enables leaving the type of linearization made for ℳ unspecified, as is usually the case in EnKF literature.

3.4 Estimating the prior increment

In contrast to the likelihood increment (Eq. 10b), the Kalman gain form of the prior increment (Eq. 10a) explicitly contains the sensitivity matrix, M_×. This issue was resolved by Bocquet and Sakov (2012) in their refinement of Sakov et al. (2012) by employing the change of variables:

$$\begin{array}{}\text{(26)}& {\displaystyle }\mathit{x}\left(\mathit{w}\right)=\stackrel{\mathrm{‾}}{\mathit{x}}+\mathbf{X}\mathit{w},\end{array}$$

where w∈ℝ^N is called the ensemble “controls” (Bannister, 2016), also known as the ensemble “weights” (Ott et al., 2004) or “coefficients” (Bocquet and Sakov, 2013).

Denote w_×, an ensemble coefficient vector, such that $\mathit{x}\left({\mathit{w}}_{\times }\right)={\mathit{x}}_{\times }$, and note that x(e_n)=x_n, where e_n is the nth column of the identity matrix. Thus, $[{\mathit{x}}_n-{\mathit{x}}_{\times }]=\mathbf{X}[{\mathit{e}}_n-{\mathit{w}}_{\times }]$, and the prior increment (Eq. 10a) with the ensemble estimates becomes

$$\begin{array}{}\text{(27)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{\Delta }}}_{\times }^{\mathrm{prior}}{\displaystyle }=(\mathbf{X}-{\stackrel{\mathrm{‾}}{\mathbf{K}}}_i{\mathbf{Y}}_i)[{\mathit{e}}_n-{\mathit{w}}_{\times }],\end{array}$$

where there is no explicit ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$, which only appears implicitly through ${\mathbf{Y}}_i={\stackrel{\mathrm{‾}}{\mathbf{M}}}_i\mathbf{X}$, as defined in Eq. (22) Alternatively, applying the subspace formula (Eq. 23) and using ${\mathbf{I}}_N={\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}({\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}{)}^{-\mathrm{1}}$ yields

$$\begin{array}{}\text{(28)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{\Delta }}}_{\times }^{\mathrm{prior}}{\displaystyle }=\mathbf{X}{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}(N-\mathrm{1})[{\mathit{e}}_n-{\mathit{w}}_{\times }].\end{array}$$

3.5 Justifying the change of variables

Lemma 1 (closure). Suppose E_i is generated by EnRML. Then, each member (column) of E_i is in the (prior) ensemble subspace. Moreover, col(X_i)⊆col(X).

Lemma 1 may be proven by noting that X is the leftmost factor in ${\stackrel{\mathrm{‾}}{\mathbf{K}}}_i$ and using induction on Eqs. (10a) and (10b). Alternatively, it can be deduced (Raanes et al., 2019) as a consequence of the implicit assumption on the prior that $\mathit{x}\sim \mathcal{N}(\stackrel{\mathrm{‾}}{\mathit{x}},{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}})$. A stronger result, namely col(X_i)=col(X), is conjectured in Appendix A, but Lemma 1 is sufficient for the present purposes: it implies that ${\mathit{w}}_{\times }\in {\mathbb{R}}^N$ exists such that $\mathit{x}\left({\mathit{w}}_{\times }\right)={\mathit{x}}_{\times }$ for any ensemble member and any iteration. Thus, the lemma justifies the change of variables (Eq. 26).

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

3.6 Simplifying the regression

Recall the definition of Eq. (22): ${\mathbf{Y}}_i={\stackrel{\mathrm{‾}}{\mathbf{M}}}_i\mathbf{X}$. Avoiding the explicit computation of ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ used in this product between the iteration-i estimate ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ and the initial (prior) X was the motivation behind the modification ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}\leftarrow {\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x},i}$ by Chen and Oliver (2013b). Here, instead, by simplifying the expression of the regression, it is shown how to compute Y_i without first computing ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$.

3.6.1 The transform matrix

Inserting the regression ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ (19b) into the definition (22),

$$\begin{array}{}\text{(29)}& {\displaystyle }{\mathbf{Y}}_i{\displaystyle }=\mathcal{M}\left({\mathbf{E}}_i\right){\mathbf{T}}_i^{+},\end{array}$$

where ${\mathbf{T}}_i^{+}={\mathbf{X}}_i^{+}\mathbf{X}$ has been defined, apparently requiring the pseudo-inversion of X_i for each i. But, as shown in Appendix A2,

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathbf{T}}_i{\displaystyle }={\mathbf{X}}^{+}{\mathbf{X}}_i,\end{array}$$

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

The symbol T has been chosen in reference to deterministic, square-root EnKFs. Indeed, multiplying Eq. (30) on the left by X and recalling Eq. (17) and Lemma 1 produces X_i=XT_i. Therefore, the “transform matrix”, T_i, describes the conditioning of the anomalies (and covariance).

Conversely, Eq. (29) can be seen as the “de-conditioning” of the posterior observation anomalies. This interpretation of 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 (Appendix A2) 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 as follows:

$$\begin{array}{}\text{(31)}& {\displaystyle }{\mathbf{E}}_i=\stackrel{\mathrm{‾}}{\mathit{x}}{\mathbf{1}}^{\mathsf{T}}+{\mathbf{XW}}_i,\end{array}$$

where the columns of ${\mathbf{W}}_i\in {\mathbb{R}}^{N\times N}$ are the ensemble coefficient vectors (Eq. 26). Multiplying Eq. (31) on the right by ${\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$ to get the anomalies produces

$$\begin{array}{}\text{(32)}& {\displaystyle }{\mathbf{X}}_i=\mathbf{X}\left({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}\right).\end{array}$$

This seems to indicate that ${\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$ is the transform matrix, T_i, discussed in the previous subsection. However, they are not fully equal: inserting X_i from Eq. (32) into Eq. (30) yields

$$\begin{array}{}\text{(33)}& {\displaystyle }{\mathbf{T}}_i{\displaystyle }={\mathbf{\Pi }}_{{\mathbf{X}}^{\mathsf{T}}}\left({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}\right),\end{array}$$

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

Appendix A3 shows that, in most conditions, this pesky projection matrix vanishes when T_i is used in Eq. (29):

$$\begin{array}{}\text{(34)}& {\displaystyle }{\mathbf{Y}}_i{\displaystyle }=\mathcal{M}\left({\mathbf{E}}_i\right)({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}\text{if}\left\{\begin{array}{l}N-\mathrm{1}\le M,\text{ or}\\ \mathcal{M}\text{is linear.}\end{array}\right.\end{array}$$

In other words, the projection ${\mathbf{\Pi }}_{{\mathbf{X}}^{\mathsf{T}}}$ can be omitted unless ℳ 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 W₀=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 w at the very beginning of the derivation of EnRML. Since X1=0, there is a redundant degree of freedom in w, meaning that there is a choice to be made in deriving its density from the original one, given by J_x,n(x) in Eq. (6). The simplest choice (Bocquet et al., 2015) results in the log posterior:

$${\displaystyle }{J}_{\mathit{w},n}\left(\mathit{w}\right){\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathit{w}-{\mathit{e}}_n{\Vert }_{\frac{\mathrm{1}}{N-\mathrm{1}}{\mathbf{I}}_N}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathcal{M}(\stackrel{\mathrm{‾}}{\mathit{x}}+\mathbf{X}\mathit{w})-{\mathit{y}}_n{\Vert }_{{\mathbf{C}}_{\mathit{\delta }}}^{\mathrm{2}}.$$

Application of the Gauss–Newton scheme with the gradients and Hessian matrix of J_w,n, followed by a reversion to 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 ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ to deal with, and (c) the statistical linearization of least-squares regression from W_i to ℳ(E_i) directly yields Eq. (34), except that there are no preconditions.

While the case of a large ensemble ($N-\mathrm{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 following questions. Why are they different? Which version is better?

The answers lie in understanding the linearization of the map $\mathit{w}\mapsto \mathcal{M}(\stackrel{\mathrm{‾}}{\mathit{x}}+\mathbf{X}\mathit{w})$ and noting that, similarly to analytical (infinitesimal) derivatives, the chain rule applies for least-squares regression. In effect, the product ${\mathbf{Y}}_i={\stackrel{\mathrm{‾}}{\mathbf{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 ℳ(x(w)). By contrast, Eq. (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, col(X), is of lower dimensions than the initial domain ($M<N-\mathrm{1}$), composite linearization results in a loss of information, manifested by the projection matrix. Therefore, the definition ${\mathbf{Y}}_i=\mathcal{M}\left({\mathbf{E}}_i\right)({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}$ is henceforth preferred to ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i\mathbf{X}$.

Numerical experiments, as in Sect. 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 1 as an eigenvector (Sakov and Oke, 2008; Livings et al., 2008). By construction, this also holds true for ${\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$, with eigenvalue 0. Now, consider adding 0=XΠ₁ to Eq. (32), yielding another valid transformation:

$$\begin{array}{}\text{(35)}& {\displaystyle }{\mathbf{X}}_i{\displaystyle }=\mathbf{X}\left(\underset{{\mathbf{\Omega }}_i}{\underbrace{{\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}+{\mathbf{\Pi }}_{\mathbf{1}}}}\right).\end{array}$$

The matrix Ω_i, in contrast to ${\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$ and T_i, has eigenvalue 1 for 1 and is thus invertible. This is used to prove Eq. (34) in Appendix A3, where Y_i is expressed in terms of ${\mathbf{\Omega }}_i^{-\mathrm{1}}$.

Numerically, the use of Ω_i in the computation (34) of Y_i was found to yield stable convergence of the new EnRML algorithm in the trivial example of ℳ(x)=αx. By contrast, the use of $(\mathbf{W}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}$ exhibited geometrically growing (in i) errors when α>1. Other formulae for the inversion are derived in Appendix A4. The one found to be the most stable is $(\mathbf{W}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}={\mathbf{W}}^{-\mathrm{1}}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$; it is therefore preferred in Algorithm 1.

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 X and/or X_i was scaled, decomposed, and truncated. If, by extreme chance or poor numerical subroutines, the matrix W_i is not invertible (this never occurred in any of the experiments except by our manual intervention; see 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 W_i.

3.7 Algorithm

To summarize, Algorithm 1 provides pseudo-code for the new EnRML formulation. The increments ${\stackrel{\mathrm{‾}}{\mathbf{\Delta }}}^{\mathrm{lklhd}}$ (Eq. 25) and ${\stackrel{\mathrm{‾}}{\mathbf{\Delta }}}^{\mathrm{prior}}$ (Eq. 28) can be recognized by multiplying line 10 on the left by 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 W. However, in case M is large, Y may be computed column by column to avoid storing E.

Line 6 is typically computed by solving ${\mathbf{Y}}^{\prime }\mathbf{W}=\mathcal{M}\left(\mathbf{E}\right)$ for Y^′ and then subtracting its column mean. Alternative formulae are discussed in Sect. 3.6.4. Line 9 may be computed using a reduced (or even truncated) SVD of ${\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}/\mathrm{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 λ>0 to (N−1) in the Hessian matrix (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; Bocquet, 2016; Chen and Oliver, 2017). Here, tapering is applied by replacing the local-domain ${\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}/\mathrm{2}}$ (implicit on lines 7 and 9) by $\mathit{\rho }\circ {\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}/\mathrm{2}}$, with ∘ being the Schur product and ρ a square matrix containing the (square-root) tapering coefficients, ${\mathit{\rho }}_{m,l}\in [\mathrm{0},\mathrm{1}]$. If the number of local domains used is large so that the number of W matrices used becomes large, then it may be more efficient to revert to the original state variables and explicitly compute the sensitivities ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$ using the local parts of ℳ(E_i) and 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) and 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-\mathrm{1}\le 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, and (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)². Both the stochastic and the deterministic (square-root) flavours of ES-MDA are included, which in the case of only one iteration (not shown), results in exactly the same ensembles as EnRML and IEnKS, respectively. Not included in the benchmark comparisons is the version of EnRML in which the prior increment is dropped (see Sect. 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”, as in Iglesias (2015), Luo et al. (2015)³, and Mandel et al. (2016)⁴, even if their background is well-tuned and their stopping condition judicious. Because they require the tangent-linear model, M_×, RML and EDA/En4DVar (Tian et al., 2008; Bonavita et al., 2012; Jardak and Talagrand, 2018) are not included. For simplicity, neither localization nor covariance hybridization will be used. Other related methods may be found in the reviews of Bannister (2016) and Carrassi et al. (2018).

4.1 Set-up

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:

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\mathrm{d}{x}_m}{\mathrm{d}t}}=\left(x_{m+\mathrm{1}}-x_{m-\mathrm{2}}\right)x_{m-\mathrm{1}}-x_m+F,\end{array}$$

for $m=\mathrm{1},\mathrm{\dots },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, x(t). Observations of the entire state vector are taken Δt_obs  time units apart with unit noise variance, meaning that $\mathit{y}\left(t\right)=\mathit{x}\left(t\right)+\mathit{\delta }\left(t\right)$, for each $t=k\cdot \mathrm{\Delta }{t}_{\text{obs}}$, with $k=\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{20}\mathrm{000}$ and C_δ=I_M.

The iterative smoothers are employed in the sequential problem of filtering, aiming to estimate x(t) as soon as y(t) comes in. In so doing, they also tackle the smoothing problem for x(t−Δt_DAW), where the length of the data assimilation window, Δt_DAW , is fixed at a near-optimal value (inferred from Figs. 3 and 4 of Bocquet and Sakov, 2013) that is also cost-efficient (i.e. short). This window is shifted by 1⋅Δt_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 the root-mean-square error (RMSE):

$$\begin{array}{}\text{(37)}& {\displaystyle }\text{RMSE}\left(t\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{M}}\Vert \mathit{x}\left(t\right)-\stackrel{\mathrm{‾}}{\mathit{x}}\left(t\right){\Vert }_{\mathrm{2}}^{\mathrm{2}}},\end{array}$$

which is recorded immediately following each analysis of the latest observation y(t). The “smoothing” error (assessed with x(t−Δt_DAW)) is also recorded. After the experiment, the instantaneous RMSE(t) values 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 (last access: 12 September 2019). 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 (Fig. 7a; 2014).

Figure 1 Benchmarks of filtering (a, b, c) and smoothing (d, e, f) accuracy, in three configurations of the Lorenz 96 system, plotted as functions of N. The y-axis changes resolution at y=1. Each of the four (coloured) iterative ensemble smoothers is plotted for 3 (hollow markers) and 10 (compact markers) iterations. It can be seen that the deterministic (i.e. square-root) methods systematically achieve lower RMSE averages. For perspective, the black lines at y=3.6 and y=0.94 show the average RMSE scores of the climatological mean and of the optimal interpolation method, respectively. The dotted lines show the scores of the stochastic (blue) and deterministic (red) EnKF.

4.2 Results

A table of RMSE averages is compiled for a range of N, and then plotted as curves for each method, in Fig. 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 Δt_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 Δt_obs ≤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 Δt_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 that 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 (Δt_obs ≥0.4) and for small N. A similar trend holds for the deterministic smoothers: the IEnKS performs better than ES-MDA for Δt_obs =0.6. The likely explanation for this result is that EnRML and IEnKS can iterate indefinitely, while ES-MDA may occasionally suffer from not “reaching” the optimum.

Furthermore, the performance of EnRML and 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 ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$, while the product ${\mathbf{Y}}_i={\stackrel{\mathrm{‾}}{\mathbf{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, ${\stackrel{\mathrm{‾}}{\mathbf{M}}}_i$, to “average sensitivity”. Section 3.6 then showed that the computation of its action on the prior anomalies, ${\mathbf{Y}}_i={\stackrel{\mathrm{‾}}{\mathbf{M}}}_i\mathbf{X}$, simplifies into a de-conditioning transformation, ${\mathbf{Y}}_i=\mathcal{M}\left({\mathbf{E}}_i\right){\mathbf{T}}_i^{+}$. Further computational gains resulted from expressing T_i in terms of the coefficient vectors, 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 ℳ 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 ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{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 Sect. 4. In the case of small ensembles and large nonlinearity, EnRML and IEnKS achieved better accuracy than the stochastic and deterministic, respectively, ES-MDA. Similarly to the trend for non-iterative filters, the deterministic smoothers systematically obtained better accuracy than the stochastic smoothers.

Appendix A: Proofs

A1 Preliminary

Proof of Theorem 1 . Assume that $\mathrm{0}<\left|{\mathbf{C}}_{\mathit{x}}\right|<\mathrm{\infty }$ and that each element of C_ℳ(x),x and 𝔼[ℳ^′(x)] is finite. Then ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}$ is a strongly consistent estimator of C_x. Likewise, ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathcal{M}\left(\mathit{x}\right),\mathit{x}}\to {\mathbf{C}}_{\mathcal{M}\left(\mathit{x}\right),\mathit{x}}$ almost surely, as N→∞. Thus, since $\stackrel{\mathrm{‾}}{\mathbf{M}}={\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathcal{M}\left(\mathit{x}\right),\mathit{x}}{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{x}}^{-\mathrm{1}}$ for sufficiently large N, Slutsky's theorem yields $\stackrel{\mathrm{‾}}{\mathbf{M}}\to {\mathbf{C}}_{\mathcal{M}\left(\mathit{x}\right),\mathit{x}}{\mathbf{C}}_{\mathit{x}}^{-\mathrm{1}}$ almost surely. The equality to 𝔼[ℳ^′(x)] follows directly from “Stein's lemma” (Liu, 1994).

Theorem 1EnKF 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 as X^a=XT^a for some ${\mathbf{T}}^a\in {\mathbb{R}}^{N\times N}$.

For a deterministic EnKF, ${\mathbf{T}}^a=\sqrt{N-\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w}}^{-\mathrm{1}/\mathrm{2}}$ for the symmetric positive definite square root of ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w}}$ or an orthogonal transformation thereof (Sakov and Oke, 2008). Hence rank(X^a)=rank(X).

For the stochastic EnKF, Eqs. (25) and (23) may be used to show that ${\mathbf{T}}^a=(N-\mathrm{1}){\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w}}\mathbf{{\rm Y}}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$, with $\mathbf{{\rm Y}}={\mathbf{I}}_N+{\mathbf{Y}}^{\mathsf{T}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}}\mathbf{D}/(N-\mathrm{1})$. Hence, for rank preservation, it will suffice to show that Υ is a.s. full rank.

We begin by writing Υ more compactly:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathbf{{\rm Y}}={\mathbf{I}}_N+{\mathbf{S}}^{\mathsf{T}}\mathbf{Z}\text{ with }\left\{\begin{array}{l}\mathbf{S}=(N-\mathrm{1}{)}^{-\mathrm{1}/\mathrm{2}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}/\mathrm{2}}\mathbf{Y},\\ \mathbf{Z}=(N-\mathrm{1}{)}^{-\mathrm{1}/\mathrm{2}}{\mathbf{C}}_{\mathit{\delta }}^{-\mathrm{1}/\mathrm{2}}\mathbf{D}.\end{array}\right.\end{array}$$

From Eqs. (4), (14), and (A1), it can be seen that column n of Z follows the law ${\mathit{z}}_n\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_P/(N-\mathrm{1}\left)\right)$. Hence, column n of Υ follows ${\mathit{\upsilon }}_n\sim \mathcal{N}({\mathit{e}}_n,{\mathbf{S}}^{\mathsf{T}}\mathbf{S}/(N-\mathrm{1}\left)\right)$ and has the sample space

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\mathcal{S}}_n=\mathit{\{}\mathit{\upsilon }\in {\mathbb{R}}^N:\mathit{\upsilon }={\mathit{e}}_n+{\mathbf{S}}^{\mathsf{T}}\mathit{z}\mathit{\}}.\end{array}$$

Now consider, for $n=\mathrm{0},\mathrm{\dots },N$, the hypothesis (H_n)

$$\begin{array}{}\text{(Hn)}& \mathrm{rank}\left(\right[{\mathbf{{\rm Y}}}_{:n},{\mathbf{I}}_{n:}\left]\right)=N,\end{array}$$

where Υ_:n denotes the first n columns of Υ, and I_n: denotes the last N−n columns of I_N. Clearly, H₀ is true. Now, suppose H_n−1 is true. Then the columns of $[{\mathbf{{\rm Y}}}_{:n-\mathrm{1}},{\mathbf{I}}_{n-\mathrm{1}:}]$ are all linearly independent. For column n, this means that ${\mathit{e}}_n\mathrm{col}\left(\right[{\mathbf{{\rm Y}}}_{:n-\mathrm{1}},{\mathbf{I}}_{n:}\left]\right)$. By contrast, from Eq. (A2), e_n∈𝒮_n. The existence of a point in ${\mathcal{S}}_n\setminus \mathrm{col}\left(\right[{\mathbf{{\rm Y}}}_{:n-\mathrm{1}},{\mathbf{I}}_{n:}\left]\right)$ means that

$$\begin{array}{}\text{(A4)}& {\displaystyle }\dim ({\mathcal{S}}_n\cap {\displaystyle }\mathrm{col}\left(\right[{\mathbf{{\rm Y}}}_{:n-\mathrm{1}},{\mathbf{I}}_{n:}\left]\right))<\dim ({\mathcal{S}}_n).\end{array}$$

Since υ_n is absolutely continuous with sampling space 𝒮_n, Eq. (A4) means that the probability that ${\mathit{\upsilon }}_n\in \mathrm{col}\left(\right[{\mathbf{{\rm Y}}}_{:n-\mathrm{1}},{\mathbf{I}}_{n:}\left]\right)$ is zero. This implies H_n a.s., establishing the induction. Identifying the final hypothesis (H_N) with rank(Υ)=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 Y) are arbitrary in Theorem 2. However, C_δ is assumed invertible. The result is therefore quite different from the topic discussed by Kepert (2004) and Evensen (2004), where rank deficiency arises due to a reduced-rank C_δ.

Conjecture 1: the rank of the ensemble is preserved by the EnRML update (a.s.), and 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), Algorithm 1 is still valid and optimal, as discussed in Sect. 3.6.3 and 3.6.4.

A2 The transform matrix

Theorem 2

$({\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 S satisfies the four properties of the Moore–Penrose characterization of the pseudo-inverse of T:

1. $\begin{array}{rl}{\displaystyle }\mathbf{TST}& {\displaystyle }=\left({\mathbf{X}}^{+}{\mathbf{X}}_i\right)\left({\mathbf{X}}_i^{+}\mathbf{X}\right)\left({\mathbf{X}}^{+}{\mathbf{X}}_i\right)\\ {\displaystyle }& {\displaystyle }={\mathbf{X}}^{+}{\mathbf{\Pi }}_{{\mathbf{X}}_i}{\mathbf{\Pi }}_{\mathbf{X}}{\mathbf{X}}_i& {\displaystyle }[{\mathbf{\Pi }}_{\mathbf{A}}={\mathbf{AA}}^{+}]\\ {\displaystyle }& {\displaystyle }={\mathbf{X}}^{+}{\mathbf{\Pi }}_{{\mathbf{X}}_i}{\mathbf{X}}_i& {\displaystyle }\left[\text{Lemma 1}\right]\\ {\displaystyle }& {\displaystyle }=\mathbf{T}& {\displaystyle }[{\mathbf{\Pi }}_{\mathbf{A}}\mathbf{A}=\mathbf{A}].\end{array}$

2. STS=S, as may be shown similarly to point 1.

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

4. The symmetry of ST 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 T_i is not symmetric, which complicates the proof in that the focus must be on X⁺X_i rather than ${\mathbf{X}}_i^{+}$ alone. Our result also shows the equivalence of S⁺ and T in general, while the additional result of the vanishing projection matrix in the case of $N-\mathrm{1}\le M$ is treated separately in Appendix A3.

A3 Proof of Eq. (34)

Lemma 2: Ω_i is invertible (provided that W_i is invertible).

Proof. We show that Ω_iu≠0 for any u≠0, where ${\mathbf{\Omega }}_i={\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}+{\mathbf{\Pi }}_{\mathbf{1}}$. For u∈col(1), Ω_iu=u. For $\mathit{u}\in \mathrm{col}(\mathbf{1}{)}^{⟂}$, ${\mathbf{\Omega }}_i\mathit{u}={\mathbf{W}}_i\mathit{u}\ne \mathrm{0}$ (Conjecture 1).

Recall that Eq. (33) was obtained by inserting X_i into the expression (Eq 30) for T_i. By contrast, the following inserts X from Eq. (35) into the expression (Eq. 29) for ${\mathbf{T}}_i^{+}$, yielding ${\mathbf{T}}_i^{+}={\mathbf{X}}_i^{+}\mathbf{X}={\mathbf{X}}_i{\mathbf{X}}_i{\mathbf{\Omega }}_i^{-\mathrm{1}}={\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}{\mathbf{\Omega }}_i^{-\mathrm{1}}={\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}{\mathbf{\Omega }}_i^{-\mathrm{1}}$, and hence

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\mathbf{Y}}_i{\displaystyle }=\left[\mathcal{M}\left({\mathbf{E}}_i\right){\mathbf{\Pi }}_{\mathbf{1}}^{⟂}\right]{\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}{\mathbf{\Omega }}_i^{-\mathrm{1}}.\end{array}$$

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

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\mathbf{Y}}_i{\displaystyle }=\left[\mathcal{M}\left({\mathbf{E}}_i\right){\mathbf{\Pi }}_{\mathbf{1}}^{⟂}\right]{\mathbf{\Omega }}_i^{-\mathrm{1}}.\end{array}$$

Thereafter, Eq. (A11) of Appendix A4 can be used to write ${\mathbf{\Omega }}_i^{-\mathrm{1}}$ in terms of $({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}$, reducing Eq. (A6) to (34).

The case of $N-\mathrm{1}\le M$

In the case of $N-\mathrm{1}\le M$, the null space of X is in the range of 1 (with probability 1; Muirhead, 1982, Theorem 3.1.4). By Lemma 2, the same applies for X_i; thus ${\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}$ in Eq. (A5) reduces to ${\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$.

The case of linearity

Let M be the matrix of the observation model ℳ, here assumed linear: ℳ(E_i)=ME_i. By Eq. (A5), ${\mathbf{Y}}_i={\mathbf{ME}}_i{\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}{\mathbf{\Omega }}_i^{-\mathrm{1}}$. But ${\mathbf{E}}_i{\mathbf{\Pi }}_{{\mathbf{X}}_i^{\mathsf{T}}}={\mathbf{X}}_i={\mathbf{E}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}$.

A4 Inverse transforms

Recall from Eq. (22) that Y_i1=0. Therefore

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}^{\pm \mathrm{1}}\mathbf{1}{\displaystyle }=(N-\mathrm{1}{)}^{\mp \mathrm{1}}\mathbf{1},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}$ is defined in Eq. (24), and the identity for ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}$ follows from that of ${\stackrel{\mathrm{‾}}{\mathbf{C}}}_{\mathit{w},i}^{-\mathrm{1}}$. Similarly, the following identities are valid also when W_i and ${\mathbf{W}}_i^{-\mathrm{1}}$ are swapped:

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathbf{W}}_i^{\mathsf{T}}\mathbf{1}=\mathbf{1},\text{(A9)}& {\displaystyle }{\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}={\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂},\text{(A10)}& {\displaystyle }({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}={\mathbf{W}}_i^{-\mathrm{1}}{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}.\end{array}$$

Equation (A8) is proven inductively (in i) by inserting Eq. (A7) into line 10 of Algorithm 1. It enables showing Eq. (A9), using ${\mathbf{\Pi }}_{\mathbf{1}}^{⟂}={\mathbf{I}}_N-{\mathbf{\Pi }}_{\mathbf{1}}$. This enables showing Eq. (A10), similarly to Theorem 3. Note that this implies that Y_i1=0 also for ${\mathbf{Y}}_i=\mathcal{M}\left({\mathbf{E}}_i\right)({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}$ and hence that the identities of this section also hold with this definition. Eqs. (A9) and (A10) can be used to show (by multiplying with Ω_i) that

$$\begin{array}{}\text{(A11)}& {\displaystyle }{\mathbf{\Omega }}_i^{-\mathrm{1}}{\displaystyle }=({\mathbf{W}}_i{\mathbf{\Pi }}_{\mathbf{1}}^{⟂}{)}^{+}+{\mathbf{\Pi }}_{\mathbf{1}}.\end{array}$$