Ensemble smoother algorithms are used to approximate a posterior distribution. Given a prior over parameters, a forward model (e.g. a petroleum reservoir model) and real-world observations: what is the posterior over the parameters? Ensemble-based algorithms describe the posterior with a set of samples, and the most common use-case is when the forward model is black-box and expensive to evaluate.

The notes are organized as follows: first we examine the Gauss-linear case, where we can compute an analytical posterior in arbitrarily high dimensions. The equation for the Gaussian posterior is the starting point for all ensemble methods. We derive the Ensemble Smoother (ES)—an algorithm that retains some properties of the Gauss-linear case. Finally we examine EnIF, which is similar to ES but works with precision matrices rather than covariances, and leans on the statistical theory of high-dimensional estimation.

Introduction

The Ensemble Kalman Filter was introduced by Evensen (1994)⁹ and the Ensemble Smoother was introduced by Van Leeuwen and Evensen (1996)¹⁰. Four years later, in Burgers et al. (1998)¹¹, it was shown that the original update equation was wrong—the covariance did not match in the Gauss-linear case.¹

New methods were proposed to deal with non-linearity. One such variant is the Iterative Ensemble Smoother (IES) of Chen and Oliver (2012¹²; 2013¹³). The idea behind IES is to move towards the maximum of the posterior in several small steps instead of a single large one, like a reduced step length optimization method. Another variant is the Ensemble Smoother with Multiple Data Assimilation (ESMDA), proposed in Emerick and Reynolds (2013)¹⁴. It also makes several small steps, but shortens step length implicitly by increasing the covariance matrix of the observations. Either way, the core idea is that linearizing, taking a small step, re-linearizing, taking another step, etc. is better than one linearization followed by one large step.

In high dimensions with few samples, empirical correlations are mostly noise and very little signal. This phenomenon is called spurious correlations, and all the methods above suffer from it. Many procedures have been proposed to alleviate it, such as setting all correlations that are below \(3 / \sqrt{r}\) to zero (Vossepoel et al., 2025)¹⁵. Another idea is to look towards the statistics of high-dimensional models, see for instance Hastie et al. (2009, Ch. 18)¹⁶. This is what Lunde (2025)¹⁷ does in the Ensemble Information Filter (EnIF), which uses graph properties to estimate sparse precision matrices and high-dimensional regression to estimate a sparse linearization of the forward model.

The papers above can be a bit tough to read: notation can be a bit inconsistent, some proofs are not as detailed as I’d like, sometimes there is no code to look at, etc. That is why I have decided to write these notes; maybe they can help clear up some confusion. If nothing else, you’ll get a second perspective. You can find well-documented and production-ready Python code for many of these algorithms at the GitHub repository github.com/equinor/iterative_ensemble_smoother/.

The Gaussian posterior

With a Gaussian prior, a Gaussian noise model and a linear forward model, we’re in the happy Gauss-linear case where an analytical expression for the posterior can be derived. In this section we establish notation and illustrate several expressions for the posterior.

Let \(\boldsymbol{x} \in \R^m\) be a vector of parameters. The forward model is denoted as \(h: \R^m \to \R^n\), and applying it to the \(m\) parameters returns a vector of \(n\) responses \(\boldsymbol{y} \in \R^n\).

The observation vector \(\boldsymbol{d} \in \R^n\) is equal to the response plus a random noise component \(\boldsymbol{\epsilon}_d\), which is often called the observation noise or measurement noise. The full equation is then \(\boldsymbol{d} = h(\boldsymbol{x}) + \boldsymbol{\epsilon}_d\), where \(\boldsymbol{\epsilon}_d \sim \mathcal{N}(\boldsymbol{0}, \Sigma_{\epsilon})\).

Think for instance of a petroleum reservoir model. The parameters \(\boldsymbol{x}\) are uncertain static properties of the subsurface, e.g., porosity. The forward model \(h\) is a reservoir flow simulator. It solves the fluid flow equations over the grid and reports predicted production data \(\boldsymbol{y}\), for instance oil production at a well. A single evaluation of \(h\) may take hours of compute, and we only have access to its inputs and outputs. The match is never exact, so the real-world observation \(\boldsymbol{d}\) departs from the simulator response by a noise term \(\boldsymbol{\epsilon}_d\). The inverse problem, called history matching, is: given the observed production \(\boldsymbol{d}\), how can we update our knowledge of the parameters \(\boldsymbol{x}\)? Characterizing the posterior distribution of \(\boldsymbol{x}\) given \(\boldsymbol{d}\) is the subject of the rest of these notes.

Exact posterior distribution \(p(\boldsymbol{x} \mid \boldsymbol{d})\)

In this section we assume that the prior distribution over parameters \(p(\boldsymbol{x})\) is multivariate Gaussian, and that the forward model \(h\) is linear. Then \(\boldsymbol{d} = H \boldsymbol{x} + \boldsymbol{\epsilon}_d = \boldsymbol{y} + \boldsymbol{\epsilon}_d\), where \(H \in \R^{n \times m}\) is a matrix. The posterior distribution \(p(\boldsymbol{x} \mid \boldsymbol{d})\) is also multivariate Gaussian.

In words, we will answer “Suppose I have a prior belief over parameters \(p(\boldsymbol{x})\) and I have observed an observation \(\boldsymbol{d}\)—how can I update my prior to a posterior?”

Suppose the prior over the parameters \(\boldsymbol{x}\) and the observation \(\boldsymbol{d}\), conditioned on the parameters, are given by multivariate Gaussians:

Then the posterior distribution over parameters, conditioned on the response, is also a multivariate Gaussian \(p(\boldsymbol{x} \mid \boldsymbol{d}) = \mathcal{N} (\boldsymbol{\mu}_{x \mid d}, \Sigma_{x \mid d})\). The posterior mean and covariance are:

The equations above are well known, and can for instance be found in Bishop (2006, §2.3.3)¹⁹.² In the equations above, the posterior mean \(\boldsymbol{\mu}_{x \mid d}\) is expressed in terms of precision matrices \(\Sigma_x^{-1}\) and \(\Sigma_\epsilon^{-1}\) rather than covariance matrices \(\Sigma_x\) and \(\Sigma_\epsilon\).

The Kalman gain form is an alternative way to express the update equations above, and we’ll spend the rest of this section deriving it. It’s more commonly used in the Kalman filter literature, and the update is given by:

In the equation above, \(H \boldsymbol{\mu}_{x}\) is the predicted mean, \((\boldsymbol{d} - H \boldsymbol{\mu}_{x})\) is the innovation term and \(K\) is the Kalman gain. In words: the posterior equals the prior, plus a linear combination of the innovation (residuals).

We will now derive the Kalman update in equation \(\eqref{eqn:posterior_mean_kalman}\). The first step is to show that \(K = \Sigma_{x \mid d}H^T \Sigma_{\epsilon}^{-1}\). To see this, we start by assuming that both expression are equal:

We right-multiply both sides with \(\left( \Sigma_\epsilon + H \Sigma_x H^T \right)\), then left-multiply both sides with \(\Sigma_{x \mid d}^{-1}\) and finally expand both sides to show that they both equal \(H^T + H^T \Sigma_\epsilon^{-1} H \Sigma_x H^T\).³

The second step is to realize that, using Remark 2, we can write

With these two results in hand, we start with equation \(\eqref{eqn:posterior_mean_precision}\) and show that it can indeed be transformed into equation \(\eqref{eqn:posterior_mean_kalman}\). The update equation for the mean becomes

In summary then, the Kalman-style update equations are:

We have seen that there are different ways to express the update equations: the precision form of equation \(\eqref{eqn:posterior_cov_precision}\), the Schur form of Remark 2 and the Kalman form in equation \(\eqref{eqn:posterior_mean_kalman}\). Which one should we use? It depends on what we wish to show or prove, or how we want to structure our computations. For instance, the Ensemble Smoother—which we discuss next—uses a Kalman style update with covariances.

The Ensemble Smoother

The Ensemble Smoother (ES) describes the posterior with an ensemble of realizations (i.e., a set of samples) instead of \(\boldsymbol{\mu}_{x \mid d}\) and \(\Sigma_{x \mid d}\). First we draw \(r\) realizations \(\{\boldsymbol{x}_j\}_{j=1}^{r}\) from the prior distribution \(p(\boldsymbol{x})\). Then we propagate them through the forward model \(h\) and obtain \(\boldsymbol{y}_j = h(\boldsymbol{x}_j)\) for each realization \(j=1, \ldots, r\). Finally, we perform data assimilation; a Kalman-style update that transforms the prior ensemble to the posterior ensemble by conditioning on the observation.

The upside of ES is that we don’t need knowledge of the forward model \(h\) to compute a posterior; \(h\) can be a black-box function. The downside is that the update equation is derived assuming (1) a linear \(h\), (2) Gaussian distributions and (3) infinitely many realizations. ES won’t actually sample the true, analytical posterior except in the happy case when conditions (1), (2) and (3) are all satisfied. Instead we approximate it.⁴ Since the forward model is expensive to compute, the number of realizations \(r\) is much smaller than the number of parameters \(m\) and responses \(n\). This leads to poor estimates of covariances due to spurious correlations.

We will now derive the update equation for ES, which computes the posterior realizations given the prior realizations. First, some notation: let \(X \in \R^{m \times r}\) be a matrix with \(r\) realizations (ensemble members); one per column.⁵ Let \(Y = h(X) = [ h(\boldsymbol{x}_1) \mid h(\boldsymbol{x}_2) \mid \cdots \mid h(\boldsymbol{x}_r)] \in \R^{n \times r}\) be a matrix of responses, where \(h\) is applied to each column.

We’ll spend the rest of this section showing that equation \(\eqref{eqn:ensemble_smoother_update}\) is the correct update equation in the following sense: when we use it on a linear model and the number of realizations \(r\) goes to infinity, the posterior mean and posterior covariance matches the analytical posterior.

To derive equation \(\eqref{eqn:ensemble_smoother_update}\), we start by examining the Gaussian update equation for the mean:

We use the results from Remark 1 and substitute \(\Sigma_x H^T = \Sigma_{xy}\) and \(H \Sigma_x H^T = \operatorname{Cov}\left[ H\boldsymbol{x} \right] = \operatorname{Cov}\left[ \boldsymbol{y} \right] = \Sigma_y\). This gives us Kalman gain in the form \(K = \Sigma_{xy} \left( \Sigma_\epsilon + \Sigma_y \right)^{-1}\).

The update equation for each realization \(j\), based on equation \(\eqref{eqn:ensemble_smoother_update}\), is given by

where \(\boldsymbol{\epsilon}_d^j \sim \mathcal{N}(\boldsymbol{0}, \Sigma_\epsilon)\). To see why this is the correct update and why we have to perturb the response \(\boldsymbol{y}^j\) with observation noise \(\boldsymbol{\epsilon}_d^j\), we have to prove that in the Gauss-linear case, this update equation recovers the same posterior mean and posterior covariance as the analytical equations in \(\eqref{eqn:posterior_mean_precision}\) and \(\eqref{eqn:posterior_cov_precision}\). To this end, suppose that \(h\) is linear, so \(\boldsymbol{y}^j = H \boldsymbol{x}^{j}_{\text{prior}}\).

We update each realization with the equation

where \(\boldsymbol{d}\) is no longer a random variable, but an actual real-world observation. The residual term \((\boldsymbol{d} - (H \boldsymbol{x}_{\text{prior}} + \boldsymbol{\epsilon}_d))\) is the actual value minus the predicted value.

First we show that the posterior mean equals equation \(\eqref{eqn:posterior_mean_kalman}\):

Second we show that the posterior covariance equals equation \(\eqref{eqn:posterior_cov_precision}\). More specifically, we show it to be equal to the Schur complement form derived in Remark 2. Below we use the fact that \(K\boldsymbol{d}\) is constant and that \(\operatorname{Cov} \left[ \boldsymbol{x}_{\text{prior}}, \boldsymbol{\epsilon}_d \right] = 0\) when we distribute the covariance over the terms.

The equalities shown in underbraces were all derived in Remark 1. From equation \(\eqref{eqn:kalman_gain}\), we observe that \(K \Sigma_{d} = (\Sigma_{xd} \Sigma_{d}^{-1}) \Sigma_{d} = \Sigma_{xd}\). The last term, \(K \Sigma_{d} K^T\), can therefore be re-written as \(\Sigma_{xd} K^T\) and cancels the second term. We are left with

where the last equality follows from Remark 2.

In summary, we have shown that if we use equation \(\eqref{eqn:ensemble_smoother_update}\) to update the ensemble, then the resulting posterior mean and covariance over the realizations matches the analytical posterior in the Gauss-linear case.

Computation

To implement ES we must compute equation \(\eqref{eqn:ensemble_smoother_update}\) given in Statement 1. Substituting empirical population estimates for \(\hat{\Sigma}_{xy}\) and \(\hat{\Sigma}_y\), the update equation becomes

This equation should never be implemented naïvely. For instance, the inverse should not be explicitly formed and neither should the huge cross-covariance matrix \(\bar{X}\bar{Y}^T\) (order of matrix multiplication matters).

Emerick and Reynolds (2012, Appendix A)²⁰ write about computation, and explain two tricks: (1) scaling the system of equations and (2) subspace inversion. Emerick (2016, Appendix B)²¹ also discusses computation. But their sentence “it is computationally more efficient to rescale based only on the diagonal elements” is not as clear as it could be: using only the diagonal elements means we don’t solve the original equation.

The best way I know to solve this system is to scale with the Cholesky factor of \(\Sigma_{\epsilon}\), then do a subspace inversion with a single SVD. This is different from both papers above, and is documented iterative_ensemble_smoother.

The Ensemble Information Filter

One issue with ES is that in equation \(\eqref{eqn:ensemble_smoother_update}\) we have to estimate \(\Sigma_{xy}\) and \(\Sigma_{y}\). When \(m\) and \(n\) are larger than \(r\), empirical covariance is a poor estimator that picks up spurious correlations—the update ends up being dominated by noise.

What if we avoid covariances and instead work with precision matrices and \(H\)? After all, there is a substantial literature on high-dimensional regression that we can lean on to estimate \(H\). Also, parameters with conditional independence structure are often better described by precision matrices \(\Lambda_x\) rather than covariance matrices \(\Sigma_{x}\), see Rue and Held (2005, Ch. 2)²². The following ideas are due to Lunde (2025)¹⁷: (1) estimate \(H\) with sparse high dimensional regression and (2) think in terms of precision rather than covariance. We’ll refer to his work as the Ensemble Information Filter (EnIF). The information terminology means that we parametrize in terms of precision rather than covariance. We’re discussing smoothers, so the filter does not exactly match this document—but we’ll stick with EnIF rather than inventing a new abbreviation.

Let us approximate \(\boldsymbol{y} = h(\boldsymbol{x})\) as \(\boldsymbol{y} = H \boldsymbol{x} + \boldsymbol{r}\). Here \(\boldsymbol{r}\) is the residual (the non-linear part of \(h\)) and \(H\) is a sparse matrix.⁷ If we model \(\boldsymbol{r} \sim \mathcal{N}(\boldsymbol{0}, \Sigma_{r})\), then a large \(\Sigma_{r}\) means the linear model is poor.⁸

We will now sketch why equation \(\eqref{eqn:ensemble_smoother_update_inf}\) is the correct update equation, in the sense that posterior mean and covariance matches under the Gauss-linear model. The analytical model given the assumption \(\boldsymbol{y} = H \boldsymbol{x} + \boldsymbol{r}\) becomes:

This is analogous to equation \(\eqref{eqn:model_def_gaussian}\), but with \(\Sigma_{\epsilon}\) replaced by \(\Sigma_{\epsilon} + \Sigma_{r}\). The posterior distribution \(p(\boldsymbol{x} \mid \boldsymbol{d}) = \mathcal{N} (\boldsymbol{\mu}_{x \mid d}, \Sigma_{x \mid d})\) over parameters is given by

From equation \(\eqref{eqn:push_through_kalman}\) and the result in Remark 1, we observe that we can write the Kalman gain as a function of precision matrices and \(H\) as

where \(\tilde{\Lambda}\) is defined as in equation \(\eqref{eqn:ensemble_smoother_update_precision_def}\) and \(\Lambda_x = \Sigma_{x}^{-1}\) by definition. At this point it’s quite easy to mimic the proof for the ES update equation, to show that equation \(\eqref{eqn:ensemble_smoother_update_inf}\) also is correct, in the sense that

The two differences between EnIF and ES, under the Gauss-linear assumption, are that (1) the Kalman gain relies on precision and \(H\) instead of covariance and (2) the covariance of \(\boldsymbol{r}\) is accounted for. We have no theoretical guarantees that EnIF is better than ES, in the sense that the update equation is equivalent but written in a different form. The hope is that if we are able to compute high-quality estimates of \(H\) and \(\Lambda_x\), the resulting update will also have high quality.

Computation

Here are some notes on computation of the EnIF update.

Subproblem 1 (Estimate \(\Lambda_x\))

• If the realizations in \(X_\text{prior}\) are known to be multivariate Gaussian, e.g. \(\boldsymbol{x}_\text{prior} \sim \mathcal{N}\left( \boldsymbol{x} \mid \boldsymbol{\mu}_x, \Sigma_x \right)\) with a known \(\Sigma_x\), then \(\Lambda_x = \Sigma_x^{-1}\) and no estimation is needed.
• If \(X_\text{prior}\) is given without information about its distribution, but the adjacency matrix \(G\) is given (\(G_{ij} = 1\) iff \((\Lambda_x)_{ij} \neq 0\)), then we have to solve the problem “estimate \(\Lambda_x\), given \(X_\text{prior}\) and \(G\).” This can for instance be done with the algorithm given in Hastie et al. (2009, §17.3.1)¹⁶ or the triangular transport approach given in Lunde (2025)¹⁷.
• If \(X_\text{prior}\) is given but no information about the graph structure \(G\) is given, then the graph structure can be estimated. A good starting point for information about this problem is Hastie et al. (2009, §17.3.2)¹⁶.

Subproblem 2 (Estimate \(H\))

The problem is to find \(H \in \mathbb{R}^{n \times m}\) such that \(Y \approx H X\). This can be formulated as \(n\) regressions, one for each response: \(X^T H^T \approx Y^T\) in matrix notation or \(X^T \boldsymbol{h}_{i, \cdot} \approx \boldsymbol{y}_{i, \cdot}\) in vector notation. For each response \(i\), we solve the regression problem to obtain \(\boldsymbol{h}_{i, \cdot}\), the \(i\)‘th row in \(H\).

This regression should be

• Sparse, since then the posterior precision \(\Lambda_x + H^T \tilde{\Lambda} H\) will retain most of the sparsity of the prior \(\Lambda_x\). Sparsity also means a more interpretable model.
• Fast, since both \(m\) and \(n\) can be huge.
• High dimensional, since \(r\) is small compared to \(m\) and \(n\).

Some kind of Lasso or boosted linear regression is a good choice. See Lunde (2025, Appendix E)¹⁷ for more details. The information criterion in Bühlmann (2006)²³ could perhaps be used.

Main computation

The algorithm inputs are \(X_{\text{prior}}\), an estimated \(H\), the observation noise covariance matrix \(\Sigma_{\epsilon}\) and responses \(Y := h(X_\text{prior})\).

1. Compute the residual matrix \(R = Y - HX_{\text{prior}}\) and estimate \(\Sigma_r\) from it. For instance by choosing \(\Sigma_r\) to be diagonal, i.e., taking the variance over each row (response-residual) in \(R\).
2. Compute \(\tilde{\Lambda} = \left( \Sigma_{\epsilon} + \Sigma_{r} \right)^{-1}\). If either \(\Sigma_{\epsilon}\) or \(\Sigma_{r}\) is dense, then do not explicitly compute the inverse like this, but instead use e.g. Cholesky factorization to solve systems of equations. If they are both diagonal, then exploit that fact in computation instead of actually forming the diagonal matrix.

3. Solve

   \begin{equation*} ( \Lambda_x + H^T \tilde{\Lambda} H ) Z = H^T \tilde{\Lambda} (\boldsymbol{d} \boldsymbol{1}^T - (Y + E)) \end{equation*}

   or, equivalently, the column version for each column \(j\)

   \begin{equation*} ( \Lambda_x + H^T \tilde{\Lambda} H ) \boldsymbol{z}_j = H^T \tilde{\Lambda} (\boldsymbol{d} - (\boldsymbol{y}_j + \boldsymbol{\epsilon}_d)) \end{equation*}

   If \(\tilde{\Lambda}\) is diagonal, then exploit that fact. Since \(( \Lambda_x + H^T \tilde{\Lambda} H )\) is sparse, large and symmetric positive definite, it can be solved with Cholesky or with Conjugate gradients. 4. Set \(X_{\text{posterior}} = X_{\text{prior}} + Z\).