Ensemble Kalman filter

Let us review first the Kalman filter. Let ${\displaystyle \mathbf {x} }$ denote the ${\displaystyle n}$-dimensional state vector of a model, and assume that it has Gaussian probability distribution with mean ${\displaystyle \mathbf {\mu } }$ and covariance ${\displaystyle Q}$, i.e., its pdf is

${\displaystyle p(\mathbf {x} )\propto \exp \left(-{\frac {1}{2}}(\mathbf {x} -\mathbf {\mu } )^{\mathrm {T} }Q^{-1}(\mathbf {x} -\mathbf {\mu } )\right).}$

Here and below, ${\displaystyle \propto }$ means proportional; a pdf is always scaled so that its integral over the whole space is one. This ${\displaystyle p(\mathbf {x} )}$, called the prior, was evolved in time by running the model and now is to be updated to account for new data. It is natural to assume that the error distribution of the data is known; data have to come with an error estimate, otherwise they are meaningless. Here, the data ${\displaystyle \mathbf {d} }$ is assumed to have Gaussian pdf with covariance ${\displaystyle R}$ and mean ${\displaystyle H\mathbf {x} }$, where ${\displaystyle H}$ is the so-called observation matrix. The covariance matrix ${\displaystyle R}$ describes the estimate of the error of the data; if the random errors in the entries of the data vector ${\displaystyle \mathbf {d} }$ are independent, ${\displaystyle R}$ is diagonal and its diagonal entries are the squares of the standard deviation (“error size”) of the error of the corresponding entries of the data vector ${\displaystyle \mathbf {d} }$. The value ${\displaystyle H\mathbf {x} }$ is what the value of the data would be for the state ${\displaystyle \mathbf {x} }$ in the absence of data errors. Then the probability density ${\displaystyle p(\mathbf {d} |\mathbf {x} )}$ of the data ${\displaystyle \mathbf {d} }$ conditional of the system state ${\displaystyle \mathbf {x} }$, called the data likelihood, is

${\displaystyle p\left(\mathbf {d} |\mathbf {x} \right)\propto \exp \left(-{\frac {1}{2}}(\mathbf {d} -H\mathbf {x} )^{\mathrm {T} }R^{-1}(\mathbf {d} -H\mathbf {x} )\right).}$

The pdf of the state and the data likelihood are combined to give the new probability density of the system state ${\displaystyle \mathbf {x} }$ conditional on the value of the data ${\displaystyle \mathbf {d} }$ (the posterior) by the Bayes theorem,

${\displaystyle p\left(\mathbf {x} |\mathbf {d} \right)\propto p\left(\mathbf {d} |\mathbf {x} \right)p(\mathbf {x} ).}$

The data ${\displaystyle \mathbf {d} }$ is fixed once it is received, so denote the posterior state by ${\displaystyle \mathbf {\hat {x}} }$ instead of ${\displaystyle \mathbf {x} |\mathbf {d} }$ and the posterior pdf by ${\displaystyle p\left(\mathbf {\hat {x}} \right)}$. It can be shown by algebraic manipulations[5] that the posterior pdf is also Gaussian,

${\displaystyle p\left(\mathbf {\hat {x}} \right)\propto \exp \left(-{\frac {1}{2}}(\mathbf {\hat {x}} -\mathbf {\hat {\mu }} )^{\mathrm {T} }{\hat {Q}}^{-1}(\mathbf {\hat {x}} -\mathbf {\hat {\mu }} )\right),}$

with the posterior mean ${\displaystyle \mathbf {\hat {\mu }} }$ and covariance ${\displaystyle {\hat {Q}}}$ given by the Kalman update formulas

${\displaystyle \mathbf {\hat {\mu }} =\mathbf {\mu } +K\left(\mathbf {d} -H\mathbf {\mu } \right),\quad {\hat {Q}}=\left(I-KH\right)Q,}$

where

${\displaystyle K=QH^{\mathrm {T} }\left(HQH^{\mathrm {T} }+R\right)^{-1}}$

is the so-called Kalman gain matrix.

The EnKF is a Monte Carlo approximation of the Kalman filter, which avoids evolving the covariance matrix of the pdf of the state vector ${\displaystyle \mathbf {x} }$. Instead, the pdf is represented by an ensemble

${\displaystyle X=\left[\mathbf {x} _{1},\ldots ,\mathbf {x} _{N}\right]=\left[\mathbf {x} _{i}\right].}$

${\displaystyle X}$ is an ${\displaystyle n\times N}$ matrix whose columns are the ensemble members, and it is called the prior ensemble. Ideally, ensemble members would form a sample from the prior distribution. However, the ensemble members are not in general independent except in the initial ensemble, since every EnKF step ties them together. They are deemed to be approximately independent, and all calculations proceed as if they actually were independent.

Replicate the data ${\displaystyle \mathbf {d} }$ into an ${\displaystyle m\times N}$ matrix

${\displaystyle D=\left[\mathbf {d} _{1},\ldots ,\mathbf {d} _{N}\right]=\left[\mathbf {d} _{i}\right],\quad \mathbf {d} _{i}=\mathbf {d} +\mathbf {\epsilon _{i}} ,\quad \mathbf {\epsilon _{i}} ~N(0,R),}$

so that each column ${\displaystyle \mathbf {d} _{i}}$ consists of the data vector ${\displaystyle \mathbf {d} }$ plus a random vector from the ${\displaystyle m}$-dimensional normal distribution ${\displaystyle N(0,R)}$. If, in addition, the columns of ${\displaystyle X}$ are a sample from the prior probability distribution, then the columns of

${\displaystyle {\hat {X}}=X+K(D-HX)}$

form a sample from the posterior probability distribution. To see this in the scalar case with ${\displaystyle H=1}$: Let ${\displaystyle x_{i}=\mu +\xi _{i},\;\xi _{i}\sim N(0,\sigma _{x}^{2})}$, and ${\displaystyle d_{i}=d+\epsilon _{i},\;\epsilon _{i}\sim N(0,\sigma _{d}^{2}).}$ Then

${\displaystyle {\hat {x}}_{i}=\left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\mu +{\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}d\right)+\left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\xi _{i}+{\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\epsilon _{i}\right)}$.

The first sum is the posterior mean, and the second sum, in view of the independence, has a variance

${\displaystyle \left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\right)^{2}\sigma _{x}^{2}+\left({\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\right)^{2}\sigma _{d}^{2}={\frac {1}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}}$,

which is the posterior variance.

The EnKF is now obtained simply by replacing the state covariance ${\displaystyle Q}$ in Kalman gain matrix ${\displaystyle K}$ by the sample covariance ${\displaystyle C}$ computed from the ensemble members (called the ensemble covariance),[6] that is: ${\displaystyle K=CH^{\mathrm {T} }\left(HCH^{\mathrm {T} }+R\right)^{-1}}$

Here we follow.[7][8] Suppose the ensemble matrix ${\displaystyle X}$ and the data matrix ${\displaystyle D}$ are as above. The ensemble mean and the covariance are

${\displaystyle E\left(X\right)={\frac {1}{N}}\sum _{k=1}^{N}\mathbf {x} _{k},\quad C={\frac {AA^{T}}{N-1}},}$

where

${\displaystyle A=X-E\left(X\right)\mathbf {e} _{1\times N}=X-{\frac {1}{N}}\left(X\mathbf {e} _{N\times 1}\right)\mathbf {e} _{1\times N},}$

and ${\displaystyle \mathbf {e} }$ denotes the matrix of all ones of the indicated size.

The posterior ensemble ${\displaystyle X^{p}}$ is then given by

${\displaystyle X^{p}=X+CH^{T}\left(HCH^{T}+R\right)^{-1}(D-HX),}$

where the perturbed data matrix ${\displaystyle D}$ is as above.

Note that since ${\displaystyle R}$ is a covariance matrix, it is always positive semidefinite and usually positive definite, so the inverse above exists and the formula can be implemented by the Cholesky decomposition.[9] In,[7][8] ${\displaystyle R}$ is replaced by the sample covariance ${\displaystyle {\tilde {D}}{\tilde {D}}^{T}/\left(N-1\right)}$ where ${\displaystyle {\tilde {D}}=D-{\frac {1}{N}}d\,\mathbf {e} _{1\times N}}$and the inverse is replaced by a pseudoinverse, computed using the singular-value decomposition (SVD) .

Since these formulas are matrix operations with dominant Level 3 operations,[10] they are suitable for efficient implementation using software packages such as LAPACK (on serial and shared memory computers) and ScaLAPACK (on distributed memory computers).[9] Instead of computing the inverse of a matrix and multiplying by it, it is much better (several times cheaper and also more accurate) to compute the Cholesky decomposition of the matrix and treat the multiplication by the inverse as solution of a linear system with many simultaneous right-hand sides.[10]

Since we have replaced the covariance matrix with ensemble covariance, this leads to a simpler formula where ensemble observations are directly used without explicitly specifying the matrix ${\displaystyle H}$. More specifically, define a function ${\displaystyle h(\mathbf {x} )}$ of the form

${\displaystyle h(\mathbf {x} )=H\mathbf {x} .}$

The function ${\displaystyle h}$ is called the observation function or, in the inverse problems context, the forward operator. The value of ${\displaystyle h(\mathbf {x} )}$ is what the value of the data would be for the state ${\displaystyle \mathbf {x} }$ assuming the measurement is exact. Then the posterior ensemble can be rewritten as

${\displaystyle X^{p}=X+{\frac {1}{N-1}}A\left(HA\right)^{T}P^{-1}(D-HX)}$

where

${\displaystyle HA=HX-{\frac {1}{N}}\left(\left(HX\right)\mathbf {e} _{N\times 1}\right)\mathbf {e} _{1\times N},}$

and

${\displaystyle P={\frac {1}{N-1}}HA\left(HA\right)^{T}+R,}$

with

${\displaystyle \left[HA\right]_{i}=H\mathbf {x} _{i}-H{\frac {1}{N}}\sum _{j=1}^{N}\mathbf {x} _{j}\ =h\left(\mathbf {x} _{i}\right)-{\frac {1}{N}}\sum _{j=1}^{N}h\left(\mathbf {x} _{j}\right).}$

Consequently, the ensemble update can be computed by evaluating the observation function ${\displaystyle h}$ on each ensemble member once and the matrix ${\displaystyle H}$ does not need to be known explicitly. This formula holds also[9] for an observation function ${\displaystyle h(\mathbf {x} )=H\mathbf {x+f} }$ with a fixed offset ${\displaystyle \mathbf {f} }$, which also does not need to be known explicitly. The above formula has been commonly used for a nonlinear observation function ${\displaystyle h}$, such as the position of a hurricane vortex.[11] In that case, the observation function is essentially approximated by a linear function from its values at ensemble members.

For a large number ${\displaystyle m}$ of data points, the multiplication by ${\displaystyle P^{-1}}$ becomes a bottleneck. The following alternative formula is advantageous when the number of data points ${\displaystyle m}$ is large (such as when assimilating gridded or pixel data) and the data error covariance matrix ${\displaystyle R}$ is diagonal (which is the case when the data errors are uncorrelated), or cheap to decompose (such as banded due to limited covariance distance). Using the Sherman–Morrison–Woodbury formula[12]

${\displaystyle (R+UV^{T})^{-1}=R^{-1}-R^{-1}U(I+V^{T}R^{-1}U)^{-1}V^{T}R^{-1},}$

with

${\displaystyle U={\frac {1}{N-1}}HA,\quad V=HA,}$

gives

${\displaystyle {\begin{aligned}P^{-1}&=\left(R+{\frac {1}{N-1}}HA\left(HA\right)^{T}\right)^{-1}\ =\\&=R^{-1}\left[I-{\frac {1}{N-1}}\left(HA\right)\left(I+\left(HA\right)^{T}R^{-1}{\frac {1}{N-1}}\left(HA\right)\right)^{-1}\left(HA\right)^{T}R^{-1}\right],\end{aligned}}}$

which requires only the solution of systems with the matrix ${\displaystyle R}$ (assumed to be cheap) and of a system of size ${\displaystyle N}$ with ${\displaystyle m}$ right-hand sides. See[9] for operation counts.