Ensemble Kalman updates are empirical Matheron updates

1 Set‑up

Let \(x\in\mathbb R^d\) be the latent state with prior \(x\sim\mathcal N(m,\Sigma)\), and \(y=Hx+\varepsilon\) the observation with \(\varepsilon\sim\mathcal N(0,R)\). The classical Gaussian conditioning gives

\[ \begin{aligned} K &= \Sigma H^\top\,(H\Sigma H^\top + R)^{-1},\\ x\mid (y=y^*) &\sim \mathcal N\Bigl(m + K(y^*-Hm),\;\Sigma - KH\Sigma\Bigr). \end{aligned} \]

Matheron’s pathwise update draws a posterior sample without forming the posterior covariance explicitly:

\[ x_{\text{post}} \stackrel{d}{=} x + \underbrace{\mathrm{Cov}(x,y)}_{C_{xy}}\;\underbrace{\mathrm{Cov}(y,y)^{-1}}_{C_{yy}^{-1}}\,(y^*-y). \]

If we replace true covariances with empirical covariances computed from a joint ensemble \((X,Y)\), we get the empirical Matheron rule:

\[ x^{(i)}_{\text{post}} = x^{(i)} + \widehat C_{xy}\,\widehat C_{yy}^{-1}\,(y^*-y^{(i)}). \]

2 Demo

The code below does three things

1. builds a tiny 1D GP prior (just a squared‑exp kernel over a grid),
2. generates one “truth” and a noisy observation \(y^*\) at a handful of points,
3. draws an empirical posterior ensemble via the Matheron/EnKF rule and checks it against the analytic Gaussian posterior.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from livingthing.matplotlib_style import set_livingthing_style, reset_default_style

set_livingthing_style()

rng = np.random.default_rng(11)

# Dimensions
d, m, N = 60, 10, 300

# Prior: squared-exponential on a 1D grid
def se_cov(n, ell=12.0, sigma=1.0, jitter=1e-8):
    x = np.arange(n)[:, None]
    S = sigma**2 * np.exp(-0.5 * (x - x.T)**2 / ell**2)
    return S + jitter*np.eye(n)

m0 = np.zeros(d)
Sigma = se_cov(d, ell=12.0, sigma=1.0)
L = np.linalg.cholesky(Sigma)

# Observation operator (m point picks) and noise
idx = np.sort(rng.choice(d, size=m, replace=False))
H = np.zeros((m, d)); H[np.arange(m), idx] = 1.0
rho = 0.15
R = (rho**2) * np.eye(m)

# One truth and its noisy observations
x_true = m0 + L @ rng.standard_normal(d)
y_star = H @ x_true + rho * rng.standard_normal(m)

# Prior ensemble
X = m0[:, None] + L @ rng.standard_normal((d, N))

# Choose EnKF flavour
stochastic = True  # True: perturbed obs; False: deterministic/square-root

if stochastic:
    Eps = rho * rng.standard_normal((m, N))
    Y = H @ X + Eps
    gamma_diag = 1e-9 * np.ones(m)   # only tiny ridge (R already inside sample cov through Y)
else:
    Y = H @ X
    gamma_diag = np.diag(R) + 1e-9

# Centered ensembles
Xc = X - X.mean(axis=1, keepdims=True)   # d x N
Yc = Y - Y.mean(axis=1, keepdims=True)   # m x N

# Woodbury factors
sqrtN1 = np.sqrt(N - 1.0)
U = Yc / sqrtN1                          # m x N
Gamma_inv = 1.0 / gamma_diag             # length-m

# residuals in obs space
Res = (y_star[:, None] - Y)              # m x N

# Build A = U^T Γ^{-1} U (N x N) and B = U^T Γ^{-1} Res (N x N) efficiently
U_over_Gamma = U * Gamma_inv[:, None]    # m x N (diag multiply)
A = U.T @ U_over_Gamma                   # N x N
B = U.T @ (Res * Gamma_inv[:, None])     # N x N

# Solve (I + A) Q = B  for Q, then ΔX = (1/sqrt(N-1)) Xc @ Q
I_N = np.eye(N)
M = np.linalg.inv(I_N + A)
Q = M @ B                                # N x N
dX = (Xc @ Q) / sqrtN1                   # d x N

X_post_enkf = X + dX

# Analytic posterior ensemble for comparison
S = H @ Sigma @ H.T + R
K = Sigma @ H.T @ np.linalg.inv(S)
mu_post = m0 + K @ (y_star - H @ m0)
Sigma_post = Sigma - K @ H @ Sigma
L_post = np.linalg.cholesky(Sigma_post + 1e-9*np.eye(d))
X_post_analytic = mu_post[:, None] + L_post @ rng.standard_normal((d, N))

# Plotting (single axes)
alpha = float(N)**(-0.75)
alpha = max(min(alpha, 0.25), 0.003)
x_axis = np.arange(d)

c_enkf     = 'tab:blue'
c_analytic = 'tab:orange'
c_obs      = 'black'

y_min = min(X_post_enkf.min(), X_post_analytic.min(), y_star.min())
y_max = max(X_post_enkf.max(), X_post_analytic.max(), y_star.max())
pad = 0.05*(y_max - y_min + 1e-12)
ylims = (y_min - pad, y_max + pad)

plt.figure(figsize=(9, 4.8), dpi=120)
for i in range(N):
    plt.plot(x_axis, X_post_enkf[:, i], color=c_enkf, alpha=alpha, linewidth=1.0)
    plt.plot(x_axis, X_post_analytic[:, i], color=c_analytic, alpha=alpha, linewidth=1.0)
plt.scatter(idx, y_star, s=32, color=c_obs, marker='o', alpha=0.9, zorder=3)

proxy_enkf = Line2D([0], [0], color=c_enkf,     lw=2, alpha=0.9, label='Empirical Matheron / EnKF (Woodbury, ΔX)')
proxy_an   = Line2D([0], [0], color=c_analytic, lw=2, alpha=0.9, label='Analytic Gaussian')
proxy_obs  = Line2D([0], [0], color=c_obs, marker='o', lw=0, label='Observations (m=10)')
plt.legend(handles=[proxy_enkf, proxy_an, proxy_obs], loc='upper right', frameon=False)

plt.title(f'Posterior samples (N={N}, {"stochastic" if stochastic else "deterministic"} EnKF)')
plt.xlabel('state index')
plt.ylabel('value')
plt.ylim(*ylims)
plt.xlim(0, d-1)
plt.tight_layout()
plt.show()

# Quick numeric check
mean_rel_err = np.linalg.norm(X_post_enkf.mean(axis=1) - mu_post) / (np.linalg.norm(mu_post) + 1e-12)
cov_rel_err  = np.linalg.norm(np.cov(X_post_enkf, bias=False) - Sigma_post, 'fro') / (np.linalg.norm(Sigma_post, 'fro') + 1e-12)
print("Relative mean error (EnKF vs analytic):", f"{mean_rel_err:.3e}")
print("Relative cov  error (EnKF vs analytic):", f"{cov_rel_err:.3e}")

Relative mean error (EnKF vs analytic): 5.756e-02
Relative cov  error (EnKF vs analytic): 8.156e-02

3 How to read it

• Everything is standard linear‑Gaussian conditioning; the only twist is replacing true covariances by ensemble covariances.
• With a noisy observation ensemble \(Y=HX+\varepsilon\), the empirical \(C_{yy}\) already contains \(R\); do not add it again.
• With a noise‑free observation ensemble \(Y=HX\), add \(R\) explicitly to \(\widehat C_{yy}\) (and you’ll typically use a square‑root transform to update the ensemble).

4 Why I care

• It gives a pathwise sampler: samples are obtained by an affine residual correction—no Cholesky of the posterior needed.
• It unifies toolboxes: localisation, inflation, and square‑root tricks from EnKF port straight over to pathwise GP sampling; conversely, autodiff and sparse‑GP ideas go the other way.
• It clarifies the zoo: stochastic EnKF, conditional GP all share the same core update; the differences are entirely in how we build \(\widehat C\) and whether we perturb \(y\).

The reason I particularly care about this is that we are moving towards the world of generative AI. These days we work in samples, not densities.

If what we actually need are posterior samples, then the Matheron/EnKF rule is the “right thing”: it’s an affine pushforward that takes a prior draw \(x\) to a posterior draw \(x' = x + C_{xy}C_{yy}^{-1}(y^\* - y)\). This avoids building or storing any \(d\times d\) covariance at all, needs only matrix–vector products with \(H\) and solves in the observation space, and gives us Monte-Carlo access to any posterior functional we care about. In other words, computing the full posterior covariance is not just expensive in large \(d\); it’s usually irrelevant. The sample-to-sample map is also a reparameterization, so it’s differentiable and streams naturally as new data arrive (update the same ensemble, don’t refactor a giant matrix). We can treat the posterior as a set of paths, not a density we must summarize, and the implementation choices make themselves. See the toy code path for a minimal demonstration of this empirical posterior sampler.

5 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 2020), 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).

@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}
}

6 Appendix: Mistakes I made while writing this up

• \(K = C_{xy}\,C_{yy}^{-1}\). Don’t sneak an extra \(H^\top\) in there.
• Don’t double‑count the noise: if we drew noisy \(Y\), we don’t add \(R\) again in \(\widehat C_{yy}\).
• Regularise \(\widehat C_{yy}\) a little (e.g. \(10^{-9}I\)) if \(m>N\) or the ensemble is small.

7 References