データ同化の数学解析と
数値解析的な課題
○竹田 航太 (京都大学) , 坂上 貴之 (京都大学)
The first author is supported by RIKEN Junior Research Associate Program and JST SPRING JPMJSP2110.
The second author is supported by JST MIRAI JPMJMI22G1.
Data Assimilation (DA)
Seamless integration of data into numerical 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}
DA
Bayesian
Variational
Control
Filtering
EnKF
...
...
...
←Today
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 in
Observation in
\mathbb{R}^{N_x}
\mathbb{R}^{N_y}
(h > 0 \text{: obs. interval})
unknown
(t_j = jh)
time
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.
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 method
- 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
\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.
(I) X_{j-1} \rightarrow \widehat{X}_j: \text{ Evolve each member by model dynamics } \Psi,
(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.
Stochastic!
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,}\\
x_j^{(k)} = (I - K_j H) \widehat{x}_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}(\widehat{X}_j),
(\bigstar)
Perturbed Observation
ETKF (deterministic EnKF)
\widehat{x}^{(k)}_j = \Psi(x^{(k)}_{j-1}), \quad k = 1, \dots, m.
(I) X_{j-1} \rightarrow \widehat{X}_j: \text{ Evolve each member by model dynamics } \Psi,
(Bishop+2001)C. H. Bishop, B. J. Etherton, and S. J. Majumdar, Adaptive Sampling with the Ensemble Transform Kalman Filter. Part I: Theoretical Aspects, Mon. Weather Rev., 129 (2001), pp. 420–436.
Deterministic!
\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.\\
\overline{x}_j = (I - K_j H)\overline{\widehat{x}}_j + K_j y_j , \quad dX_j = d\widehat{X}_j T_j.
\text{Then, update separately,}
\frac{1}{m-1} d\widehat{X}_j T_j (d\widehat{X}_j T_j)^* = (I - K_j H) \widehat{C}_j.
The explicit representation exists!
(\bigstar)
Ensemble Transform Kalman Filter
Inflation (numerical technique)
Issue: Underestimation of covariance through the DA cycle
→ Poor state estimation → (idea) inflating covariance before step (II)
(I) Prediction → inflation → (II) Analysis → ...
Rank is improved
Rank is not improved
Two basic methods of inflation
\text{Additive } (\alpha \ge 0)
\text{Multiplicative } (\alpha \ge 1)
\text{ETKF}: d\widehat{X}_j \rightarrow \alpha d\widehat{X}_j.
\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.
\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_{N_x}.
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 (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.
(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.
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 \})
"one-sided Lipschitz"
"Lipschitz"
"bounded"
Dissipative dynamical system
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}^{N_x} \in \mathbb{R}^{N_x}, 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{dx}{dt} =
B(x, x)
+
A x
+
f.
bilinear operator
linear operator
dissipative dynamical systems
Dissipative dynamical systems
The orbits of the Lorenz 63 equation.
Error analysis of PO (prev.)
\text{Let us consider the PO method with additive inflation } \alpha > 0.
\text{Let } e^{(k)}_j = x^{(k)}_j - x^\dagger_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_x}
Strong effect! (improve rank)
Theorem 1 (Kelly+2014)
\limsup_{j \rightarrow \infty} \mathbb{E} [|e^{(k)}_j|^2] \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 Assumptions (obs-1, model-1,2) hold and }\alpha \gg 1.
\text{for } k = 1, \dots, m.
Error analysis of ETKF (our)
\text{Let } e_j =\overline{x}_j - x^\dagger_j\text{ where } \overline{x}_j = \frac{1}{m} \sum_{k=1}^m x^{(k)}_j.
\text{Let us consider the ETKF with multiplicative inflation } \alpha \ge 1.
(T.+2024) K. T. & T. Sakajo (accepted 2024), SIAM/ASA Journal on Uncertainty Quantification.
\limsup_{j \rightarrow \infty} \mathbb{E}[|e_j|^2] \le \frac{N_x r^2}{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_x, \operatorname{rank}(C_0) = N_x,\alpha \gg 1, \text{ and}
x_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{x}}_j - x^\dagger_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{x}_j = (I - K_j H)\overline{\widehat{x}}_j + K_j y_j
e_j = (I - K_j H) \widehat{e}_j + K_j (y_j - H x^\dagger_j).
subtracting yields
(II)
From the Kalman update relation
(\bigstar)
x^\dagger_j
prediction error
② contraction
③ obs. noise effect
estimated by op.-norm
Sketch: Error analysis of ETKF (our)
Key relations in
Using Woodbury's lemma
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
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.
Summary 1
\limsup_{j \rightarrow \infty} \mathbb{E}[|e^{(k)}_j|^2] \le \frac{2\min\{m, N_x\}}{1-\theta}r^2,
\limsup_{j \rightarrow \infty} \mathbb{E}[|e_j|^2] \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
due to
Further/Future
This study
1. Partial observation
2. Small ensemble
m \le N_x
N_y < N_x, H \in \mathbb{R}^{N_y \times N_x}
3. Infinite-dimension
\text{state space } \mathcal{X}: \text{Hilbert space},
\text{observation space } \mathcal{Y} \subset \mathcal{X}.
Further/Future
\|(I + \widehat{C}_j H^\top R^{-1} H)^{-1}\|_{op} \ge 1.
\text{If } m \le N_x \text{ or } N_y < N_x, H \in \mathbb{R}^{N_y \times N_x}, 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
Further/Future
\|(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_x, H \in \mathbb{R}^{N_y \times N_x}, 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)
Further/Future
Infinite-dimension
Theorem (Kelly+2014) for PO + add. infl.
Theorem (T.+2024) for ETKF + multi. infl.
\text{state space } \mathcal{X}: \text{Hilbert space},
\text{observation space } \mathcal{Y} \subset \mathcal{X} \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
Theorem 3 (T.)
\limsup_{j \rightarrow \infty} \mathbb{E}[\| e^{(k)}_j \|^2] \le \frac{4\min\{m, N_x\}}{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_x = 3N \text{ for } N \in \mathbb{N}, H = [\phi_1, \phi_2, 0, \phi_4, \phi_5, 0, \phi_7, \dots, 0] \in \mathbb{R}^{N_x \times N_x},
\text{where } (\phi_i)_{i=1}^{N_x} \text{ is a standard basis of } \mathbb{R}^{N_x}, \text{ and let } \|\cdot\| = \sqrt{|\cdot|^2 + |H\cdot|^2}.
\text{Let } x_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(\widehat{C}_j + \alpha^2 I_{N_x})H^*
key 2
key 1
Further/Future
This study
Problems on
numerical analysis
other directions
1. Partial observation
2. Small ensemble
m \le N_x
N_y < N_x, H \in \mathbb{R}^{N_y \times N_x}
3. Infinite-dimension
\text{state space } \mathcal{X}: \text{Hilbert space},
\text{observation space } \mathcal{Y} \subset \mathcal{X}.
Issue 1, 2
Issue 1
Issue 1
Problems on numerical analysis
- Discretization/Model errors
- Dimension-reduction for dissipative dynamics
- Numerics for Uncertainty Quantification (UQ)
- 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) + \xi_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
4. Imbalance problems
|
Initialization （モデルを安定化する方法）→ 情報を落とすので使いたくない (1) Nonlinear Normal Mode Initialization (NNMI): C.E. Leith (1980), J. Atmos. Sci., 37, 958–968. (2) digital filter initialization (DFI): モデルにフィルターを入れる． その他数値テクニック (3) 鉛直方向だけ陰解法 (4) split explicit: dtを項ごとに変える（音波だけ細かくする等） 最近のDAでは使われていない - 4DVarではDF由来のコスト関数を入れられる - EnKFではモデルを破綻させるモードは出ないはず（subspace property） - EnKFでは，localizationで問題が起きるかも？（理論的に示されている : J. D. Kepert (2009), Q.J.R. Meteorol. Soc., 135: 1157-1176. ） |
Summary 2
1. Partial obs.
Future plans
2. Small ens.
3. Infinite-dim.
Mathematical analysis of EnKF
rank deficient
non-symmetric
DA × Numerical analysis
Model errors
Dimension-reduction
Numerical UQ
Imbalance
Other DA
Particle Filter
4DVar
データ同化の数学解析と数値解析的な課題
By kotatakeda
データ同化の数学解析と数値解析的な課題
RIMS数値解析 2024/10/24
