Data assimilation: Kalman filter, ensemble Kalman filter and 4D-Var

Before this lesson: lesson 03 (the discrete adjoint of a time stepper; here it becomes the 4D-Var gradient) and lesson 05 (chaos and Lorenz-63: why forecasts drift). From linear algebra: lesson 07 (covariance matrices, the Schur complement and Gaussian conditioning). From the probability track: the multivariate Gaussian and Bayesian inference. A gentler first look at the Kalman filter is time series lesson 07. Time: about 120 minutes, plus the notebook.

---

Why this matters #

Every weather forecast starts from an analysis: the best estimate of the current state of the atmosphere. A global model state has more than $10^8$ numbers. The observations are far fewer, and many are indirect: a satellite measures radiances, not temperatures. Data assimilation (DA) makes the analysis by Bayes' rule.

• Reanalysis is the training data of ML weather models. ERA5, the global reanalysis of ECMWF, was made with 4D-Var (Hersbach et al. 2020). ML weather models such as GraphCast are trained on ERA5 (lesson 16). So DA choices shape the "truth" that these models learn.
• ML forecasts need an initial state. ML forecast models such as GraphCast start each forecast from an analysis made by a conventional DA system (Lam et al. 2023). Learning DA itself is an active research area (section 6).
• In ML language, the Kalman filter is exact inference in a linear-Gaussian state-space model, the EnKF is a Monte Carlo approximation of it, and 4D-Var is a MAP fit of the initial condition of a differentiable simulator by backpropagation through time.

1. Predict first #

2. The picture #

One variable #

A forecast says that a temperature anomaly is $x^f$, with error variance $P^f$. A thermometer says $y$, with error variance $R$. The two errors are independent and Gaussian. Bayes' rule multiplies the two Gaussian densities in $x$:

The exponent is a quadratic in $x$, so the product is again a Gaussian. Completing the square gives its mean $x^a$ and variance $P^a$:

Read these formulas in three ways.

• The gain $K \in [0, 1]$ says how much we trust the observation. $R \to 0$ gives $K \to 1$: copy the observation. $R \to \infty$ gives $K \to 0$: keep the forecast.
• $y - x^f$ is the innovation: the part of the observation that the forecast did not predict. The analysis increment $x^a - x^f$ is the fraction $K$ of it.
• Precisions add. The precision (inverse variance) of the analysis is the sum of the two precisions. So the analysis is more certain than each source alone. $P^a$ does not depend on the observed value $y$.

Interactive figure: loading…

Two variables: the covariance carries the information #

Now the state has two numbers, $x = (x_1, x_2)$, for example the temperature at two neighbouring grid points. The forecast error has a covariance matrix $P^f$. Its off-diagonal entry says how the two errors move together. An instrument observes only $x_1$: $y = Hx + v$ with $H = (1 \;\; 0)$.

If the forecast errors are positively correlated, a forecast that is too cold at point 1 is probably too cold at point 2 as well. So when the observation says "point 1 is warmer than forecast", Bayes' rule also warms point 2. The amount is set by the correlation. With zero correlation, point 2 does not move at all.

Interactive figure: loading…

Cycling: forecast, then analysis, then forecast again #

DA runs in cycles. The forecast step moves the estimate forward with the model. The uncertainty grows: the model moves and stretches the error ellipse, and the model error adds more. The analysis step brings in the new observations. The uncertainty shrinks. Then the next forecast starts from the analysis.

Interactive figure: loading…

Three families of methods #

Method	How it represents uncertainty	What it needs	Used for
Kalman filter (KF)	full covariance $P$ ($n \times n$)	a linear model	small linear systems, tracking
ensemble Kalman filter (EnKF)	$N$ model runs; $P \approx$ sample covariance	only forward runs of the model	large nonlinear models (weather, ocean)
3D-Var / 4D-Var	a fixed background covariance $B$; the analysis is an optimum	the adjoint model (4D-Var)	large models (global weather)

Section 3 builds each one from Bayes' rule.

3. Definitions and main results #

The linear-Gaussian state-space model #

For time steps $k = 1, 2, \dots$:

with $x_0 \sim \mathcal N(x^a_0, P^a_0)$, and $x_0$, all $w_k$ and all $v_k$ independent.

• $x_k \in \R^n$ is the state, and $y_k \in \R^p$ the observation (usually $p \ll n$).
• $M$ is the model (one forecast step) and $H$ the observation operator.
• $Q$ is the model-error covariance and $R$ the observation-error covariance.
• $y_{1:k}$ means all observations $y_1, \dots, y_k$.

The superscript $f$ (forecast) marks the prior at an analysis time, and $a$ (analysis) marks the posterior. In variational DA the prior is called the background, $x^b$ with covariance $B$.

Gaussian conditioning #

Why this is true (intuition). Predict $x$ from $y$ with the best linear rule. What is left over, $x$ minus the prediction, is uncorrelated with $y$. For jointly Gaussian variables, uncorrelated means independent. So knowing $y$ tells us the prediction exactly and tells us nothing about the leftover.

Theorem 1 (Gaussian conditioning). Let $(x, y)$ be jointly Gaussian with means $\mu_x, \mu_y$ and covariance blocks $\Sigma_{xx}$, $\Sigma_{xy}$, $\Sigma_{yx} = \Sigma_{xy}^\top$, $\Sigma_{yy}$, with $\Sigma_{yy}$ invertible. Then

The covariance is the Schur complement of $\Sigma_{yy}$ in the joint covariance (linear algebra 07).

Proof.

1. Let $G = \Sigma_{xy}\Sigma_{yy}^{-1}$ and $z = x - Gy$.
2. $\operatorname{Cov}(z, y) = \Sigma_{xy} - G\Sigma_{yy} = 0$. The pair $(z, y)$ is a linear map of $(x, y)$, so it is jointly Gaussian. For jointly Gaussian variables, zero covariance means independence. So $z$ and $y$ are independent.
3. So the law of $z$ given $y$ is the same as without $y$: mean $\mu_x - G\mu_y$ and covariance $\Sigma_{xx} - G\Sigma_{yx}$. (Expand $\operatorname{Cov}(x - Gy)$ and use $G\Sigma_{yy}G^\top = \Sigma_{xy}G^\top = G\Sigma_{yx}$.)
4. $x = z + Gy$, and $Gy$ is a constant once $y$ is known. Add it to the mean. $\blacksquare$

The Kalman update #

Theorem 2 (Kalman update). Prior $x \sim \mathcal N(x^f, P^f)$. Observation $y = Hx + v$ with $v \sim \mathcal N(0, R)$ independent of $x$. Then $x \mid y \sim \mathcal N(x^a, P^a)$ with

$S$ is the innovation covariance and $K$ the Kalman gain.

Proof.

1. $(x, y)$ is a linear map of the independent Gaussians $(x, v)$, so it is jointly Gaussian.
2. Compute the blocks. $\mathbb E[y] = Hx^f$. $\operatorname{Cov}(y) = HP^fH^\top + R = S$. $\operatorname{Cov}(x, y) = \mathbb E\big[(x - x^f)(H(x - x^f) + v)^\top\big] = P^fH^\top$, because $v$ has mean $0$ and is independent of $x$.
3. Theorem 1 with $\Sigma_{xy} = P^fH^\top$ and $\Sigma_{yy} = S$ gives the mean $x^f + P^fH^\top S^{-1}(y - Hx^f)$ and the covariance $P^f - P^fH^\top S^{-1}HP^f = (I - KH)P^f$. $\blacksquare$

For one variable ($n = p = 1$, $H = 1$) this is exactly the formula of section 2. Three consequences:

• An observation never increases the uncertainty. $P^f - P^a = KSK^\top$, which is positive semidefinite.
• $K$ and $P^a$ do not depend on the observed value $y$. They can be computed before the data arrive.
• The cross-covariance moves unobserved variables. Row $i$ of $K = P^fH^\top S^{-1}$ is nonzero only if variable $i$ is correlated (in the forecast) with what $H$ observes.

Information form, Woodbury and the Joseph form

Information form (for invertible $P^f$ and $R$):

$$\begin{aligned} (P^a)^{-1} &= (P^f)^{-1} + H^\top R^{-1}H, \\ x^a &= P^a\big((P^f)^{-1}x^f + H^\top R^{-1}y\big), \\ K &= P^aH^\top R^{-1}. \end{aligned}$$

Reason: $-\log p(x \mid y) = \tfrac12 (x - x^f)^\top (P^f)^{-1}(x - x^f) + \tfrac12 (y - Hx)^\top R^{-1}(y - Hx) + \text{const}$ is quadratic in $x$. Its Hessian is the inverse posterior covariance, and its minimiser is the posterior mean. The Woodbury identity $\big((P^f)^{-1} + H^\top R^{-1}H\big)^{-1} = P^f - P^fH^\top(HP^fH^\top + R)^{-1}HP^f$ turns one form into the other. The covariance form solves a $p \times p$ system; the information form an $n \times n$ system. Use the smaller one.

Joseph form. In floating point, $(I - KH)P^f$ can lose symmetry and positive definiteness over many cycles. The Joseph form $(I - KH)P^f(I - KH)^\top + KRK^\top$ is symmetric positive semidefinite for every $K$, and it equals $P^a$ when $K$ is the Kalman gain. It is also the covariance of the analysis error $x^a - x$ for any gain $K$; minimising its trace over $K$ gives the Kalman gain again.

The Kalman filter #

Start from $x^a_0$, $P^a_0$. For $k = 1, 2, \dots$:

• Forecast step: $x^f_k = M x^a_{k-1}$ and $P^f_k = M P^a_{k-1} M^\top + Q$.
• Analysis step: Theorem 2 with prior $\mathcal N(x^f_k, P^f_k)$ and observation $y_k$, which gives $x^a_k$ and $P^a_k$.

Keep the time indices straight: $x^f_k$ uses the observations $y_{1:k-1}$, and $x^a_k$ uses $y_{1:k}$.

Theorem 3 (the Kalman filter is exact). For the linear-Gaussian model, $p(x_k \mid y_{1:k}) = \mathcal N(x^a_k, P^a_k)$ for every $k$, and $p(x_k \mid y_{1:k-1}) = \mathcal N(x^f_k, P^f_k)$.

Proof by induction on k

1. At $k = 0$ there are no data, and $p(x_0) = \mathcal N(x^a_0, P^a_0)$.
2. Forecast. Assume $x_{k-1} \mid y_{1:k-1} \sim \mathcal N(x^a_{k-1}, P^a_{k-1})$. The noise $w_k$ is independent of $x_{k-1}$ and of $y_{1:k-1}$. So, given $y_{1:k-1}$, $x_k = Mx_{k-1} + w_k$ is a linear map of independent Gaussians: Gaussian with mean $Mx^a_{k-1}$ and covariance $MP^a_{k-1}M^\top + Q$.
3. Analysis. $y_k = Hx_k + v_k$, and $v_k$ is independent of $x_k$ and of $y_{1:k-1}$. So, given $y_{1:k-1}$, the pair $(x_k, y_k)$ is exactly the setting of Theorem 2 with prior $\mathcal N(x^f_k, P^f_k)$. Conditioning on $y_k$ gives $p(x_k \mid y_{1:k})$. $\blacksquare$

Step 2 of the proof (the forecast) is where linearity matters. A nonlinear model maps a Gaussian to a non-Gaussian distribution. Then every Gaussian update (the extended Kalman filter, the EnKF) is an approximation.

Why large models do not use the KF. The KF carries $P$, an $n \times n$ matrix: $10^{16}$ numbers for $n = 10^8$. Computing $MPM^\top$ needs about $2n$ runs of the (tangent-)linear model per step. This is why large systems use ensembles or variational methods.

The Rauch–Tung–Striebel smoother

The filter uses only past observations. A smoother also uses future ones: $p(x_k \mid y_{1:K})$ for $k < K$. For the linear-Gaussian model, run the filter forward and store $x^f_k, P^f_k, x^a_k, P^a_k$. Then go backward from $x^s_K = x^a_K$, $P^s_K = P^a_K$:

$$\begin{aligned} G_k &= P^a_k M^\top (P^f_{k+1})^{-1}, \\ x^s_k &= x^a_k + G_k\,(x^s_{k+1} - x^f_{k+1}), \\ P^s_k &= P^a_k + G_k\,(P^s_{k+1} - P^f_{k+1})\,G_k^\top . \end{aligned}$$

This is the RTS smoother (Särkkä and Svensson, §12.2). It gives the exact $p(x_k \mid y_{1:K})$; the notebook of this lesson does not need it, but the same forward–backward structure appears in 4D-Var below. In ML terms, filter and smoother are the forward and backward passes of message passing on a chain.

The ensemble Kalman filter #

Keep $N$ model states $x^{(1)}, \dots, x^{(N)}$, the ensemble members; typically $N$ is 20 to 100. Let $\bar x$ be their mean, and let $A \in \R^{n \times N}$ be the anomaly matrix with columns $x^{(i)} - \bar x$. The sample covariance $P^f = AA^\top/(N - 1)$ replaces the exact covariance.

• Forecast: run every member through the full nonlinear model (plus model noise, if used). No tangent-linear model is needed.
• Analysis (stochastic EnKF): draw independent $\varepsilon_i \sim \mathcal N(0, R)$ and update each member with its own perturbed observation $y + \varepsilon_i$:

In a large model, never form $P^f$. Use $P^fH^\top = A\,(HA)^\top/(N - 1)$, where $HA$ (size $p \times N$) holds the anomalies in observation space.

Why perturb the observations. Take a linear $H$ and write $a_i$ for a forecast anomaly. The analysis anomaly is $(I - KH)\,a_i + K\varepsilon_i'$, where $\varepsilon_i'$ is the perturbation minus the mean perturbation. Its expected covariance is

the correct analysis covariance (the Joseph form above). Without the perturbations the term $KRK^\top$ is missing, and the analysis spread is too small. (Square-root filters get the correct covariance without random numbers; see the folded box below.)

Three properties.

• Linear model, $N \to \infty$: the EnKF converges to the KF.
• Nonlinear model: the forecast is a Monte Carlo run of the true model, but the analysis is still a linear, Gaussian update. So the EnKF is not exact for non-Gaussian posteriors.
• Rank. $K = P^fH^\top S^{-1} = A\,C$ with $C = (HA)^\top S^{-1}/(N - 1)$. So every column of $K$, and every increment $x^{(i)}_a - x^{(i)}$, is a combination of the columns of $A$. The anomalies sum to zero ($A\mathbf 1 = 0$), so their span has dimension at most $N - 1$. With $N = 50$ and $n = 10^8$, the analysis can change the state in only 49 directions.

Sampling error and localisation. For Gaussian data, the sample correlation of two uncorrelated variables has a standard deviation of about $1/\sqrt{N - 1}$: $0.33$ for $N = 10$. So an observation in Zurich also "corrects" a grid point in Lisbon, through a correlation that is only noise. These are spurious correlations. Localisation removes them.

• Covariance localisation: replace $P^f$ by $\rho \circ P^f$ (the elementwise product), with $\rho_{jl} = \rho(d_{jl}/c)$, where $d_{jl}$ is the distance between points $j$ and $l$, and $c$ is a length scale. The taper $\rho$ is 1 at distance 0 and decays to 0. The standard choice is the Gaspari–Cohn function, a piecewise rational function of degree 5 that is exactly 0 beyond $2c$. By the Schur product theorem, $\rho \circ P^f$ is positive semidefinite when $\rho$ and $P^f$ are. It usually has full rank, so the analysis is no longer limited to $N - 1$ directions.
• Domain localisation (used by the LETKF): analyse each grid point separately, with only the observations within a radius.
• The radius is a trade-off. Too small cuts real correlations; too large lets the noise in. More members allow a larger radius. It is tuned with twin experiments.

Inflation. Even with localisation, the ensemble is usually underdispersive: its spread is smaller than its real error. There are two reasons. First, sampling error in $K$: the same members estimate the gain and are updated with it. Second, model error that the ensemble does not represent. An overconfident ensemble gives too little weight to new observations, shrinks further, and drifts away from the truth. This is filter divergence. Multiplicative inflation

with $\lambda$ about 1.01 to 1.1, multiplies $P^f$ by $\lambda^2$. Additive inflation adds random perturbations, and relaxation methods (RTPS, RTPP) pull the analysis spread back toward the forecast spread. Inflation changes the size of $P^f$, not its span and not its correlations. So inflation alone cannot fix the rank problem or the spurious correlations. That is the job of localisation.

Interactive figure: loading…

Square-root filters: ETKF and LETKF

A square-root (deterministic) filter updates the mean with the Kalman formula and transforms the anomalies so that their sample covariance equals $(I - KH)P^f$ exactly, with no random numbers. Let $X = A/\sqrt{N - 1}$ (so $P^f = XX^\top$) and $Y = HX$. The ensemble transform Kalman filter (ETKF; Hunt, Kostelich and Szunyogh 2007) works in the $N$-dimensional space of ensemble weights:

$$\begin{aligned} \tilde P^a &= \big(I_N + Y^\top R^{-1} Y\big)^{-1}, \\ \bar x^a &= \bar x + X\tilde P^aY^\top R^{-1}\big(y - H\bar x\big), \\ X^a &= X\,(\tilde P^a)^{1/2}, \end{aligned}$$

with the symmetric square root. Then $X^a(X^a)^\top = X\tilde P^aX^\top = (I - KH)P^f$, and the new anomalies still sum to zero. The cost is dominated by $N \times N$ matrices, so it is cheap. The LETKF does this for every grid point with only the nearby observations (domain localisation), and the points are independent, so it runs in parallel. Several weather services use LETKF-type systems for their ensembles (for example the KENDA system: Schraff et al. 2016).

Variational DA: 3D-Var and 4D-Var as MAP estimation #

3D-Var. Prior $\mathcal N(x^b, B)$, observation $y = h(x) + v$ with $v \sim \mathcal N(0, R)$, and a possibly nonlinear observation operator $h$. Bayes' rule gives $-\log p(x \mid y) = J(x) + \text{const}$ with

The first term is the background term $J_b$, the second the observation term $J_o$. 3D-Var computes the MAP estimate $\arg\min_x J(x)$ with a quasi-Newton or conjugate-gradient method. It is Tikhonov regularization with weighted norms (lesson 04).

For linear $h(x) = Hx$, the 3D-Var minimum is the Kalman analysis mean. $J$ is quadratic. Set its gradient to zero: $\nabla J = B^{-1}(x - x^b) - H^\top R^{-1}(y - Hx) = 0$, so $(B^{-1} + H^\top R^{-1}H)\,x = B^{-1}x^b + H^\top R^{-1}y$. This is the information form of the Kalman mean with $x^f = x^b$, $P^f = B$. The Hessian $B^{-1} + H^\top R^{-1}H$ is $(P^a)^{-1}$. (For a Gaussian posterior the mode and the mean are the same point.)

Strong-constraint 4D-Var. Use all observations of a time window, and assume a perfect model: $x_{j+1} = \mathcal M(x_j)$, where $\mathcal M$ is one (possibly nonlinear) time step and $j = 0, \dots, s$ counts the steps. Then every state $x_j$ is a function of the initial state $x_0$, and $x_0$ is the only control variable. Observations $y_j$ exist at the steps $j$ in a set $\mathcal O$:

$J$ is $-\log p(x_0 \mid \text{all } y_j)$ up to a constant, so $\arg\min J$ is the MAP initial state. The forecast then starts from this analysis.

Theorem 4 (4D-Var gradient by the adjoint). Let $M_j = \partial\mathcal M/\partial x$ at $x_j$ (the tangent-linear model) and $H_j$ the Jacobian of $h$ at $x_j$. Set $f_j = -H_j^\top R^{-1}\big(y_j - h(x_j)\big)$ for $j \in \mathcal O$, and $f_j = 0$ otherwise. Run the adjoint recursion backward in time, with the adjoint state $\lambda_j \in \R^n$:

for $j = s - 1, \dots, 0$.

Then $\nabla J(x_0) = B^{-1}(x_0 - x^b) + \lambda_0$.

Proof.

1. Write the observation term as $J_o = \sum_{j \in \mathcal O} \phi_j(x_j)$ with $\phi_j(x) = \tfrac12 (y_j - h(x))^\top R^{-1}(y_j - h(x))$. Then $\nabla\phi_j(x_j) = f_j$.
2. $x_j$ depends on $x_0$ through $j$ model steps. By the chain rule, its Jacobian is $M_{j-1}\cdots M_1M_0$. So $\nabla_{x_0}\,\phi_j(x_j) = M_0^\top M_1^\top \cdots M_{j-1}^\top f_j$.
3. Let $\lambda_j$ be the gradient, with respect to $x_j$, of all observation terms at steps $\ge j$. The terms at steps $> j$ depend on $x_j$ only through $x_{j+1} = \mathcal M(x_j)$, so their gradient is $M_j^\top\lambda_{j+1}$. The term at step $j$ adds $f_j$. This is the recursion. At $j = 0$ it gives $\lambda_0 = \nabla_{x_0} J_o$. The recursion shares the products of step 2 between all observation times.
4. Add the gradient of the background term, $B^{-1}(x_0 - x^b)$. $\blacksquare$

Cost. One forward run (store $x_0, \dots, x_s$, or checkpoints as in lesson 03) and one backward run of the adjoint model. The adjoint run costs a small constant times a forward run, independent of $n$. Finite differences need $n + 1$ forward runs. In ML terms, each step $\lambda_j = M_j^\top\lambda_{j+1}$ is a vector–Jacobian product of reverse-mode AD. If $\mathcal M$ is one RK4 step, $M_j^\top\lambda$ is the VJP of that step (the discrete adjoint of lesson 03). It is the exact gradient of the discrete cost, so it agrees with finite differences up to their own error.

Interactive figure: loading…

Practical notes.

• Linear model, no model error ($Q = 0$): $J$ is quadratic, and the 4D-Var analysis, moved to the end of the window, equals the KF analysis there. Both are the same Gaussian posterior (worked example (c)).
• A 4D-Var analysis has no error bar by itself. It is the mode. The uncertainty needs the Hessian of $J$ (its inverse is the posterior covariance in the linear case) or an ensemble.
• Window length. In a chaotic model, small changes of $x_0$ grow like $e^{\lambda_1 t}$, where $\lambda_1$ is the largest Lyapunov exponent (lesson 05). So $J(x_0)$ gets more and more local minima as the window grows. Operational windows are 6 to 12 hours; ERA5 uses 12-hour windows (Hersbach et al. 2020).

Incremental and weak-constraint 4D-Var

Incremental 4D-Var (Courtier, Thépaut and Hollingsworth 1994; used at ECMWF) is a Gauss–Newton method. Each outer loop runs the full nonlinear model from the current $x_0$. Each inner loop minimises a quadratic cost for the increment $\delta x_0$, built with the tangent-linear model $M_j$ and its adjoint $M_j^\top$, often at a lower resolution. The quadratic inner problem is solved with conjugate gradients, preconditioned by $B^{1/2}$ (the control variable is $B^{-1/2}\delta x_0$).

Weak-constraint 4D-Var drops the perfect-model assumption. It adds model-error terms $\tfrac12\sum_j \eta_j^\top Q^{-1}\eta_j$ with $x_{j+1} = \mathcal M(x_j) + \eta_j$, so the control variable is the whole trajectory, not only $x_0$. This allows longer windows.

Ensembles inside 4D-Var. The background covariance $B$ is hard to specify. Modern systems estimate part of it from an ensemble of data assimilations (ERA5 uses 10 members; Hersbach et al. 2020), or combine 4D-Var with an ensemble covariance ("hybrid" and "4DEnVar" methods).

Answers to section 1. (1) $K = 4/5$, so $x^a = 0.8 \cdot 5 = 4$ and $P^a = 0.2 \cdot 4 = 0.8$. With $y = 50$ the mean is $40$, but the variance is still $0.8$: in the linear-Gaussian case the posterior variance does not depend on the data. (A real system would also check whether $y = 50$ is a gross error before using it.) (2) The estimate of $x_2$ also increases, by $K_2(y - x^f_1)$ with $K_2 = P^f_{21}/S > 0$. Its variance also falls. (3) Finite differences: $n + 1$ runs (forward differences) or $2n$ (central), so about $10^8$. Adjoint: one forward run plus one adjoint run, whatever $n$ is.

4. Worked example #

(a) A Kalman update in two variables. Two neighbouring grid points have forecast temperature anomalies $x^f = (0, 0)$ with $P^f = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$ (correlation $1/\sqrt{2 \cdot 2} = 0.5$). A station observes point 1 only: $H = (1\;\;0)$, $R = 1$, $y = 3$.

1. $HP^f = (2\;\;1)$ (the first row of $P^f$), and $HP^fH^\top = 2$. So $S = 2 + 1 = 3$.
2. $P^fH^\top = (2, 1)^\top$ (the first column of $P^f$). So $K = (2, 1)^\top/3 = (2/3,\ 1/3)^\top$.
3. The innovation is $y - Hx^f = 3 - 0 = 3$. So $x^a = x^f + 3K = (2, 1)$. The unobserved point 2 moves by $1$, through the correlation.
4. $KH = \begin{pmatrix} 2/3 & 0 \\ 1/3 & 0 \end{pmatrix}$, so $I - KH = \begin{pmatrix} 1/3 & 0 \\ -1/3 & 1 \end{pmatrix}$. Row 1 of $P^a = (I - KH)P^f$ is $\tfrac13 (2, 1) = (2/3, 1/3)$. Row 2 is $-\tfrac13(2, 1) + (1, 2) = (1/3, 5/3)$. So $P^a = \begin{pmatrix} 2/3 & 1/3 \\ 1/3 & 5/3 \end{pmatrix}$. The variance at point 1 drops from 2 to $2/3$; at point 2 from 2 to $5/3$.
5. Check with the information form. $(P^f)^{-1} = \frac13\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}$ and $H^\top R^{-1}H = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}$, so their sum is $\begin{pmatrix} 5/3 & -1/3 \\ -1/3 & 2/3 \end{pmatrix}$. Its determinant is $10/9 - 1/9 = 1$, and its inverse is $\begin{pmatrix} 2/3 & 1/3 \\ 1/3 & 5/3 \end{pmatrix} = P^a$. ✓
6. 3D-Var with $x^b = x^f$ and $B = P^f$ has the Hessian $(P^f)^{-1} + H^\top R^{-1}H = (P^a)^{-1}$, and its gradient at $(2, 1)$ is $B^{-1}(2, 1)^\top - H^\top(3 - 2) = (1, 0)^\top - (1, 0)^\top = 0$. So the 3D-Var minimum is $(2, 1) = x^a$. ✓ (Here $B^{-1}(2, 1)^\top = \frac13(4 - 1,\ -2 + 2)^\top = (1, 0)^\top$.)

(b) Why the stochastic EnKF perturbs the observations (one variable, $H = 1$). Let $P^f = 1$ and $R = 1$, so $K = 1/2$. A forecast anomaly $a$ becomes $(1 - K)a + K\varepsilon'$.

• With perturbed observations, the expected analysis variance is $(1 - K)^2P^f + K^2R = 0.25 + 0.25 = 0.5$. This equals $(1 - K)P^f = 0.5$, the KF value. ✓
• Without them it is $(1 - K)^2P^f = 0.25$: half the correct variance. Repeat this every cycle and the ensemble becomes far too confident.

(c) 4D-Var equals the KF at the end of the window (one variable). Model $x_{j+1} = 2x_j$ with no model error. Background $x^b = 0$, $B = 1$. One observation $y_1 = 2$ at step 1 with $R = 1$ and $h(x) = x$.

1. $x_1 = 2x_0$, so $J(x_0) = \tfrac12 x_0^2 + \tfrac12 (2 - 2x_0)^2$.
2. Adjoint: $f_1 = -(y_1 - x_1) = -(2 - 2x_0)$, $\lambda_1 = f_1$, and $\lambda_0 = M_0^\top\lambda_1 + f_0 = 2\lambda_1 + 0$. So $\nabla J = x_0 + 2(2x_0 - 2) = 5x_0 - 4$.
3. Check at $x_0 = 0.5$: the adjoint gives $\nabla J = 2.5 - 4 = -1.5$. Direct differentiation of $J$ gives $x_0 - 2(2 - 2x_0) = 0.5 - 2 = -1.5$. ✓
4. Minimum: $5x_0 - 4 = 0$ gives $x_0^\ast = 0.8$. The Hessian is $5$, so the posterior variance of $x_0$ is $1/5 = 0.2$.
5. Move it to step 1: $x_1 = 2 \cdot 0.8 = 1.6$, with variance $2^2 \cdot 0.2 = 0.8$.
6. Kalman filter: forecast $x^f_1 = 2 \cdot 0 = 0$, $P^f_1 = 2^2 \cdot 1 + 0 = 4$. Then $K = 4/5$, $x^a_1 = 0 + 0.8 \cdot 2 = 1.6$, $P^a_1 = 0.2 \cdot 4 = 0.8$. The same. ✓

5. In code #

Three short NumPy checks. First, the Kalman formulas and Bayes' rule on a grid (no Kalman formula) give the same answer for worked example (a):

import numpy as np
from scipy.stats import multivariate_normal, norm

# Worked example (a): prior N(x_f, P_f) on (x1, x2); one observation y = x1 + v, v ~ N(0, 1).
x_f, P_f, y = np.zeros(2), np.array([[2.0, 1.0], [1.0, 2.0]]), 3.0

# 1) The Kalman formulas.
H, R = np.array([[1.0, 0.0]]), np.array([[1.0]])
S = H @ P_f @ H.T + R                                # innovation covariance, here 3
K = np.linalg.solve(S, H @ P_f).T                    # K = P_f H^T S^-1 (a solve, not inv)
x_a = x_f + K @ (y - H @ x_f)
P_a = (np.eye(2) - K @ H) @ P_f
print(K.ravel().round(4), x_a.round(4))              # [0.6667 0.3333] [2. 1.]

# 2) Bayes' rule on a grid, with no Kalman formula: posterior ∝ prior × likelihood.
g = np.linspace(-8.0, 10.0, 901)
X1, X2 = np.meshgrid(g, g, indexing="ij")
prior = multivariate_normal(x_f, P_f).pdf(np.stack([X1, X2], axis=-1))
likelihood = norm.pdf(y, loc=X1, scale=1.0)           # y = x1 + v, v ~ N(0, 1)
post = prior * likelihood
post /= post.sum()
mean = np.array([(post * X1).sum(), (post * X2).sum()])
D = np.stack([X1 - mean[0], X2 - mean[1]])            # shape (2, 901, 901)
cov = np.einsum("aij,bij,ij->ab", D, D, post)
print(mean.round(4))                                  # [2. 1.]
print(cov.round(4))                                   # [[0.6667 0.3333]
                                                      #  [0.3333 1.6667]]
print(np.allclose(cov, P_a, atol=1e-4))               # True

Second, sampling error and localisation. A 10-member ensemble estimates a covariance with correlation length 3 on 60 grid points:

import numpy as np

def gaspari_cohn(r):
    """Gaspari-Cohn taper; r = distance / c. Equal to 1 at r = 0 and exactly 0 for r >= 2."""
    r = np.abs(r)
    out = np.zeros_like(r)
    a, b = r <= 1, (r > 1) & (r < 2)
    ra, rb = r[a], r[b]
    out[a] = -0.25 * ra**5 + 0.5 * ra**4 + 0.625 * ra**3 - 5 / 3 * ra**2 + 1
    out[b] = rb**5 / 12 - 0.5 * rb**4 + 0.625 * rb**3 + 5 / 3 * rb**2 - 5 * rb + 4 - 2 / (3 * rb)
    return out

rng = np.random.default_rng(0)
n, N, ell = 60, 10, 3.0                           # grid points on a line, members, correlation length
d = np.abs(np.arange(n)[:, None] - np.arange(n)[None, :]).astype(float)
B = np.exp(-d / ell)                              # true forecast-error covariance
rho = gaspari_cohn(d / (2 * ell))                 # taper with c = 2 ell: zero beyond distance 4 ell
print(np.linalg.eigvalsh(rho).min() > 0)          # True: the taper is a correlation matrix

err_raw, err_loc, far_corr = [], [], []
for _ in range(200):
    X = rng.multivariate_normal(np.zeros(n), B, size=N)   # ensemble, one member per row
    Pf = np.cov(X, rowvar=False)                           # sample covariance, divides by N - 1
    err_raw.append(np.linalg.norm(Pf - B) / np.linalg.norm(B))
    err_loc.append(np.linalg.norm(rho * Pf - B) / np.linalg.norm(B))
    far_corr.append(np.corrcoef(X, rowvar=False)[d >= 30])  # true correlation <= e^-10
print(np.linalg.matrix_rank(Pf), np.linalg.matrix_rank(rho * Pf))   # 9 60
print(round(np.std(np.concatenate(far_corr)), 3), round(1 / np.sqrt(N - 1), 3))   # 0.334 0.333
print(round(np.mean(err_raw), 2), round(np.mean(err_loc), 2))       # 1.5 0.59

The raw sample covariance has rank $N - 1 = 9$. Its relative error (1.5) is larger than the error of the zero matrix (1.0). The far correlations are pure noise of size $1/\sqrt{N - 1}$. The tapered covariance has full rank and less than half the error.

Third, 4D-Var on a linear oscillator (position and velocity; only the position is observed). cost_and_grad runs exactly the adjoint recursion of Theorem 4, and the minimiser equals the solution of one weighted least-squares problem:

import numpy as np
from scipy.linalg import expm
from scipy.optimize import minimize

# Linear oscillator: state (position, velocity); exact one-step propagator M; observe position only.
h, n_steps, every, sig_o = 0.1, 50, 10, 0.1
M = expm(h * np.array([[0.0, 1.0], [-1.0, 0.0]]))
H = np.array([[1.0, 0.0]])
x_b, B_inv, R_inv = np.array([0.5, -0.5]), np.eye(2), np.eye(1) / sig_o**2
obs_steps = range(0, n_steps + 1, every)             # observations at steps 0, 10, ..., 50

rng = np.random.default_rng(1)
x, ys = np.array([1.0, 0.5]), {}                     # the truth starts at (1, 0.5)
for j in range(n_steps + 1):
    if j in obs_steps:
        ys[j] = H @ x + sig_o * rng.standard_normal(1)
    x = M @ x

def cost_and_grad(x0):
    """Strong-constraint 4D-Var: one forward run, then one backward (adjoint) run."""
    xs = [x0]
    for j in range(n_steps):                         # forward: store the trajectory
        xs.append(M @ xs[-1])
    J = 0.5 * (x0 - x_b) @ B_inv @ (x0 - x_b)
    lam = np.zeros(2)
    for j in range(n_steps, -1, -1):                 # backward: lam_j = M^T lam_{j+1} + f_j
        if j < n_steps:
            lam = M.T @ lam
        if j in ys:
            d = ys[j] - H @ xs[j]                    # innovation at step j
            J += 0.5 * d @ R_inv @ d
            lam = lam - H.T @ R_inv @ d              # f_j = -H^T R^-1 (y_j - H x_j)
    return J, B_inv @ (x0 - x_b) + lam

v, e = np.array([0.2, 0.3]), 1e-6
fd = [(cost_and_grad(v + e * u)[0] - cost_and_grad(v - e * u)[0]) / (2 * e) for u in np.eye(2)]
print(np.abs(cost_and_grad(v)[1] - fd).max() < 1e-5)   # True: adjoint gradient = finite differences

x_4dvar = minimize(cost_and_grad, x_b, jac=True, method="L-BFGS-B").x
# The same MAP estimate as one weighted least-squares problem: y_j = H M^j x0 + noise.
G = np.vstack([H @ np.linalg.matrix_power(M, j) for j in obs_steps])
A = np.vstack([np.eye(2), G / sig_o])                # whitened rows (B = I)
b = np.concatenate([x_b, np.concatenate([ys[j] for j in obs_steps]) / sig_o])
x_lsq = np.linalg.lstsq(A, b, rcond=None)[0]
print(x_4dvar.round(3), x_lsq.round(3))              # [1.049 0.484] [1.049 0.484]

The velocity is never observed, but 4D-Var recovers it (0.484; the truth is 0.5). The model carries the information from the observed position to the velocity.

6. Where it shows up #

• Global weather prediction. ECMWF uses incremental 4D-Var together with an ensemble of data assimilations that estimates the background errors (see the IFS documentation, Part II). ERA5 was made with this system: 12-hour 4D-Var windows and a 10-member ensemble of 4D-Vars (Hersbach et al. 2020).
• Kilometre-scale ensembles. The KENDA system (Schraff et al. 2016) is an LETKF, an ensemble Kalman filter with domain localisation. It was built for the COSMO model and is used by national weather services for high-resolution ensemble forecasts.
• ML for DA (as of 2026). An ML emulator makes large ensembles cheap, and automatic differentiation gives its adjoint for 4D-Var (test it against finite differences, as here). Score-based data assimilation (Rozet and Louppe 2023) uses a diffusion-model prior over whole trajectories and adds the observations only at sampling time (lesson 13). End-to-end systems such as Aardvark Weather (Allen et al. 2025) learn the whole path from raw observations to forecasts.
• ML in general. The KF and the RTS smoother are exact inference in linear dynamical systems. 4D-Var is backpropagation through time with respect to an initial condition. Ensemble Kalman inversion uses the EnKF update as a derivative-free optimiser, also for training neural networks (Kovachki and Stuart 2019).
• Seismology. Full-waveform inversion uses the same adjoint method as 4D-Var: a forward wavefield, a backward adjoint wavefield, and the gradient with respect to the earth model (lesson 04). The control is a parameter field, not an initial state. Ensemble data assimilation has also been used to estimate and forecast the stresses on a fault in models of earthquake sequences (van Dinther, Künsch and Fichtner 2019).

7. Common mistakes #

• Using inv(S) instead of a solve with $S$, and letting $(I - KH)P$ drift away from symmetric positive semidefinite over many cycles. Symmetrise, or use the Joseph form.
• Mixing up the time indices. $x^f_k$ is the forecast for time $k$ made from $x^a_{k-1}$; the observation $y_k$ updates $x^f_k$ into $x^a_k$. Updating with $y_{k-1}$ or forecasting from $x^f_{k-1}$ is a common bug.
• A stochastic EnKF without perturbed observations, or one perturbation shared by all members. The spread becomes too small, and the filter diverges.
• $1/N$ instead of $1/(N - 1)$ in the sample covariance. With $N = 10$ this is a 10% error.
• Expecting inflation to fix $N \ll n$. Inflation keeps the span of the anomalies. Only localisation (or more members) adds directions and removes spurious correlations.
• Reading a 3D-Var or 4D-Var analysis as having no uncertainty. These methods give the mode. The uncertainty needs the Hessian (its inverse is $P^a$ in the linear case) or an ensemble.
• 4D-Var adjoint bugs: linearising around the wrong state (use the stored forward trajectory), forgetting the observation at the initial time $j = 0$, or using the continuous adjoint and expecting it to match finite differences of the discrete cost to round-off (lesson 03).

8. Check yourself #

1. Prior $\mathcal N(10, 4)$, observation $y = 13$ with $R = 2$, $H = 1$. Give $K$, $x^a$ and $P^a$.
   Answer

   $K = 4/6 = 2/3$, $x^a = 10 + \tfrac23 \cdot 3 = 12$, $P^a = \tfrac13 \cdot 4 = 4/3$. Check with precisions: $1/4 + 1/2 = 3/4 = 1/P^a$.
2. Show that $P^a \preceq P^f$ in Theorem 2: an observation never increases the uncertainty.
   Answer

   $P^f - P^a = P^fH^\top S^{-1}HP^f = KSK^\top$, because $KSK^\top = P^fH^\top S^{-1}SS^{-1}HP^f$. Since $S \succ 0$, $z^\top KSK^\top z = (K^\top z)^\top S\,(K^\top z) \ge 0$ for every $z$.
3. Your EnKF has $N = 40$ members for a model with $n = 10^6$ variables. Why can the analysis increment lie in only 39 directions, and why does localisation help where inflation does not?
   Answer

   $K = P^fH^\top S^{-1}$ with $P^f = AA^\top/(N-1)$, so $K = AC$: every increment is a combination of the anomalies, which sum to zero (rank $\le N - 1 = 39$). Inflation multiplies the anomalies by $\lambda$: same span, same correlations. Localisation $\rho \circ P^f$ removes the spurious long-range correlations (of size about $1/\sqrt{N-1} \approx 0.16$) and gives a covariance that usually has full rank, so the analysis can act locally in many more directions.
4. 4D-Var with the scalar model $x_{j+1} = 3x_j$, $x^b = 1$, $B = 1$, and one observation $y = 0$ at step 1 with $R = 9$, $h(x) = x$. (a) Compute $\nabla J(1)$ with the adjoint recursion. (b) Find the analysis $x_0^\ast$ and its posterior variance.
   Answer

   (a) $x_1 = 3$, $f_1 = -(0 - 3)/9 = 1/3 = \lambda_1$, $\lambda_0 = 3\lambda_1 = 1$, so $\nabla J(1) = (1 - 1) + 1 = 1$. (b) $J = \tfrac12(x_0 - 1)^2 + \tfrac12 (3x_0)^2/9 = \tfrac12(x_0 - 1)^2 + \tfrac12 x_0^2$, so $J' = 2x_0 - 1$ (equal to 1 at $x_0 = 1$ ✓). Then $x_0^\ast = 1/2$, and the variance is $1/J'' = 1/2$.
5. For a nonlinear model, why is the EnKF forecast step often better than the extended Kalman filter (EKF), which propagates $P$ with the tangent-linear model?
   Answer

   The EnKF runs the full nonlinear model for each member. It needs no tangent-linear or adjoint code, it needs $N$ runs instead of about $n$, and the ensemble feels the nonlinear evolution of the mean and the spread. The EKF needs about $n$ tangent-linear runs to propagate $P$, and its linearisation error grows with the forecast length.
6. Why does strong-constraint 4D-Var become hard when the window is long in a chaotic model?
   Answer

   Perturbations of $x_0$ grow like $e^{\lambda_1 t}$, where $\lambda_1$ is the largest Lyapunov exponent. So $J(x_0)$ becomes very rough, with many local minima, and the tangent-linear approximation behind incremental 4D-Var fails. A gradient method then finds a local minimum. Remedies: shorter windows (6 to 12 hours), growing the window step by step, weak-constraint 4D-Var, or ensemble methods.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. Give the Kalman gain for prior covariance $P^f$, observation operator $H$ and observation-error covariance $R$.

   Show answer

2. Give the analysis mean and covariance of the Kalman update.

   Show answer

3. Scalar case: $P^f = 4$, $R = 1$, $x^f = 0$, $y = 5$. Give $K$, $x^a$ and $P^a$.

   Show answer

4. Write the forecast step of the Kalman filter.

   Show answer

5. Why does observing $x_1$ change the estimate of the unobserved $x_2$?

   Show answer

6. Why does the stochastic EnKF perturb the observations?

   Show answer

7. What is the largest rank of the EnKF sample covariance with $N$ members, and why?

   Show answer

8. What does localisation do in an EnKF?

   Show answer

9. What does multiplicative inflation fix, and what does it not fix?

   Show answer

10. How is the 4D-Var gradient computed, and what does it cost?

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Derive the Kalman gain." Gaussian conditioning in three lines (Theorem 2). Another route: the $K$ that minimises the trace of the Joseph-form covariance $(I - KH)P(I - KH)^\top + KRK^\top$.
• "KF, EKF, EnKF, particle filter, 4D-Var: when do you use which?" KF: linear-Gaussian, small $n$. EKF: mildly nonlinear, small $n$. EnKF: large nonlinear models, near-Gaussian errors, no adjoint needed. Particle filter: strongly non-Gaussian, small $n$. 4D-Var: large models with an adjoint; it uses the dynamics over a window and gives a mode without an error bar.
• "What is 4D-Var in ML language?" A MAP fit of the initial condition of a differentiable simulator. The adjoint is backpropagation through time; memory is handled by checkpointing.
• "Your EnKF spread is much smaller than its RMSE. What do you check?" Perturbed observations (or a correct square-root update), the localisation radius, inflation, model error, and the ensemble size.
• "How would you use an ML emulator in DA?" For cheap ensembles in an EnKF, or for its autodiff adjoint in 4D-Var. Treat its error as model error, and test its gradients against finite differences.

10. Practice #

The notebook for this lesson has five exercises, each with a check cell and a folded solution:

1. Write the Kalman update and check it against Bayes' rule on a grid for a two-variable Gaussian.
2. Run the Kalman filter on a noisy oscillator with only the position observed, and check it against the exact posterior from one large Gaussian conditioning.
3. Write the 3D-Var cost and gradient, minimise it, and show that the minimum is the Kalman analysis.
4. Write the stochastic EnKF analysis with inflation, and run a twin experiment on Lorenz-63: the error must fall below the free run and below the observation error.
5. Write the strong-constraint 4D-Var gradient with the adjoint recursion for Lorenz-63 and check it with finite differences.

11. Watch and read #

Visually Explained: Kalman FiltersVisually Explained · 11 min

The forecast and update steps as Gaussian pictures; a good first view before section 2.

The Kalman Filter [Control Bootcamp]Steve Brunton · 6 min

The Kalman filter as an optimal state estimator from a control point of view.

Data Assimilation: variational data assimilation and the ensemble Kalman filterEuropean Space Agency, ESA · 56 min

A lecture on 3D-Var, 4D-Var and the EnKF for Earth-system models; it matches sections 3 and 6.

[Data Assimilation] L5: Optimal interpolation and variational (3D/4D) methodsxwu · 89 min

A full university lecture; watch the parts on 3D-Var and 4D-Var and the equivalence with the Kalman analysis.

• Evensen, Vossepoel and van Leeuwen, Data Assimilation Fundamentals (Springer 2022, open access): ch. 2–3 (DA as Bayesian inference and MAP), ch. 4 (strong-constraint 4D-Var), ch. 6 (Kalman filters and 3D-Var), ch. 8 (ensemble methods), ch. 10 (localisation and inflation).
• Carrassi, Bocquet, Bertino and Evensen, Data assimilation in the geosciences: an overview of methods, issues and perspectives (WIREs Climate Change 2018): §3 (Kalman filter and 4D-Var), §4 (EnKF, square-root filters, localisation and inflation), §5 (chaotic dynamics).
• Särkkä and Svensson, Bayesian Filtering and Smoothing, 2nd ed. (2023, free PDF from the authors): the Kalman filter, the RTS smoother and the Gaussian identities.
• Nansen Environmental and Remote Sensing Center, DA-tutorials (Jupyter notebooks): Bayesian inference, the multivariate Kalman filter, and writing your own EnKF.
• Papers cited in this lesson: Hersbach et al. 2020, "The ERA5 global reanalysis" (Q. J. R. Meteorol. Soc.); Schraff et al. 2016, "Kilometre-scale ensemble data assimilation for the COSMO model (KENDA)" (Q. J. R. Meteorol. Soc.); Hunt, Kostelich and Szunyogh 2007, "Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter" (Physica D); Rozet and Louppe 2023, Score-based data assimilation; Allen et al. 2025, Aardvark weather; Kovachki and Stuart 2019, Ensemble Kalman inversion; Lam et al. 2023, GraphCast; van Dinther, Künsch and Fichtner 2019, "Ensemble data assimilation for earthquake sequences" (Geophys. J. Int.).