ノイズを含む部分観測に対するEnsemble Kalman Filterの誤差解析

2024年度応用数学合同研究集会 2024年12月5日@龍谷大学

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

The second author is supported by JST MIRAI JPMJMI22G1.

○竹田航太，坂上貴之

\dagger

京都大学理学研究科

\dagger

noisy

partial

イントロ：気象予報と観測

(a) 観測ノイズ：統計的な誤差が含まれる．

(b) 部分観測：全格子・全変数のデータが揃っていない．

(x,y,z)\mapsto(x,y)

ノイズ

部分観測

現実の大気から得られる観測データの特徴

気象庁

→ 気象予報シミュレーションの初期値として

精度・解像度が不十分．

イントロ：データ同化

(a) 観測ノイズ

(b) 部分観測

「現実」

観測値

y_n = H u_n + \eta_n

(a)

(b)

(u_n)_n

Lorenz 96

EnKF

ある観測行列

評価

観測

一般のモデルで

H = I

の場合の結果はある．

(Kelly+2014, T.+2024)

モデル

カオス

生成

u_n = \Psi(u_{n-1})

推定誤差

推定値

(v_n)_n

データ同化

高精度

高解像度

準備：モデル（一般）

\frac{d u}{dt} = \mathcal{F}(u)

(\mathcal{F}: \R^{N_u} \rightarrow \R^{N_u}, u(0) = u_0 \in \R^{N_u})

\Psi_t (u_0) = u(t; u_0)

モデル

半群

\R^{N_u} (N_u \in \mathbb{N})

生成

(\Psi_t : \R^{N_u} \rightarrow \R^{N_u}, t \ge 0)

: 状態空間

モデル：既知，初期値：未知

真の解

未知

準備：Lorenz 96方程式

\frac{du^i}{dt} = (u^{i+1} - u^{i-2}) u^{i-1} - u^i + F,

Lorenz 96 (L96)

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

非線形・エネルギー保存

線形・散逸

外力

地球の等緯度常に分布する物理量のカオスな変動を模した現象論的モデル．

(Lorenz1996)(Lorenz+1998)

(1)

(i = 1, \dots, N_u)

u^{0} = u^{N_u}, u^{-1} = u^{N_u - 1}, u^{N_u+1} = u^1.

※ 実数の絶対値とベクトルのユークリッドノルムの両方に

|\cdot|

を用いる．

補題1

F \in \R, u_0 \in \R^{N_u}

任意の

について，(1) の一意な解が存在する．

さらに，

\rho = \sqrt{2 N_u} |F|

に対して，

B(\rho) = \{v \in \mathbb{R}^{\Nu} \mid |v| \le \rho\}

は

次を満たす：

u_0 \in B(\rho) \Rightarrow u(t; u_0) \in B(\rho).

準備：観測（一般）

離散時間モデル

\mathbb{R}^{N_u}

\Psi = \Psi_\tau, \tau > 0

y_n = H u_n + \eta_n, \quad n \in \mathbb{N}.

H \in \mathbb{R}^{N_y \times N_u}, \eta_n \sim N(0, R), R \in \mathbb{R}^{N_y \times N_y}

観測

u_n = \Psi(u_{n-1})

観測

u_1

y_1

u_2

y_2

t_1

t_2

t_3

...

\tau

\in \mathbb{R}^{N_y}

\Psi

：観測時間間隔

(N_y \le N_u)

N_y

次元観測

観測行列

観測ノイズ

準備：L96の部分観測

Lorenz 96方程式に対して，以下の部分観測を考える．

N_u = 3N

(\phi_i)_{i=1}^{N_u} \subset \mathbb{R}^{N_u}

イラスト

○

○: 観測

×: 観測しない

u^1

u^2

u^3

u \in \mathbb{R}^{3N}

N_u = 3N, N_y = 2N

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}.

に対し

とする．

ただし，

は標準基底である．

← Key 1

(12)

H = [\phi_1, \phi_2, \phi_4, \phi_5, \phi_7, \dots]^\top

準備：時系列状態推定

\mathbb{R}^{N_u}

y_1

y_2

t_1

t_2

t_3

...

\tau

y_2

y_3

u_1

u_2

u_3

\Psi

given

Y_n = \{y_i \mid i \le n\}

問題

u_n

Estimate

n \in \mathbb{N},

known:

\Psi,

観測ノイズの分布

\mathbb{P}^{u_n}({}\cdot{} | Y_n) \approx \frac{1}{m} \sum_{k=1}^m \delta_{v^{(k)}_n}

アプローチ

条件付き分布を近似．

Ensemble Kalman filter (EnKF)

\mathbb{P}^{u_n}({}\cdot{} | Y_n) \approx \frac{1}{m} \sum_{k=1}^m \delta_{v^{(k)}_n}({}\cdot{})

に対して，条件付き分布を近似．

m \in \mathbb{N}

アンサンブル数

V_n = [v^{(k)}_n]_{k=1}^m

の時間発展を計算．

各

n \in \mathbb{N}

で以下の2ステップを行い

アンサンブル

\mathbb{P}^{u_n}({}\cdot{} | Y_n)

v^{(1)}_n,\dots,v^{(m)}_n

\widehat{v}^{(k)}_n = \Psi(v^{(k)}_{n-1}), \quad k = 1, \dots, m.

：各メンバーをモデルで時間発展．

(I) 予測

V_{n-1}

\widehat{V}_n = [\hat{v}^{(k)}_n]_{k=1}^m

予測アンサンブル

\widehat{V}_n

\rightarrow

(Burgers+1998)

\widehat{V}_n

V_{n-1}

V_n

y_n

(I)

(II)

y^{(k)}_n = y_n + \eta_n^{(k)}, \quad \eta_n^{(k)} \sim N(0, R), \quad k = 1, \dots, m.

\widehat{V}_n

v_n^{(k)} = (I - K_n H) \widehat{v}_n^{(k)} + K_n y_n^{(k)}, \quad k = 1, \dots, m.

K_n = \widehat{P}_n H^\top (H \widehat{P}_n H^\top + R)^{-1}, \widehat{P}_n = \operatorname{Cov}_m(\widehat{V}_n).

(\bigstar)

：まず，観測値を複製する．

次に予測値と観測値の重みつき平均を取る．

ただし，

(II) 解析

予測誤差共分散をサンプル共分散で近似

V_n

y_n

\rightarrow

EnKFの解釈

v^{(k)}_n = (I - K_n H) \widehat{v}^{(k)}_n + K_n y^{(k)}_n.

(\bigstar)

解析値は予測値と観測値の重み付き平均．

v^{(k)}_n

\widehat{v}^{(k)}_n

y^{(k)}_n

: 観測の不確実性．既知．

\widehat{P}_n

: 予測の不確実性．アンサンブルで近似→

\operatorname{Cov}_m(\widehat{V}_n)

K_n = \widehat{P}_n H^\top (H \widehat{P}_n H^\top + R)^{-1}

重み

は以下の「比」で決まる．

\widehat{P}_n = \widehat{p}_n^2, R = r^2, H = 1

イメージ（スカラーの場合）

とする．

v^{(k)}_n = \frac{r^2}{r^2 + \widehat{p}_n^2} \widehat{v}^{(k)}_n + \frac{\widehat{p}_n^2}{r^2 + \widehat{p}_n^2} y^{(k)}_n

(\bigstar)

= \frac{1}{1 + \widehat{p}_n^2/r^2} \widehat{v}^{(k)}_n + \frac{\widehat{p}_n^2/r^2}{1 + \widehat{p}_n^2/r^2} y^{(k)}_n.

共分散インフレーション

(I) 予測 → 　　　　　　　　→ (II) 解析 → ...

アンサンブル数が少ない時に状態推定を安定させるため，

(II)解析ステップの前に共分散行列を修正．

インフレーション

(\alpha \ge 0)

\widehat{P}_n \rightarrow \widehat{P}^\alpha_n = \widehat{P}_n + \alpha^2 I_{N_u}.

加法的インフレーション

既存の方法

→ 予測の不確実性を強制的に増大

射影加法的インフレーション

\widehat{P}_n \rightarrow \widehat{P}^\alpha_n = H^\top H (\widehat{P}_n + \alpha^2 I_{N_u}) H^\top H .

提案

(\alpha \ge 0)

直交射影行列

共分散インフレーションの背景

課題: アンサンブル数少ない時，誤差共分散行列が退化する．

(I) 予測 → 　　　　　　　　→ (II) 解析 → ...

m \le N_u \Rightarrow \operatorname{rank} \widehat{P}_n \le m-1 < N_u

解決策: (II)解析ステップの前に共分散行列を修正

\left(\widehat{P}_n = \frac{1}{m-1} \sum_{k=1}^m (\widehat{v}^{(k)}_n -\overline{\widehat{v}}_n) (\widehat{v}^{(k)}_n -\overline{\widehat{v}}_n)^\top \right)

インフレーション

(\alpha \ge 0)

\widehat{P}_n \rightarrow \widehat{P}^\alpha_n = \widehat{P}_n + \alpha^2 I_{N_u}.

加法的インフレーション

既存の方法

→ 行列のランクが回復

射影加法的インフレーション

\widehat{P}_n \rightarrow \widehat{P}^\alpha_n = H^\top H (\widehat{P}_n + \alpha^2 I_{N_u}) H^\top H .

提案

(\alpha \ge 0)

直交射影行列

データ同化

(a) 観測ノイズ

(b) 部分観測

モデル

「現実」

観測値

データ同化

推定値

推定誤差

カオス

生成

y_n = H u_n + \eta_n

u_n = \Psi(u_{n-1})

(a)

(b)

(u_n)_n

(v_n)_n

✓Lorenz 96

✓EnKF

✓ある観測行列

評価

(これから)

観測

一般のモデルで

H = I

の場合の結果はある．

(Kelly+2014, T.+2024)

主結果：誤差評価

各メンバーの誤差を

\delta^{(k)}_n = v_n^{(k)} - u_n \ (k=1, \dots, m)

とおき，

\|v\| = \sqrt{|v|^2 + |Hv|^2}

ノルム

を用いて評価する．

(r > 0)

定理1

\limsup_{n \rightarrow \infty} \mathbb{E}[\|\delta^{(k)}_n \|^2] \le \frac{4 N_y}{1-\theta}r^2

(k = 1, \dots, m).

v_n^{(k)} \in B(\rho) \quad (n \in \N, k = 1, \dots, m).

L96方程式 (1)，観測行列 (12)及び観測誤差共分散

R = r^2 I_{N_y}

に対し，EnKF及び射影加法的インフレーションを適用する．

任意の観測時間間隔

\tau > 0

アンサンブル数

m \in \mathbb{N}

インフレーションパラメータ

\alpha > 0

に対して，

以下を仮定する．

このとき，十分小さい

\tau > 0

と十分大きい

\alpha > 0

に対し，

\theta \in (0, 1)

が存在して，以下が成立．

予測誤差を抑える

予測の重みを小さく

観測ノイズの分散

主結果：誤差評価

定理1の詳細な評価

\mathbb{E}[\| \delta^{(k)}_n \|^2] \le \theta^n \mathbb{E}[\| \delta^{(k)}_0 \|^2] + \frac{4 N_y r^2 ( 1 - \theta^n)}{1-\theta},

(k = 1, \dots, m).

各メンバーの誤差を

\delta^{(k)}_n = v_n^{(k)} - u_n

とおき，

\|v\| = \sqrt{|v|^2 + |Hu|^2}

ノルム

を用いて評価する．

十分小さい

\tau > 0

と十分大きい

\alpha > 0

に対し，

\theta \in (0, 1)

が存在して，以下が成立．

exponential decay

initial error

定理1の仮定の下で

主結果：証明の流れ

定理1

射影インフレーション

EnKF

部分観測 (12)

Lorenz 96方程式

誤差拡大(補題2)

(I)予測ステップ

誤差縮小(補題3)

(II)解析ステップ

Key1

Key2

誤差のノルム

．．．

誤差拡大(補題2)

誤差縮小(補題3)

(I)予測

(II)解析

主結果：予測ステップでの誤差拡大

L96方程式に対して，良い部分観測 (12)を考えると以下の誤差拡大評価を得る．

補題2

観測行列 (12)を考える．

u_0, v_0 \in B(\rho), t \ge 0

について，

\delta(t) = \Psi_t(v_0) - \Psi_t(u_0)

とおく．このとき，

t \ge 0

で以下が成立．

は予稿を参照．

A_1(\tau), B_1(\tau), B_2(\tau)

(15)

(16)

注意：

A_1(\tau), B_1(\tau), B_2(\tau) < 1

を十分小さくと取ると

t = \tau > 0

とできる．

※ 補題2の証明 → Gronwallの不等式に基づく．

(20)

Key1

主結果：解析ステップでの誤差縮小

EnKFに射影加法的インフレーションを適用すると観測成分での誤差縮小を得る．

\Theta = \Theta(\alpha)

注意：

は

\alpha

を大きく取ることで，

\Theta(\alpha) \searrow 0 \quad (\alpha \rightarrow \infty)

いくらでも小さくできる．

補題3

定理1の仮定の下，

\Theta = \left(\frac{r^2}{r^2 + \alpha^2}\right)^2

とおく．

について，以下が成立．

k = 1, \dots, m

このとき，

\mathbb{E}[|H\delta^{(k)}_n|^2] \le \Theta \mathbb{E}[|H\widehat{\delta}^{(k)}_{n-1}|^2] + 2N_y r^2,

\mathbb{E}[|\delta^{(k)}_n|^2] \le \mathbb{E}[|\widehat{\delta}^{(k)}_{n-1}|^2] + 2N_y r^2.

(23)

(24)

(25)

Key2

予測誤差

観測誤差

縮小

※ 作用素ノルムの評価で射影が重要．

主結果：解析ステップでの誤差縮小

補題3の証明

\mathbb{E}[|H\delta^{(k)}_n|^2] \le \Theta \mathbb{E}[|H\widehat{\delta}^{(k)}_{n-1}|^2] + 2N_y r^2,

\mathbb{E}[|\delta^{(k)}_n|^2] \le \mathbb{E}[|\widehat{\delta}^{(k)}_{n-1}|^2] + 2N_y r^2.

(23)

(24)

(23)の右辺第１項，(24)右辺第2項の評価にそれぞれ

|(I_{N_y} + HG_n)^{-1}|_{op} = \Theta^{1/2},

|(I_{N_u} + G_n)^{-1}G_n|_{op} \le 1

という評価を用いる．

(27)

(26)

G_n = r^{-2}\widehat{P}_n^{\alpha} H^\top H

ただし，

である．

(目標)

誤差縮小

観測ノイズの影響

(27)では，

G_n

の対称半正定値性が重要．

→ インフレーション時の射影から従う．

(26)では，

G_n

の最小固有値が正であることが重要．

→ 加法的インフレーションから従う．

\lambda_{min}(G_n) = \lambda_{min}(r^{-2}H^\top H (\widehat{P}_n + \alpha^2 I) H^\top H) = r^{-2}\alpha^2 > 0

|\cdot|_{op}

は作用素ノルム．

主結果：証明の流れ

定理1

射影インフレーション

EnKF

部分観測 (12)

Lorenz 96方程式

誤差拡大(補題2)

誤差縮小(補題3)

誤差のノルム

．．．

誤差拡大(補題2)

でコントロール

誤差縮小(補題3)

でコントロール

\alpha

\tau

既存研究との比較

(a) ＼ (b)	全観測	部分観測
ノイズなし	-	Synchronization [6, 7]
ノイズあり	EnKF [3, 4]	EnKF (本研究)

(a) 観測ノイズ：統計的な誤差が含まれる．

(b) 部分観測：全格子・全変数のデータが揃っていない．

モデルは「一般」

[3] (Kelly+2014) [4] (T.+2024) [6] (Law+2016) [7] (Azouani+2014)

L96方程式, Navier-Stokes方程式 (NSE)

L96方程式

※ Synchronization

ノイズのない部分観測から時間漸近的に真の解を再構成

\mathcal{F} \circlearrowright u(t) \overset{H u(t)}{\rightarrow} v(t) \circlearrowleft \mathcal{F} \Rightarrow |u(t) - v(t)| \rightarrow 0 \ (t \rightarrow \infty)

観測の設定と誤差評価

予想

拡張

Synchronizationの既存結果からHを使う

Lorenz 63 (Hayden+2011), 2D-NSE (Aouzani+2014)

抽象論

散逸力学系 → 有限次元アトラクター，inertial manifold

→ 有限ランクの観測作用素 Hの存在

与えられた観測に対して時間大域評価が可能か判定することや，観測に対する制約のもとで最低ランクのものを見つけることは難しい

一方で，

(Hayden+2011)

数値結果: L96 (定理1の検証)

② add, projで挙動は変わらない．→ addでも誤差評価の可能性．

EnKFに2種類のインフレーションを3種類ので適用．

\alpha

加法的(add)，射影加法的(proj)

その他パラメータ：

N_u = 60, N_y = 40, \tau = 0.01, r=1, m=10.

① 誤差は

\alpha=0.0 > 2.0 > 0.5.

→ 　が大きければ良いわけではない．

\alpha

\alpha = 0.0

\alpha = 2.0

\alpha = 0.5

定理1の評価

観察

数値結果: Lorenz 63方程式

\left\{ \begin{array}{cc} \dot{x} & = -\sigma x + \sigma y,\\ \dot{y} & = \varrho x - y - xz,\\ \dot{z} & = - bz + xy. \end{array} \right.

典型的なパラメータは

\sigma = 10, b = 8/3, \varrho = 28.

3変数の常微分方程式

まとめ

1. Synchronizationが成り立つような良い部分観測を取り，

　予測誤差共分散を射影することで，

　EnKFの時間大域的な誤差評価を得た．

2. NSEなど他の散逸力学系でも同様の方針で示せそう．

　射影がない加法的インフレーションでも示せるかも．

近い将来

3. 超高次元の問題への応用上は，メモリ負荷の理由から

　加法的インフレーションが非現実的

　→ この場合はより精密な共分散行列の解析が必要．

さらにその先．．．

Key1: L96の部分観測 → 誤差拡大のコントロール

Key2: 射影加法的インフレーション → 行列の対称性

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.

References

(Law+2016) K. J. H. Law et al., Filter accuracy for the Lorenz 96 model: Fixed versus adaptive observation operators, Phys. Nonlinear Phenom., 325, 1–13.
(Azouani+2014) A. Azouani, E. Olson, and E. S. Titi, Continuous Data Assimilation Using General Interpolant Observables, J. Nonlinear Sci., 24, 277–304.
(Hayden+2011) K. Hayden, E. Olson, and E. S. Titi, Discrete data assimilation in the Lorenz and 2D Navier-Stokes equations, Phys. -Nonlinear Phenom., 240, pp. 1416– 1425.