Mathematical analysis of
data assimilation
Kota Takeda
(Kyoto University / RIKEN)
& related topics
Kota Takeda
PhD Student at Kyoto University
Junior Research Associate at RIKEN
Mathematics
Data Assimilation
Position
Kota Takeda 2
Activities
Organized MS at EASIAM 2024 Macao!
Advisor for high-school students
Volunteer at NPO clack
Cooperation with Leave a Nest Co., Ltd.
Teach programming
Assist students' research activities
See my website
Marathon
Finished Kobe marathon!
- Self-introduction ✓
- Data Assimilation/My research
- Related topics
This talk
Data Assimilation (DA)
The seamless integration of data and models
Model
Dynamical systems and its numerical simulations
ODE, PDE, SDE,
Numerical analysis, HPC,
Fluid mechanics, Physics
chaos
error
Data
Obtained from experiments or measurements
Statistics, Data analysis,
Radiosonde/Radar/Satellite observation
noise
limited
×
∂tu+(u⋅∇)u=−ρ1∇p+ν△u+f
\partial_t \bm{u} + (\bm{u} \cdot \nabla) \bm{u} = - \frac{1}{\rho} \nabla p + \nu \triangle \bm{u} + \bm{f}
Each has Uncertainty! (incomplete)
Data Assimilation (DA)
-
Control
- Nudging/Continuous Data Assimilation
-
Bayesian
- Particle Filter (PF)
- Ensemble Kalman filter (EnKF)
-
Variational
- Four-dimensional Variation method (4DVar)
Several formulations & algorithms of DA:
current study
→ next slide
Numerical Weather Prediction
typical application
×
Setup: Model
dtdu=F(u)
\frac{d u}{dt} = \mathcal{F}(u)
(F:U→U,u(0)=u0∈U)
(\mathcal{F}: \mathcal{U} \rightarrow \mathcal{U}, u(0) = u_0 \in \mathcal{U})
u(t)=Ψt(u0)
u(t) = \Psi_t (u_0)
Model dynamics in
Semi-group in
U
\mathcal{U}
U
\mathcal{U}
generates
U
\mathcal{U}
(Ψt:U→U,t≥0)
(\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
t
A small perturbation can
grow exponentially fast.
unknown orbit
Setup: Observation
Discrete-time model (at obs. time)
U
\mathcal{U}
(Ψ=Ψh,h>0: observation interval)
(\Psi = \Psi_h, h > 0 \text{: observation interval})
We have noisy observations in
t
t
yj=Huj+ξj,j∈N.
y_j = H u_j + \xi_j, \quad j \in \mathbb{N}.
(H:L(U,Y),ξj∈Y noise)
(H: \mathcal{L}(\mathcal{U}, \mathcal{Y}), \xi_j \in \mathcal{Y} \text{ noise})
Noisy observation
uj=Ψ(uj−1)
u_j = \Psi(u_{j-1})
observe
u1
u_1
y1
y_1
u2
u_2
y2
y_2
t1
t_1
t2
t_2
t3
t_3
...
h
h
bounded linear
∈Y
\in \mathcal{Y}
Y
\mathcal{Y}
(another Hilbert space).
∈Y
\in \mathcal{Y}
Aim of Data Assimilation
uj∈U
u_j \in \mathcal{U}
Unkown state
Given obs.
Yj={yi∣i≤j}⊂Y
Y_j = \{y_i \mid i \le j\} \subset \mathcal{Y}
vj
v_j
Construct
uj
u_j
estimate
Quality Assessment
E[∣vj−uj∣2]
\mathbb{E}\left[|v_j - u_j|^2\right]
yj
y_j
U
\mathcal{U}
Y
\mathcal{Y}
uj−1
u_{j-1}
uj
u_j
yj−1
y_{j-1}
Given
Model
Obs.
Obs.
unknown
(∣⋅∣: norm of U,E: expectation)
(|\cdot|: \text{ norm of } \mathcal{U}, \mathbb{E}: \text{ expectation})
(for j∈N)
(\text{for } j \in \mathbb{N})
(filtering problem)
Bayesian Data Assimilation
Pvj−1
\mathbb{P}^{v_{j-1}}
Pv^j
\mathbb{P}^{\hat{v}_j}
Pvj
\mathbb{P}^{v_j}
(I) Prediction
by Bayes' rule
by model
(II) Analysis
yj
y_j
We can compute it iteratively (2 steps at each time step).
Prediction
...
(⇒Pvj=Puj(⋅∣Yj))
(\Rightarrow \mathbb{P}^{v_j} = \mathbb{P}^{u_j}({}\cdot{} | Y_j))
U
\mathcal{U}
t
t
y1
y_1
t1
t_1
t2
t_2
t3
t_3
h
h
t0
t_0
(I) Prediction
Ψ
\Psi
(II)Analysis
Puj(⋅∣Yj)
\mathbb{P}^{u_j}({}\cdot{} | Y_j)
Yj={yi∣i≤j}: observations up to time j.
Y_j = \{y_i \mid i \le j\}: \text{ observations up to time } j.
Goal: compute the conditional distribution
y2
y_2
Need approximations!
Repeat (I) & (II)
→ next slide
- Two major implementations of EnKF:
- Perturbed Observation (PO): stochastic → sampling error
- Ensemble Transform Kalman Filter (ETKF): deterministic
Ensemble filters
U
\mathcal{U}
t
t
y1
y_1
y2
y_2
t1
t_1
t2
t_2
t3
t_3
...
h
h
t0
t_0
Repeat (I) & (II)
(II)Analysis
(I) Prediction
Ψ
\Psi
v0(1)
v_0^{(1)}
v0(m)
v_0^{(m)}
v0(2)
v_0^{(2)}
v1(1)
\widehat{v}_1^{(1)}
v1(1)
v_1^{(1)}
⋮
\vdots
Particle Filter(PF)
→ full distribution
Ensemble Kalman filter (EnKF)
→ mean & Covariance
(Burgers+1998)
(Bishop+2001)
A set of state vectors
(m→∞)
(m \rightarrow \infty)
Mathematical analysis of EnKF
Assumptions on model
(Reasonable)
Assumption on observation
(Strong)
We focus on finite-dimensional state spaces
(U=RNu for Nu∈N)
(\mathcal{U} = \R^{N_u} \text{ for } N_u \in \mathbb{N})
(T.+2024) K. T. & T. Sakajo (accepted 2024), SIAM/ASA Journal on Uncertainty Quantification.
Assumptions for analysis
Ny=Nu,H=INu,ξ∼N(0,r2INu)
N_y = N_u, H = I_{N_u}, \xi \sim \mathcal{N}(0, r^2 I_{N_u})
Assumption (obs-1)
"full observation"
Strong!
∃ρ>0 s.t. ∀u0∈B(ρ),∀t>0,u(t;u0)∈B(ρ).
\exist \rho > 0 \text{ s.t. } \forall u_0 \in B(\rho), \forall t > 0, u(t ; u_0) \in B(\rho).
∃β>0 s.t. ∀u∈B(ρ),∀v∈RNu,⟨F(u)−F(v),u−v⟩≤β∣u−v∣2.
\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.
∃β>0 s.t. ∀u,v∈B(ρ),∣F(u)−F(v)∣≤β∣u−v∣.
\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')
∀u0∈RNu,u˙(t)=F(u(t)) has an unique solution, and
\forall u_0 \in \R^{N_u}, \dot{u}(t) = F(u(t)) \text{ has an unique solution, and}
Reasonable
(B(ρ)={u∈RNu∣∣u∣≤ρ})
(B(\rho) = \{u \in \R^{N_u} \mid |u| \le \rho \})
"one-sided Lipschitz"
"Lipschitz"
"bounded"
or
Dissipative dynamical systems
Dissipative dynamical systems
Lorenz 63, 96 equations, widely used in geophysics,
dtdui=(ui+1−ui−2)ui−1−ui+f(i=1,…,Nu),
\frac{du^i}{dt} = (u^{i+1} - u^{i-2}) u^{i-1} - u^i + f \quad (i = 1, \dots, N_u),
Example (Lorenz 96)
u=(ui)i=1Nu∈RNu,f∈R,
\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
dtdu=
\frac{du}{dt} =
B(u,u)
B(u, u)
+
+
Au
A u
+
+
f.
f.
bilinear operator
linear operator
dissipative dynamical systems
Mimics chaotic variation of
physical quantities at equal latitudes.
(Lorenz1996) (Lorenz+1998)
Dissipative dynamical systems
The orbits of the Lorenz 63 equation.
Literature review: analysis of EnKF
H=INu
H = I_{N_u}
Error analysis for dissipative dynamical systems
-
Full observation ( )
-
PO + add. inflation (Kelly+2014)
-
ETKF + multi. inflation (T.+2024)
-
-
Partial observation → future
asmp obs1
asmp model1
asmp model2
asmp model2'
How EnKF approximate
the true state?
E∣error∣2≤O(robs2).
\mathbb{E}|error|^2 \le O(r_{obs}^2).
→ current
-
Consistency (Mandel+2011, Kwiatkowski+2015)
-
Sampling errors (Al-Ghattas+2024)
How EnKF approximate KF?
μjm,Cjm→μj,Cj
\mu^{m}_j, C^{m}_j \rightarrow \mu_j, C_j
Analysis for linear-Gaussian systems
Error analysis of ETKF (our)
Let ej=vj−uj where vj=m1k=1∑mvj(k).
\text{Let } e_j =\overline{v}_j - u_j\text{ where } \overline{v}_j = \frac{1}{m} \sum_{k=1}^m v^{(k)}_j.
Let us consider the ETKF with multiplicative inflation α≥1.
\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.
j→∞limsupE[∣ej∣2]≤1−θNur2+(1−θ1−Θ−1)D.
\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.
Theorem 2 (T.+2024)
where Θ=Θ(β,ρ,h,r,m,α),D=D(β,ρ)>0.
\text{where } \Theta = \Theta(\beta, \rho, h, r, m, \alpha) , D = D(\beta, \rho) > 0.
Then, ∃θ=θ(β,ρ,h,r,m,α)<1 s.t.
\text{Then, } \exist \theta = \theta(\beta, \rho, h, r, m, \alpha) < 1\text{ s.t.}
Let Assumptions (obs-1, model-1,2’) hold.
\text{Let Assumptions (obs-1, model-1,2') hold.}
Suppose m>Nu,rank(C0)=Nu,α≫1, and
\text{Suppose } m > N_u, \operatorname{rank}(C_0) = N_u,\alpha \gg 1, \text{ and}
vj(k)∈B(ρ) for k=1,…,m and j∈N.
v_j^{(k)} \in B(\rho) \text{ for } k = 1, \dots, m \text{ and } j \in \N.
Sketch: Error analysis of ETKF
<e−βh
< e^{-\beta h}
=1+α2r−2λmin(Cj)1
= \frac{1}{1 + \alpha^2 r^{-2} \lambda_{min}(\widehat{C}_j)}
①
②
Key result
∵
\because
(1) We establish a lower bound of eigenvalues of Cj by α≫1.
(1) \text{ We establish a lower bound of eigenvalues of } \widehat{C}_j \text{ by } \alpha \gg 1.
(2) If Cj is positive-difinite, we have a bound by α≫1.
(2) \text{ If } \widehat{C}_j \text{ is positive-difinite, we have a bound by } \alpha \gg 1.
Data Assimilation (DA)
-
Control
- Nudging/Continuous Data Assimilation
-
Bayesian
- Particle Filter (PF)
- Ensemble Kalman filter (EnKF)
- Quantum DA
-
Variational
- Four-dimensional Variation method (4DVar)
PDE
Optimal Transport
Metrics
Machine
Learning
Diffusion model
Back Prop.
Supervised Learning
Lyapunov
analysis
Signature
Wasserstein
PH
Operator
Algebra
Koopman operator
Optimization
Dynamical system
Numerical Analysis
structure-preserving
Numerical error
SDE
Dynamics
Uncertainty Quantification
Probability/Statistics
& related topic
Stochastic Galerkin
Random ODE
Other formulations
Control formulation
Problem
Keywords/Mathematics
(Olson+2003), (Azouani+2011)
Nudging/Continuous Data Assimilation
Squeezing property, Determining mode/node
AOT-algorithm
dtdu=F(u),
\frac{du}{dt} = \mathcal{F}(u),
v=Pu+q,dtdq=(I−P)F(v).
v = \mathcal{P} u + q, \quad \frac{dq}{dt} = (I - \mathcal{P})\mathcal{F}(v).
If we have noiseless & continuous time observations
through an orthogonal projection
→ copy obs. to the simulated state
(true)
(simulate)
P:U→Y.
\mathcal{P}: \mathcal{U} \rightarrow \mathcal{Y}.
Estimate the smallest rank of P s.t.
\text{Estimate the smallest rank of } \mathcal{P} \text{ s.t.}
t→∞lim∣v(t)−u(t)∣→0.
\lim_{t \rightarrow \infty}|v(t) -u(t)| \rightarrow 0.
PDE (Partial Differential Equations)
Control theory/Synchronization
Incomplete atmospheric data
(Charney+1969)
u=∣k∣≤λ∑ukϕk+∣k∣>λ∑ukϕk
u = \sum_{|k| \le \lambda} u_k \phi_k + \sum_{|k| >\lambda} u_k \phi_k
v=∣k∣≤λ∑ukϕk+∣k∣>λ∑vkϕk
v = \sum_{|k| \le \lambda} u_k\phi_k + \sum_{|k| >\lambda} v_k \phi_k
Pλu
\mathcal{P}_{\lambda} u
copy
converge
(t→∞)
(t\rightarrow\infty)
λ>0: threshold for wavelength
\lambda > 0: \text{ threshold for wavelength}
Variational formulation
Problem
Keywords/Mathematics
(Law+2015), (Korn+2009)
4DVar
Adjoint method to compute gradient
Optimization, Iteration method
Well-posedness of a mininization problem
\text{Well-posedness of a mininization problem}
v0=u0argminJ(u0;θ)
v_0 = \argmin_{u_0} J(u_0; \theta)
To estimate the initial state
from observation data
y1,⋯,yJ∈Y
y_1, \cdots, y_J \in \mathcal{Y}
J(u0;θ)=j=1∑J∥R−1/2(yj−H(Ψhj(u0;θ)))∥2
J(u_0;\theta) = \sum_{j=1}^J \|R^{-1/2} (y_j - H(\Psi_{hj}(u_0;\theta)))\|^2
u0∈U
u_0 \in \mathcal{U}
minimize a cost function.
∇u0J(u0)
\nabla_{u_0} J(u_0)
u0
u_0
y1
y_1
y2
y_2
y3
y_3
y4
y_4
Supervised Learning
(θ: model parameters)
(\theta: \text{ model parameters})
θ=θargminJ(u0;θ)
\widehat{\theta} = \argmin_{\theta} J(u_0; \theta)
Back propagation
DA & Numerical Analysis
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.
(Reich+2015)
Question
- Can we introduce the errors as random noises in models?
- Can minor modifications to the existing error analysis work?
Ψt(x)=Ψt(x)+ξt?
\widetilde{\Psi}_t(x) = \Psi_t(x) + \xi_t ?
Math: Numerical analysis, UQ
Dimension-reduction for dissipative dynamics
Motivation
- Defining "essential dimension" and relaxing the condition:
- should be independent of discretization.
(Tokman+2013)
Ness≤Nu
N_{ess} \le N_u
m>Nu→m>Ness
m > N_u \rightarrow m > N_{ess}
Ness
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
Math: Dynamical system
Numerics for UQ
Motivation
- Computing the Lyapunov spectrum of atmospheric models.
- Estimating unstable modes (e.g., singular vectors of Jacobian).
(Bridges+2001) (Ginelli+2007) (Saiki+2010)
Question
- Can we approximate them only using ensemble simulations?
- Can we estimate their approximation errors?
Math: Lyapunov analysis, Numerical analysis
Imbalance problems
Motivation
- Assimilating observations can break the balance relations inherited from dynamics. → numerical divergence
(e.g., geostrophic balance and gravity waves）
- Explaining the mechanism (sparsity or noise of observations?).
- Only ad hoc remedies (Initialization, IAU).
|
(Kalney 2003), (Hastermann+2021) |
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
Math: Structure-preserving scheme, constrained optimization
Related topics
Stochastic Galerkin
Problem
Keywords/Mathematics
(Sullivan2015) (Janjić+2023)
Weak formulation, Galerkin method, Lax-Milgram,
Generalized Polynomial Chaos Expansions
Karhunen–Lo`eve Expansions
Spectral expansion for sample space.
U(t,θ)=k∈N∑uk(t)ϕk(θ)
U(t, \theta) = \sum_{k \in \mathbb{N}} u_k(t) \phi_k(\theta)
θ(ω)∼N(0,1),
\theta (\omega) \sim N(0, 1),
E[U(t)]=u1(t),
\mathbb{E} [U(t)] = u_1(t),
V[U(t)]=k∈N∑uk2⟨ϕk2⟩.
\mathbb{V} [U(t)] = \sum_{k \in \mathbb{N}} u_k^2 \langle \phi_k^2 \rangle.
U(t,ω)
U(t, \omega)
Random field
is represented by
using a stochastic germ
and the associated polynomials
(ϕk(θ))k∈N.
(\phi_k(\theta))_{k \in \mathbb{N}}.
We can compute cumulants using the coefs.
U¨=−ΩU,U(0)=1,U˙(0)=0.
\ddot{U} = - \Omega U, U(0) = 1, \dot{U}(0) = 0.
logΩ∼N(0,σ2).
\log \Omega \sim N(0, \sigma^2).
Quantum DA & Operator Algebra
Problem
Keywords/Mathematics
Koopman operator
Generalizing data assimilation by embedding it in a nonabelian operator algebra based on von Neuman algebra
B=B(H).
\mathfrak{B} = B(H).
Quantum mechanics
Quantum computer
A=L∞(X,μ).
\mathfrak{A} = L^\infty(X, \mu).
(H=L2(X,μ))
(H = L^2(X, \mu))
S(A):
S(\mathfrak{A}):
Observable
π:A∋f↦πf∈B,
\pi: \mathfrak{A} \ni f \mapsto \pi f \in \mathfrak{B},
Embedding
(πf)g=fg∈H, ∀g∈H.
(\pi f) g = fg \in H, \ \forall g \in H.
State
continuous linear positive functionals on
A.
\mathfrak{A}.
Metrics on DA
Problem
Keywords/Mathematics
Persistent Homology
In the variational formulation, the analysis state can be understood as the barycenter of the prediction and observation.
(N. Feyeux et al. 2018)
Signature / Rough path
Wasserstein / OT
v=uargmind1(u,v)2+d2(Hu,y)2.
v = \argmin_u d_1(u, \widehat{v})^2 + d_2(Hu, y)^2.
(N. Feyeux et al. 2018)
Different ways to define a barycenter.
Limitation of the Euclid average.
where d1,d2: distance-like functions.
\text{where } d_1, d_2 \text{: distance-like functions.}
PF & Optimal Transport
Problem
Keywords/Mathematics
(Reich2013)
Optimal Transport, Sinkhorn algorithm
Monge-Kantorovich, Benamou-Brenier
u(1)
\widehat{u}^{(1)}
prior
posterior
Aim: obtain equal-weighted posterior samples
OT
Gradient flow
Mean-field, Fokker-Planck
Diffusion model
Schrödinger bridge
stochastic-ver.
Z={u^1,…,u^m}
\mathcal{Z} = \{\hat{u}_1, \dots, \hat{u}_m\}
Prior samples
Seek a coupling
between two discrete random variables
Za:Ω→Z
Z^a: \Omega \rightarrow \mathcal{Z}
Zf:Ω→Z
Z^f: \Omega \rightarrow \mathcal{Z}
pf=(1/m,…,1/m)T
p^f = (1/m, \dots, 1/m)^T
pa=(w1,…,wm)T
p^a = (w_1, \dots, w_m)^T
T∈Rm×m
T \in \mathbb{R}^{m \times m}
tij=(T)ij≥0.
t_{ij} = (T)_{ij} \ge 0.
i=1∑mtij=1/m,j=1∑mtij=wi
\sum_{i=1}^m t_{ij} = 1/m, \quad \sum_{j=1}^m t_{ij} = w_i
u(m)
\widehat{u}^{(m)}
TminEZf,Za[d(zf,za)]=Tmini,j∑mtijd(xif,xjf).
\min_T \mathbb{E}_{Z^f, Z^a}[d(z^f, z^a)] = \min_T \sum_{i, j}^m t_{ij} d(x^f_i, x^f_j).
⋮
\vdots
u(1)
u^{(1)}
u(m)
u^{(m)}
⋮
\vdots
(Jarrah+2023)
Mathematical analysis of data assimilation Kota Takeda Applied Mathematics Freshman Seminar 2024 , 2024/11/2-4, Kyoto University (Kyoto University / RIKEN) & related topics
Mathematical analysis of data assimilation and related topics
By kotatakeda
Mathematical analysis of data assimilation and related topics
Applied Mathematics Freshman Seminar 2024 Kyoto
- 137
