Uniform error bounds of

the ensemble square root filter for

chaotic dynamics with

multiplicative covariance inflation

The 1st MMS Workshop for Young Researchers 2024/11/20-21 @Kyoto University

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

The second author is supported by JST MIRAI JPMJMI22G1.

Kota Takeda  ,   Takashi Sakajo

\dagger
\ddagger
\dagger

 Kyoto University,    RIKEN R-CCS

\dagger
\ddagger
  • What's Data Assimilation (DA)?
  • Setting: model, observation
  • Problem: Sequential state estimation (filtering)
  • Approach & Algorithms: EnKF
  • Results (Error analysis): prev. & our
  • Discussion/Further analysis

Contents

Data Assimilation (DA)

Seamless integration of data and models

Model

Dynamical systems and its numerical simulations

ODE, PDE, SDE,

Numerical analysis, HPC,

Fluid mechanics, Physics

Data

Obtained from experiments or measurements

Statistics, Data analysis,

Radiosonde/Radar/Satellite/ observation

×

\partial_t \bm{u} + (\bm{u} \cdot \nabla) \bm{u} = - \frac{1}{\rho} \nabla p + \nu \triangle \bm{u} + \bm{f}

Numerical Weather Prediction

typical application

×

Data Assimilation (DA)

  1. Control
    • ​Nudging/Continuous Data Assimilation
  2. Bayesian
    • Particle Filter (PF)
    • ​Ensemble Kalman filter (EnKF)
  3. Variational
    • ​​Four-dimensional Variation method (4DVar)

Several formulations & algorithms of DA:

focus on

Setup: Model

\frac{d u}{dt} = \mathcal{F}(u)
(\mathcal{F}: \mathcal{U} \rightarrow \mathcal{U}, u(0) = u_0 \in \mathcal{U})
u(t) = \Psi_t (u_0)

Model dynamics in

Semi-group in

\mathcal{U}
\mathcal{U}

generates

\mathcal{U}
(\Psi_t : \mathcal{U} \rightarrow \mathcal{U}, t \ge 0)

(assume)

: separable Hilbert space

: separable Hilbert space (state space)

The known model generates

the unknown true chaotic orbit.

t

A small perturbation can

grow exponentially fast.

unknown orbit

unknown

Chaotic dynamical systems

The orbits of the Lorenz 63 equation.

Setup: Observation

Discrete-time model (at obs. time)

\mathcal{U}
(\Psi = \Psi_h, h > 0 \text{: observation interval})

We have noisy observations in

t
y_j = H u_j + \xi_j, \quad j \in \mathbb{N}.
(H: \mathcal{L}(\mathcal{U}, \mathcal{Y}), \xi_j \sim N(0, R))

Noisy observation

u_j = \Psi(u_{j-1})

observe

u_1
y_1
u_2
y_2
t_1
t_2
t_3

...

h

bounded linear

\in \mathcal{Y}
\mathcal{Y}

(another Hilbert space).

\in \mathcal{Y}
\Psi

Sequential state estimation

\mathcal{U}
t
y_1
y_2
t_1
t_2
t_3

...

h

given

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

Problem

u_j

estimate

\text{For } j \in \mathbb{N},

known:

\Psi,
H,

obs. noise distribution.

y_2
y_3
u_1
u_2
u_3
\Psi
\Psi
  • What's DA? ✓
  • Setting: model, observation ✓
  • Problem: Sequential state estimation (filtering) ✓
  • Approach & Algorithms:
    • Bayesian approach → Gaussian approximation (KF) → Monte Carlo (EnKF)
    • EnKF: PO (Stoc.), ETKF (Det.) + inflation (technique)
  • Results (Error analysis)
  • Discussion/Further analysis

Contents

Bayesian Data Assimilation

Approach: Bayesian Data Assimilation

\mathbb{P}^{v_j} = \mathbb{P}^{u_j}({}\cdot{} | Y_j)

construct a random variable      so that

v_j

given

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

Problem

u_j

estimate

Sequential Data Assimilation

(I) Prediction

\Psi

(II)Analysis

We can construct it iteratively (2 steps at each time step)

(I) Prediction

(II) Analysis

Prediction

\mathbb{P}^{v_{j-1}}
\mathbb{P}^{\hat{v}_j}
\mathbb{P}^{v_j}

...

by model

by Bayes' rule

y_j
\mathcal{U}
y_1
y_2
t
t_0
t_1
t_2
t_3
h

with auxiliary variable

\hat{v}_j

Repeat (I) & (II)

\Rightarrow \mathbb{P}^{v_j} = \mathbb{P}^{u_j}({}\cdot{} | Y_j)

Proposition

\text{Assume } \mathbb{P}^{v_0} = \mathbb{P}^{u_0}

Need approximations!

→ next slide

Data Assimilation Algorithms

Only focus on Approximate Gaussian Algorithms

KF & EnKF

Kalman filter (KF) [Kalman 1960]

  • All Gaussian distribution
  • Update of mean and covariance.

Assume: linear F, Gaussian noises

\mathbb{P}^{v_j} = \mathcal{N}(m_j, C_j).
\mathbb{P}^{v_j} \approx \frac{1}{m} \sum_{k=1}^m \delta_{v_j^{(k)}}.

Ensemble Kalman filter (EnKF)

Non-linear extension by Monte Carlo method

  •  
  • Update of ensemble
V_j = [v_j^{(1)},\dots,v_j^{(m)}] \in \mathbb{R}^{N_x \times m}

(Reich+2015) S. Reich and C. Cotter (2015), Probabilistic Forecasting and Bayesian Data Assimilation, Cambridge University Press, Cambridge.

(Law+2015) K. J. H. Law, A. M. Stuart, and K. C. Zygalakis (2015), Data Assimilation: A Mathematical Introduction, Springer.

KF & EnKF

Kalman filter(KF)

Assume: linear F, Gaussian noises

m_j, C_j
\widehat{m}_j, \widehat{C}_j
m_{j-1}, C_{j-1}
y_j

Linear F

Bayes' rule

predict

analysis

Ensemble Kalman filter (EnKF) 

Non-linear extension by Monte Carlo method

V_j
V_{j-1}
\widehat{V}_j
y_j
V_j = [v_j^{(k)}]_{k=1}^m, \widehat{V}_j = [\widehat{v}_j^{(k)}]_{k=1}^m.
\text{nonlinear } \mathcal{F}
\widehat{m}_j, \widehat{C}_j
m_j, C_j

approx.

ensemble

\approx
\approx
m_j = \overline{v}_j = \frac{1}{m} \sum_{k=1}^m v^{(k)}_j
C_j = \operatorname{Cov}_m(V_j) = \frac{1}{m-1} \sum_{k=1}^m (v^{(k)}_j - \overline{v}_j) \otimes (v^{(k)}_j - \overline{v}_j)

Two variants

→ next slide

  • Two major variants of EnKF:
    • Perturbed Observation (PO): stochastic → sampling error
    • Ensemble Transform Kalman Filter (ETKF): deterministic

Two EnKF variants

\mathcal{U}
t
y_1
y_2
t_1
t_2
t_3

...

h
t_0

Repeat (I) & (II)

(II)Analysis

(I) Prediction

\Psi
v_0^{(1)}
v_0^{(m)}
v_0^{(2)}
\widehat{v}_1^{(1)}
v_1^{(1)}
\vdots

(Burgers+1998)

(Bishop+2001)

Two variants

PO (stochastic EnKF)

\widehat{v}^{(k)}_j = \Psi(v^{(k)}_{j-1}), \quad k = 1, \dots, m.
(I) V_{j-1} \rightarrow \widehat{V}_j: \text{ Evolve each member by model dynamics } \Psi,

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

Stochastic!

y^{(k)}_j = y_j + \xi_j^{(k)}, \quad \xi_j^{(k)} \sim N(0, R), \quad k = 1, \dots, m.
(II) \widehat{V}_j, y_j \rightarrow V_j: \text{Replicate observations,}\\
v_j^{(k)} = (I - K_j H) \widehat{v}_j^{(k)} + K_j y_j^{(k)}, \quad k = 1, \dots, m,
\text{Then, update each member,}
\text{where } K_j = \widehat{C}_j H^\top (H \widehat{C}_j H^\top + R)^{-1}, \widehat{C}_j = \operatorname{Cov}_m(\widehat{V}_j),
(\bigstar)
y^{(k)}_j = y_j + \xi_j^{(k)}, \quad \xi_j^{(k)} \sim N(0, R), \quad k = 1, \dots, m.
(II) \widehat{V}_j, y_j \rightarrow V_j: \text{Replicate observations,}\\
v_j^{(k)} = (I - K_j H) \widehat{v}_j^{(k)} + K_j y_j^{(k)}, \quad k = 1, \dots, m,
\text{Then, update each member,}
\text{where } K_j = \widehat{C}_j H^\top (H \widehat{C}_j H^\top + R)^{-1}, \widehat{C}_j = \operatorname{Cov}_m(\widehat{V}_j),
(\bigstar)

Perturbed Observation

\cdot^\top

: A transpose (adjoint) of a matrix (linear operator)

ETKF (deterministic EnKF)

\widehat{v}^{(k)}_j = \Psi(v^{(k)}_{j-1}), \quad k = 1, \dots, m.
(I) V_{j-1} \rightarrow \widehat{V}_j: \text{ Evolve each member by model dynamics } \Psi,

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

Deterministic!

\text{with the transform matrix } T_j \in \mathbb{R}^{m \times m} \text{ s.t. }
(II) \widehat{V}_j, y_j \rightarrow V_j: \text{Decompose } \widehat{V}_j =\overline{\widehat{v}}_j + d\widehat{V}_j.\\
\overline{v}_j = (I - K_j H)\overline{\widehat{v}}_j + K_j y_j , \quad dV_j = d\widehat{V}_j T_j.
\text{Then, update separately,}
\frac{1}{m-1} d\widehat{V}_j T_j (d\widehat{V}_j T_j)^\top = (I - K_j H) \widehat{C}_j.

The explicit representation exists!

(\bigstar)

Ensemble Transform Kalman Filter

ETKF (deterministic EnKF)

Ensemble Transform Kalman Filter

\text{The transform matrix } T_j \in \mathbb{R}^{m \times m} \text{ satisfying}
\frac{1}{m-1} d\widehat{V}_j T_j (d\widehat{V}_j T_j)^\top = (I - K_j H) \widehat{C}_j

Lemma

\text{is uniquely given by}
T_j = \left(I_m + \frac{1}{m-1} d\widehat{V}_j^\top H^\top R^{-1} H d\widehat{V}_j \right)^{-\frac{1}{2}}.
\text{(up to unitary transformation)}

ETKF is one of a class of deterministic EnKF known as ESRF (Ensemble Square Root Filters).

※ESRF (deterministic EnKF)

ESRF

ETKF

\frac{1}{m-1} d\widehat{V}_j T_j (d\widehat{V}_j T_j)^\top = (I - K_j H) \widehat{C}_j.
T_j \in \mathbb{R}^{m \times m}

EAKF (Ensemble Adjustment Kalman Filter)

\frac{1}{m-1} A_j d\widehat{V}_j (A_j d\widehat{V}_j)^\top = (I - K_j H) \widehat{C}_j.
A_j \in \mathcal{L}(\mathcal{U})

: bounded linear

(Anderson 2001) J. L. Anderson (2001), An Ensemble Adjustment Kalman Filter for Data Assimilation, Mon. Weather Rev., 129, 2884–2903.

Key perspective of EnKF (KF)

v_j = (I - K_j H) \widehat{v}_j + K_j y_j.
(\bigstar)
\text{where } K_j = \widehat{C}_j H^\top (H \widehat{C}_j H^\top + R)^{-1}.

The analysis state     is the weighted average of

the prediction     and observation

v_j
\widehat{v}_j
y_j

The weighted       is determined by estimated uncertainties

K_j
R

: uncertainty in observation (prescribed)

\widehat{C}_j

: uncertainty in prediction (approx. by ensemble)

Inflation (numerical technique)

Issue: Underestimation of covariance       through the DA cycle

→ Overconfident in the prediction → Poor state estimation

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

Rank is improved

Rank is not improved

\widehat{C}_j

Idea: inflating the covariance before (II)

Two basic methods of inflation

\text{Additive } (\alpha \ge 0)
\text{Multiplicative } (\alpha \ge 1)
\text{ETKF}: d\widehat{V}_j \rightarrow \alpha d\widehat{V}_j.
\text{where } d\widehat{V}_j = [\widehat{v}^{(k)}_j - \overline{\widehat{v}}_j]_{k=1}^m, \overline{\widehat{v}}_j = \frac{1}{m} \sum_{k=1}^m \widehat{v}^{(k)}_j.
\text{PO}: \widehat{C}_j \rightarrow \alpha^2\widehat{C}_j,
\text{Only for PO}: \widehat{C}_j \rightarrow \widehat{C}_j + \alpha^2 I_{\mathcal{U}}.

※depending on the implementation

※

PO

Two EnKF variants

Simple

Stochastic

Require many samples

ETKF

Complicated

Deterministic

Require few samples

Practical !

Inflation: add., multi.

Inflation: only multi.

※depending on the implementation

※

  • What's DA? ✓
  • Setting: model, observation ✓
  • Problem: Sequential state estimation (filtering) ✓
  • Approach & Algorithms: ✓
    • Bayesian approach → Gaussian approximation (KF) → Monte Carlo (EnKF)
    • EnKF: PO (Stoc.), ETKF (Det.) + inflation (technique)
  • Results (Error analysis)
  • Discussion/Further analysis

Contents

Error analysis of EnKF

Assumptions on model

(Reasonable)

Assumption on observation

(Strong)

We focus on finite-dimensional state spaces

(\mathcal{U} = \R^{N_u} \text{ for } N_u \in \mathbb{N})

(T.+2024) K. T. & T. Sakajo (2024), SIAM/ASA Journal on Uncertainty Quantification, 12(4), 1315–1335.

Assumptions for analysis

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

Assumption (obs-1)

"full observation"

Strong!

\exist \rho > 0 \text{ s.t. } \forall u_0 \in B(\rho), \forall t > 0, u(t ; u_0) \in B(\rho).
\exist \beta > 0 \text{ s.t. } \forall u \in B(\rho), \forall v \in \mathbb{R}^{N_u}, \langle F(u) - F(v), u - v \rangle \le \beta | u - v|^2.
\exist \beta > 0 \text{ s.t. } \forall u, v \in B(\rho), |F(u) - F(v)| \le \beta |u - v|.

Assumption (model-1)

Assumption (model-2)

Assumption (model-2')

\forall u_0 \in \R^{N_u}, \frac{du}{dt}(t) = F(u(t)) \text{ has an unique solution, and}

Reasonable

(B(\rho) = \{u \in \R^{N_u} \mid |u| \le \rho \})

"one-sided Lipschitz"

"Lipschitz"

"bounded"

or

Dissipative dynamical system

Dissipative dynamical system

Dissipative dynamical systems

Lorenz 63, 96 equations, widely used in geoscience,

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

Example (Lorenz 96)

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

non-linear conserving

linear dissipating

forcing

Assumptions (model-1, 2, 2')

hold

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

Lorenz 63 eq. can be written in the same form.

→ We can discuss ODE and PDE in a unified framework.

can be generalized

\frac{du}{dt} =
B(u, u)
+
A u
+
f.

bilinear operator

linear operator

dissipative dynamical systems

Mimics chaotic variation of

physical quantities at equal latitudes.

(Lorenz1996) (Lorenz+1998)

u

Dissipative dynamical systems

The orbits of the Lorenz 63 equation.

Literature review

  • Consistency (Mandel+2011, Kwiatkowski+2015)

  • Sampling errors (Al-Ghattas+2024)

  • Stability (Tong+2016a,b)

  • Error analysis (full observation                )

    • PO + add. inflation (Kelly+2014)

    • ETKF + multi. inflation (T.+2024) → current

    • EAKF + multi. inflation (T.)

  • Error analysis (partial observation) → future  (partially solved, T.)

H = I_{N_u}

How EnKF approximate KF?

How EnKF approximate

the true state?

※ Continuous time formulation (EnKBF)

  • EnKBF with full-obs. (de Wiljes+2018)

Error analysis of PO (prev.)

\text{Let us consider the PO method with additive inflation } \alpha > 0.
\text{Let } e^{(k)}_j = v^{(k)}_j - u_j \text{ for } k = 1, \dots, m.

(Kelly+2014) D. Kelly et al. (2014), Nonlinearity, 27(10), 2579–2603.

\widehat{C}_j \rightarrow \widehat{C}_j + \alpha^2 I_{N_u}

Strong effect! (improve rank)

Theorem 1 (Kelly+2014)

\mathbb{E} [|e^{(k)}_j|^2] \le \theta^j \mathbb{E} [|e^{(k)}_0|^2] + 2m r^2 \frac{1 - \theta^j}{1-\theta}
\text{Then, } \exist \theta = \theta(\beta, h, r, \alpha) \in (0, 1), \text{ s.t. }
\text{Let Assumptions (obs-1, model-1,2) hold and }\alpha \gg 1.
\text{for } j \in \mathbb{N}, k = 1, \dots, m.

exponential decay

initial error

\mathbb{E} [\cdot] \text{ : expectation w.r.t. obs. noise and initial uncertainty}

Error analysis of PO (prev.)

\text{Let us consider the PO method with additive inflation } \alpha > 0.
\text{Let } e^{(k)}_j = v^{(k)}_j - u_j \text{ for } k = 1, \dots, m.

variance of observation noise

uniformly in time

(Kelly+2014) D. Kelly et al. (2014), Nonlinearity, 27(10), 2579–2603.

\widehat{C}_j \rightarrow \widehat{C}_j + \alpha^2 I_{N_u}

Strong effect! (improve rank)

Corollary 1 (Kelly+2014)

\limsup_{j \rightarrow \infty} \mathbb{E} [|e^{(k)}_j|^2] \le \frac{2 m}{1-\theta}r^2
\text{Then, } \exist \theta = \theta(\beta, h, r, \alpha) \in (0, 1), \text{ s.t. }
\text{Let Assumptions (obs-1, model-1,2) hold and }\alpha \gg 1.
\text{for } k = 1, \dots, m.
\mathbb{E} [\cdot] \text{ : expectation w.r.t. obs. noise and initial uncertainty}

Error analysis of ETKF (our)

\text{Let } e_j =\overline{v}_j - u_j\text{ where } \overline{v}_j = \frac{1}{m} \sum_{k=1}^m v^{(k)}_j.
\text{Let us consider the ETKF with multiplicative inflation } \alpha \ge 1.

(T.+2024) K. T. & T. Sakajo (2024), SIAM/ASA Journal on Uncertainty Quantification, 12(4), 1315–1335.

\mathbb{E}[|e_j|^2] \le \theta^j \mathbb{E}[|e_0|^2] + \frac{N_u r^2(1 - \theta^j)}{1-\theta} + \left(\frac{1 - \Theta}{1 - \theta} - 1\right)D.

Theorem 2 (T.+2024)

\text{where } \Theta = \Theta(\beta, \rho, h, r, m, \alpha) , D = D(\beta, \rho) > 0.
\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_u, \operatorname{rank}(C_0) = N_u,\alpha \gg 1, \text{ and}
v_j^{(k)} \in B(\rho) \text{ for } k = 1, \dots, m \text{ and } j \in \N.

Error analysis of ETKF (our)

\text{Let } e_j =\overline{v}_j - u_j\text{ where } \overline{v}_j = \frac{1}{m} \sum_{k=1}^m v^{(k)}_j.
\text{Let us consider the ETKF with multiplicative inflation } \alpha \ge 1.

(T.+2024) K. T. & T. Sakajo (2024), SIAM/ASA Journal on Uncertainty Quantification, 12(4), 1315–1335.

\limsup_{j \rightarrow \infty} \mathbb{E}[|e_j|^2] \le \frac{N_u r^2}{1-\theta} + \left(\frac{1 - \Theta}{1 - \theta} - 1\right)D.

Corollary 2 (T.+2024)

\text{where } \Theta = \Theta(\beta, \rho, h, r, m, \alpha) , D = D(\beta, \rho) > 0.
\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_u, \operatorname{rank}(C_0) = N_u,\alpha \gg 1, \text{ and}
v_j^{(k)} \in B(\rho) \text{ for } k = 1, \dots, m \text{ and } j \in \N.

Sketch: Error analysis of ETKF (our)

Strategy: estimate the changes in DA cycle

\widehat{e}_j
(\widehat{e}_j = \overline{\widehat{v}}_j - u_j)
(I)
(II)

expanded

by dynamics

contract &

noise added

e_{j-1}
e_j

Applying Gronwall's lemma to the model dynamics,

(I)
| \widehat{e}_j | \le e^{\beta h} |e_{j-1}|

+ {additional term}.

①

\overline{v}_j = (I - K_j H)\overline{\widehat{v}}_j + K_j y_j
e_j = (I - K_j H) \widehat{e}_j + K_j (y_j - H u_j).

subtracting          yields

(II)

From the Kalman update relation

(\bigstar)
u_j

prediction error

② contraction

③ obs. noise effect

estimated by op.-norm

Sketch: Error analysis of ETKF (our)

Key estimates in

Using Woodbury's lemma, we have

I - K_j H = (I + \widehat{C}_j H^\top R^{-1} H)^{-1}.

② contraction

③ obs. noise

(II)
\| \cdot \|_{op}: \text{operator norm}
\|I - K_j H\|_{op} = \|(I + \widehat{C}_j H^\top R^{-1} H)^{-1}\|_{op}

Contraction (see the next slide)

= \|(I + \alpha^2 r^{-2}\widehat{C}_j)^{-1}\|_{op}.
H = I, R = r^2 I

assumption

\|K_j\|_{op} = \|(I + \alpha^2 r^{-2}\widehat{C}_j)^{-1} \alpha^2 r^{-2}\widehat{C}_j\|_{op} \le 1.
(\forall \alpha, r \ge 0, \widehat{C}_j \text{: positive-definite})

Sketch: Error analysis of ETKF (our)

< e^{-\beta h}
= \frac{1}{1 + \alpha^2 r^{-2} \lambda_{min}(\widehat{C}_j)}

①

②

Key result

\because
(1) \text{ We establish a lower bound of eigenvalues of } \widehat{C}_j \text{ by } \alpha \gg 1.
(2) \text{ If } \widehat{C}_j \text{ is positive-difinite, we have a bound by } \alpha \gg 1.

Sketch: Error analysis of ETKF (our)

\text{Let } \widehat{\lambda}_j = \lambda_{min}(\widehat{C}_j), \lambda_j = \lambda_{min}(C_j).
(1) \text{ We establish a lower bound of eigenvalues of } \widehat{C}_j \text{ by } \alpha \gg 1.
\widehat{\lambda}_{j-1} \rightarrow \lambda_j \rightarrow \widehat{\lambda}_j
\Rightarrow \widehat{\lambda}_j \ge \frac{e^{-ch} \alpha^2 \widehat{\lambda}_{j-1}}{1 + \alpha^2 r^{-2} \widehat{\lambda}_{j-1}}
\therefore \liminf_{j \rightarrow \infty} \widehat{\lambda}_j \ge \frac{r^2 (e^{-ch} \alpha^2 - 1)}{\alpha^2} > 0.

estimate the change

\widehat{\lambda}_{j-1} \rightarrow \lambda_j

(inflation & analysis)

\lambda_j = \frac{\alpha^2 \widehat{\lambda}_{j-1}}{1 + \alpha^2 r^{-2} \widehat{\lambda}_{j-1}}
\lambda_j \rightarrow \widehat{\lambda}_j

(prediction)

\widehat{\lambda}_j \ge e^{-ch} \lambda_j
\left(c = \frac{8m \beta\rho^2}{m-1}\right)

Comparison

\limsup_{j \rightarrow \infty} \mathbb{E}[|e^{(k)}_j|^2] \le \frac{2 m}{1-\theta}r^2,
\limsup_{j \rightarrow \infty} \mathbb{E}[|e_j|^2] \le \frac{N_u r^2}{1-\theta} + \left(\frac{1 - \Theta}{1 - \theta} - 1\right)D.
\text{Suppose } m > N_u \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

due to

Numerical result

Discussion

(T.+2024)

1. Partial observation

2. Small ensemble

m \le N_u
N_y < N_u, H \in \mathbb{R}^{N_y \times N_u}

3. Infinite-dimension

\text{state space } \mathcal{U}: \text{Hilbert space},
\text{observation space } \mathcal{Y} \subset \mathcal{U}.

→ Issue 1, 2

→ Issue 1

→ Issue 1

partially solved (T.)

Discussion

\|(I + \widehat{C}_j H^\top R^{-1} H)^{-1}\|_{op} \ge 1.
\text{If } m \le N_u \text{ or } N_y < N_u, H \in \mathbb{R}^{N_y \times N_u}, R \in \mathbb{R}^{N_y \times N_y},

Issue 1 (rank deficient)

\text{then,}

② Error contraction in (II)

1. Partial observation

2. Small ensemble

Discussion

\|(I + \widehat{C}_j H^\top R^{-1} H)^{-1} \widehat{C}_j H^\top R^{-1} H\|_{op} \le 1.
\text{Let } N_y < N_u, H \in \mathbb{R}^{N_y \times N_u}, R \in \mathbb{R}^{N_y \times N_y}.

Issue 2 (non-symmetric)

\text{If } \widehat{C}_j H^\top R^{-1} H \text{ is symmetric then,}

③ bound of obs. noise

If non-symmetric, it doesn't hold in general.

Example

A = \left[ \begin{array}{cc} a & 0 \\ b & 0 \end{array} \right] \Rightarrow \|(I + A)^{-1} A\|_{op} = \frac{\sqrt{a^2 + b^2}}{1+a} > 1.
(a > 0, b^2 \ge 1 + 2a)

Discussion

Infinite-dimension

Theorem (Kelly+2014) for PO + add. infl.

Theorem (T.+2024) for ETKF + multi. infl.

\text{state space } \mathcal{U}: \text{Hilbert space},
\text{observation space } \mathcal{Y} \subset \mathcal{U} \text{: Hilbert space}.

(     improve rank)

\because

still holds

doesn't hold

(due to Issue 1: rank deficient)

※with full-observations

(ex: 2D-NSE with the Sobolev space)

Further result for partial obs.

The projected PO + add. infl. with partial obs. for Lorenz 96

\text{Let } N_u = 3N \text{ for } N \in \mathbb{N}, H = [\phi_1, \phi_2, \phi_4, \phi_5, \phi_7, \dots]^\top \in \mathbb{R}^{N_y \times N_u},
\text{where } (\phi_i)_{i=1}^{N_u} \text{ is a standard basis of } \mathbb{R}^{N_u}.

Illustration (partial observations for Lorenz 96 equation)

○

○

×

○

○

×

○

○

×

○

○

×

○: observed

×: unobserved

u^1
u^2
u^3
u \in \mathbb{R}^{3N}
N_u = 3N, N_y = 2N

the special partial obs. for Lorenz 96

Further result for partial obs.

The projected PO + add. infl. with partial obs. for Lorenz 96

Theorem 3 (T.)

\limsup_{j \rightarrow \infty} \mathbb{E}[\| e^{(k)}_j \|^2] \le \frac{4 \min\{m, N_u\}}{1-\theta}r^2
\text{For } h \ll 1 \text{ and } \alpha \gg 1, \exist \theta = \theta(\beta, h, r, \alpha) \in (0, 1), \text{ s.t. }
\text{for } k = 1, \dots, m.
\text{Let } N_u = 3N \text{ for } N \in \mathbb{N}, H = [\phi_1, \phi_2, \phi_4, \phi_5, \phi_7, \dots]^\top \in \mathbb{R}^{N_y \times N_u},
\text{where } (\phi_i)_{i=1}^{N_u} \text{ is a standard basis of } \mathbb{R}^{N_u}, \text{ and let } \|\cdot\| = \sqrt{|\cdot|^2 + |H\cdot|^2}.
\text{Let } v_j^{(k)} \in B(\rho) \text{ for } k = 1, \dots, m, \text{ and } j \in \N.

The projected add. infl.:

\widehat{C}_j \rightarrow H^\top H(\widehat{C}_j + \alpha^2 I_{N_u})H^\top H

key 2

key 1

Summary

What's DA?

Sequential state estimation using model and data.

Results

Error analysis of a practical DA algorithm (ETKF)

under an ideal setting (T.+2024).

Discussion

Two issues with generalizing our analysis.

(rank, non-symmetric)

References

  • (T.+2024) K. T. & T. Sakajo, SIAM/ASA Journal on Uncertainty Quantification, 12(4), 1315–1335.
  • (T.) Kota Takeda, Error Analysis of the Ensemble Square Root Filter for Dissipative Dynamical Systems, PhD Thesis, Kyoto University, 2025.
  • (Kelly+2014) D. T. B. Kelly, K. J. H. Law, and A. M. Stuart (2014), Well-posedness and accuracy of the ensemble Kalman filter in discrete and continuous time, Nonlinearity, 27, pp. 2579–260.
  • (Al-Ghattas+2024) O. Al-Ghattas and D. Sanz-Alonso  (2024), Non-asymptotic analysis of ensemble Kalman updates: Effective dimension and localization, Information and Inference: A Journal of the IMA, 13.
  • (Tong+2016a) X. T. Tong, A. J. Majda, and D. Kelly (2016), Nonlinear stability and ergodicity of ensemble based Kalman filters, Nonlinearity, 29, pp. 657–691.
  • (Tong+2016b) X. T. Tong, A. J. Majda, and D. Kelly (2016), Nonlinear stability of the ensemble Kalman filter with adaptive covariance inflation, Comm. Math. Sci., 14, pp. 1283–1313.
  • (Kwiatkowski+2015) E. Kwiatkowski and J. Mandel (2015), Convergence of the square root ensemble Kalman filter in the large ensemble limit, Siam-Asa J. Uncertain. Quantif., 3, pp. 1–17.
  • (Mandel+2011) J. Mandel, L. Cobb, and J. D. Beezley (2011), On the convergence of the ensemble Kalman filter, Appl.739 Math., 56, pp. 533–541.

References

  •  (de Wiljes+2018) J. de Wiljes, S. Reich, and W. Stannat (2018), Long-Time Stability and Accuracy of the Ensemble Kalman-Bucy Filter for Fully Observed Processes and Small Measurement Noise, Siam J. Appl. Dyn. Syst., 17, pp. 1152–1181.

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

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

  • (Anderson 2001) J. L. Anderson (2001), An Ensemble Adjustment Kalman Filter for Data Assimilation, Mon. Weather Rev., 129, 2884–2903.

  • (Dieci+1999) L. Dieci and T. Eirola (1999), On smooth decompositions of matrices, Siam J. Matrix Anal. Appl., 20, pp. 800–819.

  • (Reich+2015) S. Reich and C. Cotter (2015), Probabilistic Forecasting and Bayesian Data Assimilation, Cambridge University Press, Cambridge.

  • (Law+2015) K. J. H. Law, A. M. Stuart, and K. C. Zygalakis (2015), Data Assimilation: A Mathematical Introduction, Springer.

Appendix

Problems on numerical analysis

  1. Discretization/Model errors
  2. Dimension-reduction for dissipative dynamics
  3. Numerics for Uncertainty Quantification (UQ)
  4. Imbalance problems

1. Discretization/Model errors

Motivation

- Including discretization errors in the error analysis.

- Consider more general model errors and ML models.

- Error analysis with ML/surrogate models.

cf.: S. Reich & C. Cotter (2015), Probabilistic Forecasting and Bayesian Data Assimilation, Cambridge University Press.

 

Question

- Can we introduce the errors as random noises in models?

- Can minor modifications to the existing error analysis work?

\widetilde{\Psi}_t(x) = \Psi_t(x) + \eta_t ?

2. Dimension-reduction for dissipative dynamics

Motivation

- Defining "essential dimension"                      and relaxing the condition:

-          should be independent of discretization.

cf.: C. González-Tokman & B. R. Hunt (2013), Ensemble data assimilation for hyperbolic systems, Physica D: Nonlinear Phenomena, 243(1), 128-142.

N_{ess} \le N_x
m > N_x \rightarrow m > N_{ess}
N_{ess}

Question

- Dissipativeness, dimension of the attractor？

- Does the dimension depend on time and state? (hetero chaos)

- Related to the SVD and the low-rank approximation?

Attractor

3. Numerics for UQ

Motivation

- Evaluating physical parameters in theory.

- Computing the Lyapunov spectrum of atmospheric models.

- Estimating unstable modes (e.g., singular vectors of Jacobian).

cf.: T. J. Bridges & S. Reich (2001), Computing Lyapunov exponents on a Stiefel manifold. Physica D, 156, 219–238.

F. Ginelli et al. (2007), Characterizing Dynamics with Covariant Lyapunov Vectors, Phys. Rev. Lett. 99(13), 130601.

Y. Saiki & M. U. Kobayashi (2010), Numerical identification of nonhyperbolicity of the Lorenz system through Lyapunov vectors, JSIAM Letters, 2, 107-110.

Question

- Can we approximate them only using ensemble simulations?

- Can we estimate their approximation errors?

4. Imbalance problems

Motivation

- Assimilating observations can break the balance relations inherited from dynamics. → numerical divergence

(e.g., the geostrophic balance and gravity waves）

- Explaining the mechanism (sparsity or noise of observations?).

- Only ad hoc remedies (Initialization, IAU).

cf.: E. Kalney (2003), Atmospheric modeling, data assimilation and predictability, Cambridge University Press.

Hastermann et al. (2021), Balanced data assimilation for highly oscillatory mechanical systems, Communications in Applied Mathematics and Computational Science, 16(1), 119–154.

Question

- Need to define a space of balanced states?

- Need to discuss the analysis state belonging to the space?

- Can we utilize structure-preserving schemes?

balanced states