Data assimilation: 3DVAR, 4DVAR, and ensemble filters

For each integer $N\ge 2$, let $(X_{N,i},U_{N,i}{)}_{i=1}^N$ be a jointly exchangeable row of real random pairs: its law is invariant when the same permutation is applied to every index. Assume that $U_{N,1},\dots ,U_{N,N}$ are pairwise uncorrelated, all have the law of a fixed square-integrable random variable $U$, and $\mathrm{\parallel }{X}_{N,1}-U_{N,1}{\mathrm{\parallel }}_2:=(\mathbb{E}\mathrm{\mid }{X}_{N,1}-U_{N,1}{\mathrm{\mid }}^2{)}^{1\mathrm{/}2}⟶0.$ Here $\mathrm{Cov}(V,W):=\mathbb{E}[(V-\mathbb{E}V)(W-\mathbb{E}W)]$. Then the empirical mean ${\bar{X}}_N=N^{-1}{\sum }_{i=1}^N{X}_{N,i}$ converges to $\mathbb{E}U$ in $L^2$, meaning $\mathbb{E}\mathrm{\mid }{\bar{X}}_N-\mathbb{E}U{\mathrm{\mid }}^2\to 0$, and hence in probability.

Let $1\le q<s<\mathrm{\infty }$. Suppose random vectors $Z_n,Z$ in a fixed finite-dimensional Euclidean space satisfy $Z_n\to Z$ in probability, $\mathbb{E}\mathrm{\parallel }Z{\mathrm{\parallel }}^s<\mathrm{\infty }$, and ${\sup }_n\mathbb{E}\mathrm{\parallel }{Z}_n{\mathrm{\parallel }}^s<\mathrm{\infty }$. Then $\mathbb{E}\mathrm{\parallel }{Z}_n-Z{\mathrm{\parallel }}^q\to 0$.

Let $d\ge 1$. For each $N\ge 2$, let $(X_{N,i},U_{N,i}{)}_{i=1}^N$ be jointly exchangeable random pairs in ${\mathbb{R}}^d$. Assume that $U_{N,1},\dots ,U_{N,N}$ are iid with the law of a fixed $U$, $\mathbb{E}\mathrm{\parallel }U{\mathrm{\parallel }}^4<\mathrm{\infty }$, and $\mathrm{\parallel }{X}_{N,1}-U_{N,1}{\mathrm{\parallel }}_4\to 0$. Define ${\bar{X}}_N=N^{-1}{\sum }_i{X}_{N,i},C_N=N^{-1}{\sum }_i(X_{N,i}-{\bar{X}}_N)(X_{N,i}-{\bar{X}}_N{)}^T\mathrm{.}$ Here $\mathrm{Cov}(U):=\mathbb{E}[(U-\mathbb{E}U)(U-\mathbb{E}U{)}^T]$. Then ${\bar{X}}_N\to \mathbb{E}U$ and $C_N\to \mathrm{Cov}(U)$ in probability, entry by entry.

Fix a finite horizon $T\ge 0$, state dimensions $d_t\ge 1$, observation dimensions $p_t\ge 1$ for $1\le t\le T$, and forecast maps ${\mathrm{\Psi }}_{t-1}:{\mathbb{R}}^{d_{t-1}}\to {\mathbb{R}}^{d_t}$ and observation maps $G_t:{\mathbb{R}}^{d_t}\to {\mathbb{R}}^{p_t}$ that are measurable for the Borel sigma-algebras: inverse images of Borel sets are Borel (A measurable function between measurable spaces, The Borel sigma-algebra of a topological space). A deterministic model trajectory satisfies $x_t={\mathrm{\Psi }}_{t-1}(x_{t-1})$; an additive-error model instead has $X_t={\mathrm{\Psi }}_{t-1}(X_{t-1})+q_t$. Observations have the form $Y_t=G_t(X_t)+r_t$. A forecast at time $t$ is information about $X_t$ formed from the model and observations through time $t-1$, before using $Y_t$; an analysis at time $t$ is information after using $Y_t$. In the state-space terminology of Filtering, predictive, and smoothing laws in a finite-horizon state-space model, a filtering law is the conditional law of $X_t$ given $Y_{1:t}$ and a predictive law is the conditional law of $X_t$ given $Y_{1:t-1}$. A fixed-interval smoothing law conditions a state, or the full state path, on observations through the terminal time, as in Fixed-interval smoothing law in a linear-Gaussian model. Under that page's standard-Borel hypotheses, the filtering, predictive, and full-path smoothing conditional kernels exist by Filtering and finite-horizon smoothing laws exist; conditional versions may differ on null observation histories.

For a symmetric positive-definite matrix $M$, write $\mathrm{\parallel }v{\mathrm{\parallel }}_{M^{-1}}^2:=v^T{M}^{-1}v$. In a finite-dimensional Gaussian 3DVAR update, the forecast mean is $m^f$, its fixed background covariance is $B\succ 0$, the linear observation map is $H$, the datum is $y$, and the observation-error covariance is $\mathrm{\Gamma }\succ 0$. The objective is $J_{3D}(u)=\frac{1}{2}\mathrm{\parallel }u-m^f{\mathrm{\parallel }}_{B^{-1}}^2+\frac{1}{2}\mathrm{\parallel }y-Hu{\mathrm{\parallel }}_{{\mathrm{\Gamma }}^{-1}}^2\mathrm{.}$ It returns one analysis point; the fixed background covariance is not an analysis posterior covariance.

Fix an ensemble size $N\ge 2$, forecast map $\mathrm{\Psi }$, observation map $G$, datum $y\in {\mathbb{R}}^p$, and observation covariance $\mathrm{\Gamma }\succ 0$. The input is a state ensemble $x^1,\dots ,x^N\in {\mathbb{R}}^d$, forecast-error draws ${\xi }^i$ that are independent across members and of the input ensemble, each with the caller-specified additive model-error law, and observation perturbations ${\eta }^i\sim N_p(0,\mathrm{\Gamma })$ under Multivariate normal law, including singular covariance. The observation perturbations are independent across members and independent of the forecast ensemble and forecast-error draws.

1. Forecast each member: ${\widehat{x}}^i=\mathrm{\Psi }(x^i)+{\xi }^i$.

2. Form predicted observations $z^i=G({\widehat{x}}^i)$, means $\bar{x}=N^{-1}{\sum }_i{\widehat{x}}^i$, $\bar{z}=N^{-1}{\sum }_i{z}^i$, and empirical covariance blocks $C_{xz}=N^{-1}{\sum }_i({\widehat{x}}^i-\bar{x})(z^i-\bar{z}{)}^T,C_{zz}=N^{-1}{\sum }_i(z^i-\bar{z})(z^i-\bar{z}{)}^T\mathrm{.}$ The divisor-$N$ convention is part of this algorithm.

3. Form $S_N=C_{zz}+\mathrm{\Gamma }$. Since $C_{zz}$ is positive semidefinite and $\mathrm{\Gamma }$ is positive definite, $S_N$ is positive definite. Factor it by Cholesky and solve $S_N{K}_N^T=C_{xz}^T$ by triangular solves.

4. Update every member using its perturbed datum: $x^{a,i}={\widehat{x}}^i+K_N(y+{\eta }^i-z^i)\mathrm{.}$

Return the updated ensemble and its empirical mean and covariance. When $G(x)=Hx$, $C_{xz}=C_x{H}^T$ and $C_{zz}=H{C}_x{H}^T$, giving the usual sample-covariance Kalman gain. This item specifies the perturbed-observation variant; deterministic square-root updates are different finite-ensemble procedures. Finite ensembles in nonlinear models are approximations and carry no general posterior-exactness claim.

Fix $T\ge 1$, positive dimensions $d_0,\dots ,d_T$ and $p_1,\dots ,p_T$, open sets $U_t\subseteq {\mathbb{R}}^{d_t}$, and continuously differentiable maps ${\mathrm{\Psi }}_{t-1}:U_{t-1}\to U_t$ and $G_t:U_t\to {\mathbb{R}}^{p_t}$. Let $m_b\in {\mathbb{R}}^{d_0}$, $B\in {\mathbb{R}}^{d_0\times d_0}$ be symmetric positive definite, $y_t\in {\mathbb{R}}^{p_t}$, and ${\mathrm{\Gamma }}_t\in {\mathbb{R}}^{p_t\times p_t}$ be symmetric positive definite. In weak-constraint mode also require symmetric positive-definite $Q_t\in {\mathbb{R}}^{d_t\times d_t}$. The strong-constraint decision domain is $U_0$; the weak-constraint decision domain is $U_0\times \cdots \times U_T$. Both objectives are continuously differentiable on their stated open decision domains. For a positive-definite covariance $M$, write $\mathrm{\parallel }v{\mathrm{\parallel }}_{M^{-1}}^2=v^T{M}^{-1}v$.

In strong-constraint 4DVAR, choose the initial state $u_0$ and generate $u_t={\mathrm{\Psi }}_{t-1}(u_{t-1})$ for $1\le t\le T$. Minimize the trajectory objective $J_s(u_0)=\frac{1}{2}\mathrm{\parallel }{u}_0-m_b{\mathrm{\parallel }}_{B^{-1}}^2+\frac{1}{2}{\sum }_{t=1}^T\mathrm{\parallel }{y}_t-G_t(u_t){\mathrm{\parallel }}_{{\mathrm{\Gamma }}_t^{-1}}^2\mathrm{.}$

In weak-constraint 4DVAR, the decision variable is the full path $(u_0,\dots ,u_T)$. Given model-error covariances $Q_t\succ 0$, minimize $J_w(u_0,\dots ,u_T)=\frac{1}{2}\mathrm{\parallel }{u}_0-m_b{\mathrm{\parallel }}_{B^{-1}}^2+\frac{1}{2}{\sum }_{t=1}^T\mathrm{\parallel }{y}_t-G_t(u_t){\mathrm{\parallel }}_{{\mathrm{\Gamma }}_t^{-1}}^2+\frac{1}{2}{\sum }_{t=1}^T\mathrm{\parallel }{u}_t-{\mathrm{\Psi }}_{t-1}(u_{t-1}){\mathrm{\parallel }}_{Q_t^{-1}}^2\mathrm{.}$

Select one objective, an initial decision vector in its stated decision domain, a finite iteration limit $M\ge 0$, and either positive fixed steps or the backtracking parameters required by Gradient descent. Run that algorithm on the selected objective for at most $M$ updates and return its last iterate (or report its domain failure). In strong mode return the trajectory generated from the last initial-state iterate; in weak mode return the last path iterate. This is a finite numerical procedure, not a guarantee that a global or local minimizer is reached. The page's BCS-20 hypotheses for gradient descent supply no general guarantee for these objectives.

For integers $d,p\ge 1$, given a forecast mean $m^f\in {\mathbb{R}}^d$, fixed background covariance $B\in {\mathbb{R}}^{d\times d}$, linear observation matrix $H\in {\mathbb{R}}^{p\times d}$, datum $y\in {\mathbb{R}}^p$, and observation-error covariance $\mathrm{\Gamma }\in {\mathbb{R}}^{p\times p}$, require $B$ and $\mathrm{\Gamma }$ symmetric positive definite. Define the 3DVAR objective $J(u)=\frac{1}{2}(u-m^f{)}^T{B}^{-1}(u-m^f)+\frac{1}{2}(y-Hu{)}^T{\mathrm{\Gamma }}^{-1}(y-Hu)\mathrm{.}$ Form $S=HB{H}^T+\mathrm{\Gamma }$. It is positive definite because for every $z\mathrm{\ne }0$, $z^{T}Sz=(H^{T}z{)}^{T}B(H^{T}z)+z^T\mathrm{\Gamma }z>0.$ Factor $S=L{L}^T$ with the unique positive-diagonal Cholesky factor from A matrix admits a Cholesky factorisation with positive diagonal exactly when it is Hermitian positive definite, and that factor is unique. Solve $LW=HB$ and $L^T{K}^T=W$; thus $K^T=S^{-1}HB$ without explicitly forming $S^{-1}$. Return the analysis point $m^a=m^f+K(y-H{m}^f)\mathrm{.}$ For a sequential 3DVAR filter, first set $m_{j+1}^f={\mathrm{\Psi }}_j(m_j^a)$, then apply this update with the fixed matrices $B,H,\mathrm{\Gamma }$ and observation $y_{j+1}$. The routine returns a point estimate only. Its unique-minimizer and Gaussian-MAP properties are established separately; no analysis covariance is returned.

Assume AC and the finite-horizon linear-Gaussian model of Linear-Gaussian state-space models with independent noises, with $R_k\succ 0$ for $1\le k\le T$. Fix the observation values $y_{1:T}$. For each $N$, initialize $U_{N,i}^{a,0}$ iid with the model's $N(m_0,P_0)$ prior. At each time use independent per-member noises $q_{N,i}^k\sim N(0,Q_k)$ and ${\eta }_{N,i}^k\sim N(0,R_k)$, independent across $i,k$ and from the preceding ensemble, and form $U_{N,i}^{f,k}=F_k{U}_{N,i}^{a,k-1}+q_{N,i}^k,U_{N,i}^{a,k}=U_{N,i}^{f,k}+K_k(y_k+{\eta }_{N,i}^k-H_k{U}_{N,i}^{f,k}),$ where $K_k=P_{k\mathrm{\mid }k-1}{H}_k^T(H_k{P}_{k\mathrm{\mid }k-1}{H}_k^T+R_k{)}^{-1}$ is the deterministic gain from Exact linear-Gaussian Kalman filter. Then, for every $k\le T$, the analysis members $U_{N,1}^{a,k},\dots ,U_{N,N}^{a,k}$ are iid with law $N(m_{k\mathrm{\mid }k},P_{k\mathrm{\mid }k})$, the exact Kalman filtering law at the fixed data $y_{1:k}$.

Fix the perturbed-observation update of Perturbed-observation ensemble Kalman filter. A row of random inputs is exchangeable if its joint law is invariant under every common permutation of the member indices. If the forecast states, forecast errors, observation perturbations, and any reference-member variables form an exchangeable input row, then the updated states together with those reference variables are exchangeable.

Fix a finite horizon linear-Gaussian model from Linear-Gaussian state-space models with independent noises, with each observation covariance $R_k\succ 0$. Run the perturbed-observation procedure of Perturbed-observation ensemble Kalman filter with iid Gaussian initial members and independent Gaussian process and observation perturbations having the model's laws. For every finite $p\ge 1$, there are constants $c_{k,p}$ for $0\le k\le T$ and $b_{k,p}$ for $1\le k\le T$, independent of ensemble size $N\ge 2$ and member index $i$, such that ${\sup }_{N,i}\mathbb{E}\mathrm{\parallel }{X}_{N,i}^{a,k}{\mathrm{\parallel }}^p\le c_{k,p}(0\le k\le T),{\sup }_N\mathbb{E}\mathrm{\parallel }{K}_N^k{\mathrm{\parallel }}_{\mathrm{o}\mathrm{p}}^p\le b_{k,p}(1\le k\le T)\mathrm{.}$ Here $\mathrm{\parallel }\cdot {\mathrm{\parallel }}_{\mathrm{o}\mathrm{p}}$ is the Euclidean operator norm. The claim is for each fixed finite horizon, not a bound uniform in $k$ as the horizon grows.

Let $d,p\ge 1$, $m^f\in {\mathbb{R}}^d$, $y\in {\mathbb{R}}^p$, $H\in {\mathbb{R}}^{p\times d}$, and symmetric positive-definite $B\in {\mathbb{R}}^{d\times d}$ and $\mathrm{\Gamma }\in {\mathbb{R}}^{p\times p}$. Put $S=HB{H}^T+\mathrm{\Gamma }$, $K=B{H}^T{S}^{-1}$, and $J(u)=\frac{1}{2}(u-m^f{)}^T{B}^{-1}(u-m^f)+\frac{1}{2}(y-Hu{)}^T{\mathrm{\Gamma }}^{-1}(y-Hu)\mathrm{.}$ The 3DVAR analysis $m^a=m^f+K(y-H{m}^f)$ from Finite-dimensional 3DVAR analysis is the unique global minimizer of $J$. If $m^f,B$ are the mean and covariance of a Gaussian forecast and the observation model is $Y=HX+r$ with independent $r\sim N_p(0,\mathrm{\Gamma })$, then $m^a$ is the unique maximum-a-posteriori point of the posterior density. It is also the mean update of the exact one-step Kalman filter with identity state transition, zero process covariance, and datum $y$.

Let $d,p\ge 1$ be integers, let $\mathrm{\Psi }:{\mathbb{R}}^d\to {\mathbb{R}}^d$, $H\in {\mathbb{R}}^{p\times d}$, and fixed symmetric positive-definite $B\in {\mathbb{R}}^{d\times d}$, $\mathrm{\Gamma }\in {\mathbb{R}}^{p\times p}$. Define $S=HB{H}^T+\mathrm{\Gamma }$, $K=B{H}^T{S}^{-1}$, and $F=(I-KH)\circ \mathrm{\Psi }$. Suppose $F$ is globally Lipschitz with constant $a$ in a norm $\mathrm{\parallel }\cdot \mathrm{\parallel }$ on ${\mathbb{R}}^d$, where $0\le a<1$. Let the true path satisfy $v_{j+1}^{†}=\mathrm{\Psi }(v_j^{†})$, let $y_{j+1}=H{v}_{j+1}^{†}+{ϵ}_j$, and assume ${\sup }_{j\ge 0}\mathrm{\parallel }{ϵ}_j{\mathrm{\parallel }}_2\le ϵ<\mathrm{\infty }$. Starting from any $m_0\in {\mathbb{R}}^d$, run the 3DVAR update $m_{j+1}=(I-KH)\mathrm{\Psi }(m_j)+K{y}_{j+1}\mathrm{.}$ Let $e_1,\dots ,e_p$ be the standard basis of ${\mathbb{R}}^p$ and set $\kappa :={\sum }_{r=1}^p\mathrm{\parallel }K{e}_r\mathrm{\parallel }<\mathrm{\infty }$. Then ${\mathrm{lim\hspace{0.17em}sup}}_{j\to \mathrm{\infty }}\mathrm{\parallel }{m}_j-v_j^{†}\mathrm{\parallel }\le \frac{\kappa }{1-a}ϵ\mathrm{.}$

Assume AC and a finite-horizon model from Linear-Gaussian state-space models with independent noises with $X_k=F_k{X}_{k-1}+q_k$, $Y_k=H_k{X}_k+r_k$, positive-semidefinite $P_0,Q_k$, and positive-definite observation covariances $R_k$. Fix data values $y_{1:T}$. For each $N\ge 2$, initialize the perturbed-observation EnKF of Perturbed-observation ensemble Kalman filter with iid $N(m_0,P_0)$ members, and at each $k$ use independent per-member draws $q_{N,i}^k\sim N(0,Q_k)$ and ${\eta }_{N,i}^k\sim N(0,R_k)$, independent across times and members and of the initial ensemble. Use the divisor-$N$ sample covariance specified in that algorithm. Write ${\bar{X}}_N^{a,k}$ and $C_N^{a,k}$ for the updated empirical mean and covariance, and ${\pi }_k(\cdot ;y_{1:k})=N(m_{k\mathrm{\mid }k},P_{k\mathrm{\mid }k})$ for the exact Kalman filtering law. Then, for each fixed $k\le T$, ${\bar{X}}_N^{a,k}\stackrel{\mathbb{P}}{\to }{m}_{k\mathrm{\mid }k},C_N^{a,k}\stackrel{\mathbb{P}}{\to }{P}_{k\mathrm{\mid }k}\mathrm{.}$ Moreover, for every fixed bounded Lipschitz function $\phi :{\mathbb{R}}^{d_k}\to \mathbb{R}$, $N^{-1}{\sum }_{i=1}^N\phi (X_{N,i}^{a,k})\stackrel{\mathbb{P}}{\to }\int \phi (x){\pi }_k(dx;y_{1:k})\mathrm{.}$ The horizon is fixed while $N\to \mathrm{\infty }$; no uniform-in-time error bound is asserted. The exact conclusion is restricted to linear-Gaussian dynamics and observations. It does not make nonlinear finite-ensemble EnKF exact.