Ensemble Kalman updates are empirical Matheron updates

1 Introduction

The Ensemble Kalman Filter (EnKF) (Evensen 2003, 2009a) is a cornerstone in data assimilation for large-scale dynamical systems due to its computational efficiency and scalability. The EnKF approximates the state estimation problem by evolving an ensemble of state vectors through the model dynamics and updating them using observational data.

Separately, the Matheron update provides a sample-based method for conditioning Gaussian random variables on observations(Doucet 2010; Wilson et al. 2020, 2021). We transform prior samples into posterior samples, which is sometimes much easier than calculating posterior densities. This approach is well-established in geostatistics and spatial statistics, but the connection to ensemble methods in data assimilation seems not to be well-known.

In this work, we establish that the ensemble update step in the EnKF is equivalent to an empirical Matheron update by putting them on a common probabilistic footing. By explicitly representing the ensemble mean and covariance using empirical approximations, we demonstrate this equivalence. Recognising this connection provides an alternative probabilistic foundation for the EnKF and suggests potential improvements in ensemble data assimilation techniques by leveraging the properties of the Matheron update. Conversely, the analytic Matheron updates in the literature could benefit from the computational optimisations of the data assimilation community.

2 Matheron Update

The Matheron update is a technique for sampling from the conditional distribution of a Gaussian random variable given observations, without explicitly computing the posterior covariance (Doucet 2010; Wilson et al. 2020, 2021) but instead updating prior samples. Spoiler: the EnKF update is the same thing for a particular kind of Gaussian model where the posterior is generated by a simulator.

Suppose we have a jointly Gaussian vector $$\begin{array}{}\text{(1)}& \left[\begin{array}{c}\mathsf{x}\\ \mathsf{y}\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}{\mathit{m}}_{\mathsf{x}}\\ {\mathit{m}}_{\mathsf{y}}\end{array}\right],\left[\begin{array}{cc}{\mathrm{C}}_{\mathsf{x}\mathsf{x}}& {\mathrm{C}}_{\mathsf{x}\mathsf{y}}\\ {\mathrm{C}}_{\mathsf{y}\mathsf{x}}& {\mathrm{C}}_{\mathsf{y}\mathsf{y}}\end{array}\right]).\end{array}$$

Then the conditional distribution $(\mathsf{x}|\mathsf{y}={\mathit{y}}^{\ast })$ is again Gaussian with mean and covariance: $$\begin{array}{}\text{(2)}& \mathbb{E}[\mathsf{x}\mid \mathsf{y}={\mathit{y}}^{\ast }]={\mathit{m}}_{\mathsf{x}}+{\mathrm{C}}_{\mathsf{x}\mathsf{y}}{\mathrm{C}}_{\mathsf{y}\mathsf{y}}^{-1}({\mathit{y}}^{\ast }-{\mathit{m}}_{\mathsf{y}})\end{array}$$

and $$\begin{array}{}\text{(3)}& {\mathrm{C}}_{\mathsf{x}\mid \mathsf{y}}={\mathrm{C}}_{\mathsf{x}\mathsf{x}}-{\mathrm{C}}_{\mathsf{x}\mathsf{y}}{\mathrm{C}}_{\mathsf{y}\mathsf{y}}^{-1}{\mathrm{C}}_{\mathsf{y}\mathsf{x}}.\end{array}$$

The Matheron update (Doucet 2010) says that drawing a sample from $(\mathsf{x}\mid \mathsf{y}={\mathit{y}}^{\ast })$ is equivalent to $$\begin{array}{}\text{(4)}& (\mathsf{x}\mid \mathsf{y}={\mathit{y}}^{\ast })\stackrel{\mathrm{d}}{=}\mathsf{x}+{\mathrm{C}}_{\mathsf{x}\mathsf{y}}{\mathrm{C}}_{\mathsf{y}\mathsf{y}}^{-1}({\mathit{y}}^{\ast }-\mathsf{y}).\end{array}$$

3 Kalman Filter

We begin by recalling that in state filtering the goal is to update our estimate of the system state when a new observation becomes available. In a general filtering problem, the objective is to form the posterior distribution $${\mu }_{{\mathsf{x}}_{t+1}\mid {\mathsf{x}}_t,{\mathit{y}}^{\ast }}$$ from the prior $${\mu }_{{\mathsf{x}}_{t+1}\mid {\mathsf{x}}_t}$$ by incorporating the new information in the observation ${\mathit{y}}^{\ast }$.

We assume a known observation operator $\mathcal{H}$ such that observations $\mathsf{y}$ are a priori related to the state $\mathsf{x}$ by $$\mathsf{y}=\mathcal{H}(\mathsf{x}).$$ Thus, there exists a joint law for the prior random vector, $$\begin{array}{}\text{(5)}& \left[\begin{array}{c}\mathsf{x}\\ \mathsf{y}\end{array}\right]=\left[\begin{array}{c}\mathsf{x}\\ \mathcal{H}(\mathsf{x})\end{array}\right],\end{array}$$ which is determined by the prior state distribution ${\mu }_{\mathsf{x}}$ and the observation operator $\mathcal{H}$.

The analysis step—in state filtering parlance—is the update at time $t$ of ${\mu }_{{\mathsf{x}}_t}$ into the posterior $${\mu }_{{\mathsf{x}}_t\mid ({\mathsf{y}}_t={\mathit{y}}_t^{\ast })},$$ i.e. incorporating the likelihood of the observation ${\mathit{y}}_t^{\ast }={\mathit{y}}_t$. Although recursive updating in $t$ is the focus of the classic Kalman filter, here we are concerned only with the observational update step. Henceforth we suppress the time index $t$ and consider an individual analysis update.

Suppose that the state and observation noise are independent, and all variates are defined over a finite-dimensional real vector space, with $$\mathit{x}\in {\mathbb{R}}^{D_{\mathsf{x}}},\mathit{y}\in {\mathbb{R}}^{D_{\mathsf{y}}}.$$ Furthermore, suppose that at the update step our prior belief about $\mathsf{x}$ is Gaussian with mean ${\mathit{m}}_{\mathsf{x}}$ and covariance ${\mathrm{C}}_{\mathsf{xx}}$, that the observation noise is centred Gaussian with covariance $\mathrm{R}$, and that the observation operator is linear with matrix $\mathrm{H}$ so that the observation is related to the state via $$\mathit{y}=\mathrm{H}\mathit{x}+\mathit{\epsilon },\mathit{\epsilon }\sim \mathcal{N}(0,\mathrm{R}).$$

Then the joint distribution of the prior state and observation is Gaussian, and equation Equation 5 implies that it is $$\left[\begin{array}{c}\mathsf{x}\\ \mathsf{y}\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}{\mathit{m}}_{\mathsf{x}}\\ {\mathit{m}}_{\mathsf{y}}\end{array}\right],\left[\begin{array}{cc}{\mathrm{C}}_{\mathsf{xx}}& {\mathrm{C}}_{\mathsf{xy}}\\ {\mathrm{C}}_{\mathsf{yx}}& {\mathrm{C}}_{\mathsf{yy}}\end{array}\right])$$ which can be rewritten as $$\mathcal{N}(\left[\begin{array}{c}{\mathit{m}}_{\mathsf{x}}\\ \mathrm{H}{\mathit{m}}_{\mathsf{x}}\end{array}\right],\left[\begin{array}{cc}{\mathrm{C}}_{\mathsf{xx}}& {\mathrm{C}}_{\mathsf{xx}}{\mathrm{H}}^{\mathrm{\top }}\\ \mathrm{H}{\mathrm{C}}_{\mathsf{xx}}& \mathrm{H}{\mathrm{C}}_{\mathsf{xx}}{\mathrm{H}}^{\mathrm{\top }}+\mathrm{R}\end{array}\right]).$$

When an observation ${\mathit{y}}^{\ast }$ is obtained, we apply the standard formulae for Gaussian conditioning (see, for example, equations for the conditional Gaussian) to obtain $$(\mathsf{x}\mid \mathsf{y}={\mathit{y}}^{\ast })\sim \mathcal{N}({\mathit{m}}_{\mathsf{x}\mid \mathsf{y}},{\mathrm{C}}_{\mathsf{x}\mid \mathsf{y}})$$ with the conditional mean $${\mathit{m}}_{\mathsf{x}\mid \mathit{y}}={\mathit{m}}_{\mathsf{x}}+\mathrm{K}({\mathit{y}}^{\ast }-{\mathit{m}}_{\mathsf{y}})={\mathit{m}}_{\mathsf{x}}+\mathrm{K}({\mathit{y}}^{\ast }-\mathrm{H}{\mathit{m}}_{\mathsf{x}}),$$ and the conditional covariance $${\mathrm{C}}_{\mathsf{x}\mid \mathsf{y}}={\mathrm{C}}_{\mathsf{xx}}-\mathrm{K}{\mathrm{C}}_{\mathsf{yx}}={\mathrm{C}}_{\mathsf{xx}}-\mathrm{K}\mathrm{H}{\mathrm{C}}_{\mathsf{xx}}.$$

The Kalman gain is defined by $$\begin{array}{}\text{(6)}& \mathrm{K}:-{\mathrm{C}}_{\mathsf{xy}}{\mathrm{H}}^{\mathrm{\top }}{\left({\mathrm{C}}_{\mathsf{yy}}\right)}^{-1},\end{array}$$ or equivalently, $$\mathrm{K}:-{\mathrm{C}}_{\mathsf{xx}}{\mathrm{H}}^{\mathrm{\top }}{(\mathrm{H}{\mathrm{C}}_{\mathsf{xx}}{\mathrm{H}}^{\mathrm{\top }}+\mathrm{R})}^{-1}.$$ For simplicity we assume a constant diagonal $\mathrm{R}$, so that $$\mathrm{R}={\rho }^2{\mathrm{I}}_{D_{\mathsf{y}}}.$$

4 Ensemble Kalman Filter

In high-dimensional or nonlinear settings, directly computing these posterior updates is often intractable. The Ensemble Kalman Filter (EnKF) addresses this by representing the belief about the state empirically via an ensemble of $N$ state vectors sampled from the prior distribution, $$\mathrm{X}=\left[\begin{array}{cccc}{\mathit{x}}^{(1)}& {\mathit{x}}^{(2)}& \cdots & {\mathit{x}}^{(N)}\end{array}\right],$$ and working with the empirical measure ${\hat{\mu }}_{\mathrm{X}}\approx {\mu }_{}$.

Since Gaussian measures are completely specified by their first two moments, we construct empirical measures that match the target in these moments. For convenience, we introduce the following notation for the matrix mean and deviations: $$\begin{array}{}\text{(7)}& \overline{\mathrm{X}}:-\frac{1}{N}\sum _{i=1}^N{\mathit{x}}^{(i)},\breve{\mathrm{X}}:-\frac{1}{\sqrt{N-1}}(\mathrm{X}-\overline{\mathrm{X}}{\mathbf{1}}^{\mathrm{\top }}),\end{array}$$ where ${\mathbf{1}}^{\mathrm{\top }}$ is a row vector of $N$ ones.

The ensemble mean (see Equation 8) and covariance (see Equation 9) are computed from the empirical measure: $$\begin{array}{}\text{(8)}& \hat{\mathbb{E}}[\mathsf{x}]:-{\mathbb{E}}_{\mathsf{x}\sim {\hat{\mu }}_{\mathrm{X}}}[\mathsf{x}]={\hat{\mathit{m}}}_{\mathsf{x}}=\overline{\mathrm{X}},\end{array}$$ and $$\begin{array}{}\text{(9)}& \begin{array}{rl}\widehat{\mathrm{Var}}(\mathsf{x}):-{\mathrm{Var}}_{\mathsf{x}\sim {\hat{\mu }}_{\mathrm{X}}}(\mathsf{x})& ={\hat{\mathrm{C}}}_{\mathsf{xx}}:-\frac{1}{N-1}\sum _{i=1}^N({\mathit{x}}^{(i)}-{\hat{\mathit{m}}}_{\mathsf{x}})({\mathit{x}}^{(i)}-{\hat{\mathit{m}}}_{\mathsf{x}}{)}^{\mathrm{\top }}\\ & =\breve{\mathrm{X}}{\breve{\mathrm{X}}}^{\mathrm{\top }}+({\xi }^2{\mathrm{I}}_{D_{\mathsf{x}}}),\end{array}\end{array}$$ where ${\xi }^2$ is a scalar constant that introduces regularization to ensure invertibility.

Abusing notation, we associate a Gaussian distribution with the empirical measure: $${\hat{\mu }}_{\mathrm{X}}\approx \mathcal{N}(\hat{\mathbb{E}}[\mathsf{x}],\widehat{\mathrm{Var}}(\mathsf{x}))=\mathcal{N}(\overline{\mathrm{X}},\breve{\mathrm{X}}{\breve{\mathrm{X}}}^{\mathrm{\top }}+{\xi }^2{\mathrm{I}}_{D_{\mathsf{x}}}).$$

We overload the observation operator to apply to ensemble matrices, writing $$\mathrm{Y}:-\mathcal{H}\mathrm{X}:-\left[\begin{array}{cccc}\mathcal{H}({\mathit{x}}^{(1)})& \mathcal{H}({\mathit{x}}^{(2)})& \dots & \mathcal{H}({\mathit{x}}^{(N)})\end{array}\right].$$ This also induces an empirical joint distribution $$\begin{array}{}\text{(10)}& {\hat{\mu }}_{\left[\begin{array}{c}\mathrm{X}\\ \mathrm{Y}\end{array}\right]}:-\mathcal{N}(\left[\begin{array}{c}\overline{\mathrm{X}}\\ \overline{\mathrm{Y}}\end{array}\right],\left[\begin{array}{cc}\breve{\mathrm{X}}{\breve{\mathrm{X}}}^{\mathrm{\top }}+{\xi }^2{\mathrm{I}}_{D_{\mathsf{x}}}& \breve{\mathrm{X}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}\\ \breve{\mathrm{Y}}{\breve{\mathrm{X}}}^{\mathrm{\top }}& \breve{\mathrm{Y}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}+{\upsilon }^2{\mathrm{I}}_{D_{\mathsf{y}}}\end{array}\right]).\end{array}$$

Here the regularization term ${\upsilon }^2{\mathrm{I}}_{D_{\mathsf{y}}}$ is introduced to ensure the empirical covariance is invertible when the ensemble of observations is rank deficient (typically when $D_{\mathsf{y}}>N$).

The Kalman gain in the ensemble setting is constructed by plugging the empirical estimates from Equation 10 into Equation 6. This yields $$\begin{array}{}\text{(11)}& \hat{\mathrm{K}}={\hat{\mathrm{C}}}_{\mathsf{xx}}{\mathcal{H}}^{\mathrm{\top }}{(\mathcal{H}{\hat{\mathrm{C}}}_{\mathsf{xx}}{\mathcal{H}}^{\mathrm{\top }}+{\upsilon }^2{\mathrm{I}}_{D_{\mathsf{y}}}+{\rho }^2{\mathrm{I}}_{D_{\mathsf{y}}})}^{-1}.\end{array}$$ Here, ${\rho }^2{\mathrm{I}}_{D_{\mathsf{y}}}$ is the observation error covariance matrix and ${\upsilon }^2{\mathrm{I}}_{D_{\mathsf{y}}}$ regularizes the empirical observation covariance. Combining these two, we define $${\gamma }^2:-{\upsilon }^2+{\rho }^2,$$ so that the effective covariance in the observation space is $${\mathrm{C}}_{\mathsf{yy}}=\breve{\mathrm{Y}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}+{\gamma }^2{\mathrm{I}}_{D_{\mathsf{y}}}.$$ We also define $${\mathrm{Y}}^{\ast }={\mathit{y}}^{\ast }{\mathbf{1}}^{\mathrm{\top }},$$ so that each column of ${\mathrm{Y}}^{\ast }$ equals the observation ${\mathit{y}}^{\ast }$. Then the analysis update for the ensemble is given by $$\begin{array}{}\text{(12)}& {\mathrm{X}}^{\prime }=\mathrm{X}+\hat{\mathrm{K}}({\mathrm{Y}}^{\ast }-\mathcal{H}\mathrm{X}).\end{array}$$ In other words, each ensemble member is updated as $${\mathit{x}}^{(i)}{}^{\prime }\leftarrow {\mathit{x}}^{(i)}+\hat{\mathrm{K}}({\mathit{y}}^{\ast }-\mathcal{H}({\mathit{x}}^{(i)})).$$ Equating moments, we see that $${\hat{\mu }}_{\mathsf{x}\sim {\mathrm{X}}^{\prime }}\approx {\mu }_{\mathsf{x}\mid \mathsf{y}={\mathit{y}}^{\ast }},$$ which justifies the EnKF update as an empirical approximation to the Gaussian posterior update.

5 Empirical Matheron Update is Equivalent to Ensemble Kalman Update

Under the substitution $$\begin{array}{rlrl}{\mathit{m}}_{\mathsf{x}}& \to \overline{\mathrm{X}},& {\mathrm{C}}_{\mathsf{xx}}& \to \breve{\mathrm{X}}{\breve{\mathrm{X}}}^{\mathrm{\top }}+{\xi }^2{\mathrm{I}}_{D_{\mathsf{x}}},\\ {\mathit{m}}_{\mathsf{y}}& \to \overline{\mathrm{Y}},& {\mathrm{C}}_{\mathsf{yy}}& \to \breve{\mathrm{Y}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}+{\gamma }^2{\mathrm{I}}_{D_{\mathsf{y}}},\\ {\mathrm{C}}_{\mathsf{xy}}& \to \breve{\mathrm{X}}{\breve{\mathrm{Y}}}^{\mathrm{\top }},& {\mathrm{C}}_{\mathsf{yx}}& \to \breve{\mathrm{Y}}{\breve{\mathrm{X}}}^{\mathrm{\top }},\end{array}$$ the Matheron update (Equation 4 from the earlier section) becomes identical to the ensemble update given in Equation 12.

6 Computational Complexity

As is well understood in the EnKF literature, one of the main advantages of the EnKF is that the analysis update can be performed using only empirical ensemble statistics, thereby avoiding the direct computation and inversion of full, high-dimensional covariance matrices. The naive Kalman update must operate on these full covariance matrices, which can be computationally prohibitive.

In practice, we calculate the gain from Equation 11 as follows: $$\begin{gathered} \begin{array}{rlr}\hat{\mathrm{K}}& ={\hat{\mathrm{C}}}_{\mathsf{xx}}{\mathcal{H}}^{\mathrm{\top }}{(\mathcal{H}{\hat{\mathrm{C}}}_{\mathsf{xx}}{\mathcal{H}}^{\mathrm{\top }}+{\gamma }^2{\mathrm{I}}_{D_{\mathsf{y}}})}^{-1 \\ }& =\breve{\mathrm{X}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}{(\breve{\mathrm{Y}}{\breve{\mathrm{Y}}}^{\mathrm{\top }}+{\gamma }^2{\mathrm{I}}_{D_{\mathsf{y}}})}^{-1 \\ }\\ & ={\gamma }^{-2}{\mathrm{I}}_{D_{\mathsf{y}}}-{\gamma }^{-4}\breve{\mathrm{Y}}{({\mathrm{I}}_N+{\gamma }^{-2}{\breve{\mathrm{Y}}}^{\mathrm{\top }}\breve{\mathrm{Y}})}^{-1}{\breve{\mathrm{Y}}}^{\mathrm{\top }}.\end{array} \end{gathered}$$ Here, ${\gamma }^2={\upsilon }^2+{\rho }^2$ aggregates the regularization term for the empirical observation covariance and the intrinsic observation noise.

Thus, the overall computational cost of the ensemble update is $$\mathcal{O}(D_{\mathsf{y}}{N}^2+N^3+D_{\mathsf{x}}{D}_{\mathsf{y}}N).$$ Since in practical applications the ensemble size $N$ is typically much smaller than the state and observation dimensions $D_{\mathsf{x}}$ and $D_{\mathsf{y}}$, these costs are significantly lower than those incurred by the naive Kalman update, which is dominated by $\mathcal{O}(D_{\mathsf{y}}^3)$ operations. By carefully ordering operations to work with the $N\times N$ system, the EnKF update is both computationally efficient and scalable to high-dimensional state spaces.

7 Implications

Understanding the EnKF as an empirical Matheron update opens up several possibilities. Both the Matheron update and Ensemble Kalman methods have been extensively studied in their respective fields. In recent years, numerous developments in Bayesian machine learning have advanced both approaches; see the many developments in Ensemble Kalman methods in the Ensemble Kalman notebook and especially the Neural networks meet Kalman filters notebook.

If we interpret the pathwise Matheron update in the context of the empirical measures used in the EnKF, it suggests the potential to import developments from one field to the other. The vast literature on the EnKF indicates that many ideas developed there might be applicable to Gaussian Process regression via the Matheron update. For instance, the clear probabilistic interpretation of the ensemble update could inspire improved regularization and localization strategies. One example is the Local Ensemble Transform Kalman Filter (Bocquet, Farchi, and Malartic n.d.), which provides state-of-the-art performance in high-dimensional filtering and system identification settings and could be adapted to the Matheron update (homework exercise).

Where to go from here? I don’t know; but cite me if this makes your life easier.

@misc{MacKinlay2025Ensemble,
  title = {The {{Ensemble Kalman Update}} Is an {{Empirical Matheron Update}}},
  author = {MacKinlay, Dan},
  year = {2025},
  month = feb,
  number = {arXiv:2502.03048},
  eprint = {2502.03048},
  publisher = {arXiv},
  doi = {10.48550/arXiv.2502.03048},
  archiveprefix = {arXiv}
}

8 References