Uniform error bounds of

the ensemble transform Kalman filter

for chaotic dynamics

with multiplicative covariance inflation

Kota Takeda , Takashi Sakajo

1. Department of Mathematics, Kyoto University

2. Data assimilation research team, RIKEN R-CCS

MSJ Autumn Meeting 2024 at Osaka University, Sep. 6, 2024

1,2

The first author is supported by RIKEN Junior Research Associate Program and JST SPRING JPMJSP2110.

The second author is supported by JST MIRAI JPMJMI22G1.

Setup

\frac{d x^\dagger}{dt} = \mathcal{F}(x^\dagger)

(\mathcal{F}: \mathbb{R}^{N_x} \rightarrow \mathbb{R}^{N_x})

y_j = H x^\dagger_j + \xi_j

(H: \mathbb{R}^{N_y \times N_x}, \xi_j \in \mathbb{R}^{N_y} \text{ noise})

x^\dagger_j = x^\dagger(jh)

Model dynamics

Observation

\mathbb{R}^{N_x}

\mathbb{R}^{N_y}

(h > 0 \text{ obs. interval})

unknown

Aim

x^\dagger_j

Unkown

Given

Y_j = \{y_i \mid i \le j\}

x_j

Construct

x^\dagger_j

estimate

Quality Assessment

\mathbb{E}\left[|x_j - x^\dagger_j|^2\right]

(\text{i.e. } \mathbb{P}^{x_j} \approx \mathbb{P}^{x^\dagger_j}({}\cdot{} | Y_j))

y_j

\mathbb{R}^{N_x}

\mathbb{R}^{N_y}

x^\dagger_{j-1}

x^\dagger_j

y_{j-1}

Given

Model

Obs.

unknown

Sequential Data Assimilation

\mathbb{P}^{x_{j-1}}

\mathbb{P}^{\hat{x}_j}

\mathbb{P}^{x_j}

(I) Prediction

by Bayes' rule

by model

(II) Analysis

y_j

We can construct the exact distribution by 2 steps at each time step.

Prediction

\mathbb{P}^{x_1}

\mathbb{P}^{\hat{x}_1}

\mathbb{P}^{x_0}

...

Data Assimilation Algorithms

Only focus on Approximate Gaussian Algorithms

Data Assimilation Algorithms

\mathbb{P}^{x_j}　\approx \frac{1}{m} \sum_{k=1}^m \delta_{x_j^{(k)}}.

Kalman filter(KF) [Kalman 1960]

All Gaussian distribution
Update of mean and covariance.

Assume: linear F, Gaussian noises

Ensemble Kalman filter(EnKF)

Non-linear extension by Monte Carlo

Update of ensemble

\mathbb{P}^{x_j} = \mathcal{N}(m_j, C_j).

X_j = [x_j^{(1)},\dots,x_j^{(m)}].

Data Assimilation Algorithms

Kalman filter(KF)

Assume: linear F, Gaussian noises

Ensemble Kalman filter (EnKF)

Non-linear extension by Monte Carlo method

m_j, C_j

X_j

\widehat{m}_j, \widehat{C}_j

m_{j-1}, C_{j-1}

y_j

X_{j-1}

\widehat{X}_j

y_j

X_j = [x_j^{(k)}]_{k=1}^m, \widehat{X}_j = [\widehat{x}_j^{(k)}]_{k=1}^m.

Linear F

\text{nonlinear } \mathcal{F}

Bayes' rule

\widehat{m}_j, \widehat{C}_j

m_j, C_j

approx.

ensemble

\approx

predict

analysis

Two major implementations:
- Perturbed Observation (PO): stochastic
- Ensemble Transform Kalman Filter (ETKF): deterministic

EnKF (illustration)

PO (stochastic EnKF)

\widehat{x}^{(k)}_j = \Psi(x^{(k)}_{j-1}), \quad k = 1, \dots, m.

y^{(k)}_j = y_j + \xi_j^{(k)}, \quad \xi_j^{(k)} \sim N(0, R), \quad k = 1, \dots, m.

(II) \widehat{X}_j, y_j \rightarrow X_j: \text{Replicate observations,}\\

(I) X_{j-1} \rightarrow \widehat{X}_j: \text{ Evolve each member by model dynamics } \Psi,

x_j^{(k)} = \widehat{x}_j^{(k)} + K_j(y_j^{(k)} - H \widehat{x}_j^{(k)}), \quad k = 1, \dots, m,

\text{Then, update each member,}

Stochastic!

(Burgers+1998) G. Burgers, P. J. van Leeuwen, and G. Evensen, Analysis Scheme in the Ensemble Kalman Filter,684, Mon. Weather Rev., 126 (1998), pp. 1719–1724.

\text{where } K_j = \widehat{C}_j H^\top (H \widehat{C}_j H^\top + R)^{-1}, \widehat{C}_j = \operatorname{Cov}(\widehat{X}_j),

ETKF (deterministic EnKF)

\widehat{x}^{(k)}_j = \Psi(x^{(k)}_{j-1}), \quad k = 1, \dots, m.

\text{with the transform matrix } T_j \in \mathbb{R}^{m \times m} \text{ s.t. }

(II) \widehat{X}_j, y_j \rightarrow X_j: \text{Decompose } \widehat{X}_j =\overline{\widehat{x}}_j + d\widehat{X}_j.\\

(I) X_{j-1} \rightarrow \widehat{X}_j: \text{ Evolve each member by model dynamics } \Psi,

\overline{x}_j = \overline{\widehat{x}}_j^{(k)} + K_j(y_j^{(k)} - H \overline{\widehat{x}}_j^{(k)}), \quad dX_j = d\widehat{X}_j T_j.

\text{Then, update separately,}

Deterministic!

(Bishop+2001)C. H. Bishop, B. J. Etherton, and S. J. Majumdar, Adaptive Sampling with the Ensemble Transform681
Kalman Filter. Part I: Theoretical Aspects, Mon. Weather Rev., 129 (2001), pp. 420–436.

\frac{1}{m-1} d\widehat{X}_j T_j (d\widehat{X}_j T_j)^* = (I - K_j H) \widehat{C}_j.

Inflation (numerical technique)

Issue: Underestimation of covariance through the DA cycle.

→ Poor state estimation. → (idea) inflating covariance before step (II).

Two basic methods of inflation

\text{Additive } (\alpha \ge 0): \widehat{C}_j \rightarrow \widehat{C}_j + \alpha^2 I_{N_x} \text{ (only for PO)}.

\text{Multiplicative } (\alpha \ge 1): \widehat{C}_j \rightarrow \alpha^2\widehat{C}_j \text{ (PO)},

d\widehat{X}_j \rightarrow \alpha d\widehat{X}_j \text{ (ETKF)}.

\text{where } d\widehat{X}_j = [\widehat{x}^{(k)}_j - \overline{\widehat{x}}_j]_{k=1}^m, \overline{\widehat{x}}_j = \frac{1}{m} \sum_{k=1}^m \widehat{x}^{(k)}_j.

(I) Prediction → inflation → (II) Analysis → ...

Literature review

Consistency (Mandel+2011, Kwiatkowski+2015)
Sampling errors (Al-Ghattas+2024)
Error analysis (full observation )
- PO + add. inflation (Kelly+2014)
- ETKF + multi. inflation (T.+2024) → current
Error analysis (partial observation) → future

(Kelly+2014) D. T. B. Kelly, K. J. H. Law, and A. M. Stuart, Well-posedness and accuracy of the ensemble716
Kalman filter in discrete and continuous time, Nonlinearity, 27 (2014), pp. 2579–260.

(Al-Ghattas+2024) O. Al-Ghattas and D. Sanz-Alonso, Non-asymptotic analysis of ensemble Kalman updates: Effective669
dimension and localization, Information and Inference: A Journal of the IMA, 13 (2024).

(Kwiatkowski+2015) E. Kwiatkowski and J. Mandel, Convergence of the square root ensemble kalman filter in the large en-719
semble limit, Siam-Asa J. Uncertain. Quantif., 3 (2015), pp. 1–17.

(Mandel+2011) J. Mandel, L. Cobb, and J. D. Beezley, On the convergence of the ensemble Kalman filter, Appl.739 Math., 56 (2011), pp. 533–541.

H = I_{N_x}

How EnKF approximate KF?

How EnKF approximate

the true state?

Mathematical analysis of EnKF

Assumptions on model

(Reasonable)

Assumption on observation

(Strong)

Assumptions for analysis

N_y = N_x, H = I_{N_x}, \xi \sim \mathcal{N}(0, r^2 I_{N_x})

Assumption (obs-1)

full observation

Strong!

\exist \rho > 0 \text{ s.t. } \forall x_0 \in B(\rho), \forall t > 0, x(t ; x_0) \in B(\rho).

\exist \beta > 0 \text{ s.t. } \forall x \in B(\rho), \forall x' \in \mathbb{R}^{N_x}, \langle F(x) - F(x'), x - x' \rangle \le \beta | x - x'|^2.

\exist \beta > 0 \text{ s.t. } \forall x, x' \in B(\rho), |F(x) - F(x')| \le \beta | x - x'|.

Assumption (model-1)

Assumption (model-2)

Assumption (model-2')

\forall x_0 \in \R^{N_x}, \dot{x}(t) = F(x(t)) \text{ has an unique solution, and}

Reasonable

(B(\rho) = \{x \in \R^{N_x} \mid |x| \le \rho \})

Dissipative dynamical system

Dissipative dynamical systems

Lorenz' 63, 96 equations, widely used in geoscience,

\frac{dx^i}{dt} = (x^{i+1} - x^{i-2}) x^{i-1} - x^i + f \quad (i = 1, \dots, 40),

Example (Lorenz 96)

\bm{x} = (x^i)_{i=1}^{40} \in \mathbb{R}^{40}, f \in \R,

non-linear conserving

linear dissipating

forcing

Assumption (model-1, 2, 2')

hold

https://kotatakeda.github.io/lorenz-webgl/

Example (Lorenz 63)

(incompressible 2D Navier-Stokes equations in infinite-dim. setting)

Error analysis of PO (prev.)

Theorem (Kelly+2014)

\text{Let us consider the PO method with additive inflation } \alpha > 0.

\limsup_{j \rightarrow \infty} e^{(k)}_j \le \frac{2\min\{m, N_x\}}{1-\theta}r^2

\text{Then, } \exist \theta = \theta(\beta, h, r, \alpha) \in (0, 1), \text{ s.t. }

\text{Let } e^{(k)}_j = \mathbb{E}\left[|x^{(k)}_j - x^\dagger_j|^2\right] \text{ for } k = 1, \dots, m.

\text{Let Assumptions (obs-1, model-1,2) hold and }\alpha \gg 1.

variance of observation noise

uniformly in time

\text{for } k = 1, \dots, m.

Error analysis of ETKF (our)

\limsup_{j \rightarrow \infty} e_j \le \frac{N_x r^2}{1-\theta} + \left(\frac{1 - \Theta}{1 - \theta} - 1\right)D.

Theorem (T.+2024)

\text{Let } e_j = \mathbb{E}\left[|\overline{x}_j - x^\dagger_j|^2\right] \text{ where } \overline{x}_j = \frac{1}{m} \sum_{k=1}^m x^{(k)}_j.

\text{where } \Theta = \Theta(\beta, \rho, h, r, m, \alpha) , D = D(\beta, \rho).

\text{Let us consider the ETKF with multiplicative inflation } \alpha \ge 1.

\text{Then } \exist \theta = \theta(\beta, \rho, h, r, m, \alpha) < 1\text{ s.t.}

\text{Let Assumptions (obs-1, model-1,2') hold.}

\text{Suppose } m > N_x \text{ so that } \operatorname{rank}([x^{(k)}_0]_{k=1}^m) = N_x \text{ and } \alpha \gg 1,

Sketch: Error analysis of ETKF (our)

(1) \text{ If } \widehat{C}_j \text{ is positive-difinite, we have a bound by } \alpha \gg 1.

(2) \text{ We also have a lower bound of eigenvalues of } \widehat{C}_j \text{ by } \alpha \gg 1.

Summary

\limsup_{j \rightarrow \infty} e^{(k)}_j \le \frac{2\min\{m, N_x\}}{1-\theta}r^2,

\limsup_{j \rightarrow \infty} e_j \le \frac{N_x r^2}{1-\theta} + \left(\frac{1 - \Theta}{1 - \theta} - 1\right)D.

\text{Suppose } m > N_x \text{ and } \alpha \gg 1,

\text{Suppose } \alpha \gg 1,

ETKF with multiplicative inflation (our)

PO with additive inflation (prev.)

= O(r^2)

r \rightarrow 0

accurate observation limit

= O(r^2)

O(r^4)

\text{for } k = 1, \dots, m.

Two approximate Gaussian filters (EnKF) are assessed in terms of the state estimation error.

due to

Future

\|(I + \widehat{C}_j H^\top R^{-1} H)^{-1}\|_{op} \ge 1.

\text{If } m <= N_x \text{ or } N_y < N_x, H \in \mathbb{R}^{N_y}, R \in \mathbb{R}^{N_y \times N_y},

Issue

\text{then,}

Error contraction in (II)

Issue: partial observation

y = (x^1, x^2) + \xi

x = (x^1, x^2, x^3)

Partially observed

H = I_{N_x}

In general, we can obtain partial information about the state.

Example

(full observation) in Assumption (obs-1) is unrealistic

(H \text{ is a projection.})

true state

observation