#概要
複雑な（またはカオス的な）偏微分方程式を解くと、初期値の設定によって予測結果が大きく異なる場合があります。昔から気象などの分野では、天気予報の精度を上げるため、観測データを使い、データ同化の手法（変分法）を用いて数値モデルの修正を行なっています。
近年では、IoTやビックデータ活用などによりデータ同化が注目を集めています。データ同化は、大きく２つの手法があり、カルマンフィルタと変分法があります。
カルマンフィルタの特徴は、観測されたデータをもとに、逐次的に数値モデルを修正する方法です。変分法の特徴は、数値モデルが観測結果と一致するように、初期値や境界条件を修正する方法です。
今回の記事では、カルマンフィルタの基礎となる最適内挿法と３次元変分法について紹介します。４次元変分法との違いは、方程式が時間発展するかしないかの違いです。
簡単な推定問題
変数$x$を、ある時刻と位置における物理量とする。２つの値（観測値または予測値などの値）を$x_1,x_2$とし、誤差を$\epsilon_1,\epsilon_2$とする。変数$x$の真の値を$x^t$とし、この値に近づくように推定値$x^a$を求めることを考える。もちろん、真の値$x^t$と誤差$\epsilon$は未知である。また$x_1,x_2$は
\begin{align}
x_1=x^t+\epsilon_1 \\
x_2=x^t+\epsilon_2 \\
\end{align}
と書ける。
推定値$x^a$を求めるために以下２つの条件を課す。
（１）$x_1,x_2,x^a$は不偏推定量である。
平均の不偏推定量とは、母平均の平均と標本の平均が一致することである。つまり、
\begin{align}
\langle x_1\rangle =x^t \rightarrow \langle\epsilon_1 \rangle =0 \\
\langle x_2\rangle =x^t \rightarrow \langle\epsilon_2 \rangle =0 \\
\end{align}
を意味する。また、推定値$x^a$も不偏推定値とする。
\begin{align}
\langle x^a\rangle =x^t
\end{align}
（２）$x_1$と$x_2$の誤差は無相関である。
つまり、
\begin{align}
\langle (x_1-x^t)(x_2-x^t) \rangle =\langle\epsilon_1 \epsilon_2\rangle =0 \\
\end{align}
である。また$x_1$と$x_2$の誤差分散を$\sigma_1=\langle\epsilon_1^2\rangle,\sigma_2=\langle\epsilon_2^2\rangle$とする。
推定値$x^a$を$x_1,x_2$を使い
\begin{align}
x^a=ax_1+bx_2
\end{align}
と表す。
まずは、推定値$x^a$の平均を計算する。推定値$x^a$は不偏推定値なので
\begin{align}
\langle x^a\rangle &=x^t\\
&= a\langle x_1\rangle +b\langle x_2\rangle\\
&= (a+b)x^t\\
\mbox{よって}&\Rightarrow b=1-a
\end{align}
となる。次に、推定値$x^a$の誤差分散を計算する。
\begin{align}
\sigma_a^2 &=\left\langle (x^a-x^t)^2 \right\rangle \\
&= \left\langle \left(ax_1+(1-a)x_2 -x^t\right)^2 \right\rangle \\
&=a^2\left\langle \epsilon_1^2 \right\rangle +a(1-a)\left\langle \epsilon_1 \epsilon_2 \right\rangle +(1-a)^2\left\langle \epsilon_2^2 \right\rangle \\
&=a^2\sigma_1^2 +(1-a)^2\sigma_2^2 \\
\end{align}
最後の行は、条件（２）を使った。この式を平方完成して、誤差分散$\sigma_a^2 $が最小となる条件を求めると
\begin{align}
\sigma_a^2 &=(\sigma_1^2+\sigma_2^2)a^2 - 2\sigma_2^2 a+\sigma_2^2 \\
&= (\sigma_1^2+\sigma_2^2) \left(a^2 - \frac{2\sigma_2^2}{\sigma_1^2+\sigma_2^2} a+\frac{\sigma_2^2}{\sigma_1^2+\sigma_2^2} \right) \\
&= (\sigma_1^2+\sigma_2^2) \left(a - \frac{\sigma_2^2}{\sigma_1^2+\sigma_2^2} \right)^2 +\frac{\sigma_1^2\sigma_2^2}{\sigma_1^2+\sigma_2^2} \\
\end{align}
となり、誤差分散$\sigma_a^2 $が最小となる条件は
\begin{align}
a = \frac{\sigma_2^2}{\sigma_1^2+\sigma_2^2}
\end{align}
である。最終的に推定値$x^a$と誤差分散$\sigma_a^2 $は
\begin{align}
x^a&=x_1+\frac{\sigma_1^2}{\sigma_1^2+\sigma_2^2}(x_2-x_1)\\
\sigma_a^2&= \frac{\sigma_1^2\sigma_2^2}{\sigma_1^2+\sigma_2^2}
\end{align}
と求まる。
最適内挿法
３次元変分法を説明する前に、最適内挿法について説明する。最適内挿法と3次元変分法は関係がることが後々わかる。
先ほど求めた式は、
\begin{align}
x^a&=x_1+\frac{\sigma_1^2}{\sigma_1^2+\sigma_2^2}(x_2-x_1)
\end{align}
であるが、$x_1$を予測値として、$x_2$を観測値とすれば、数値モデルの予測値$x_1$を右辺の２項目が修正量と解釈できる。この式を拡張して、$x^a$を推論値、$x^b$を予測値、$y$を観測値として
\begin{align}
x^a_i&=x_i^b+W_{ij}\Delta y_j\\
\Delta y_j &= y_j - H_{jk}x_k^b
\end{align}
とかく。$H_{ij}$は観測行列であり、モデルによって決まる行列である。推論値$x^a$を求めるには行列$W_{ij}$を求める必要がある。$W_{ij}$は、簡単な推定問題で$a,b$を求めた手順で計算できる。
最初に、$x_i^t$を真の値とし、予測誤差$\epsilon^b_i$と観測誤差$\epsilon^o_i$の平均はゼロとし（条件１）、無相関とする（条件２）。
\begin{align}
\langle x_i^b-x_i^t \rangle &=\langle\epsilon_i^b\rangle =0 \\
\langle x_i^o - x_i^t \rangle &=\langle \epsilon_i^o\rangle =0 \\
\langle\epsilon_i^b \epsilon_j^o\rangle &=0
\end{align}
ただし、誤差分散行列は、対称行列であり
\begin{align}
B_{ij}=\langle\epsilon_i^b\epsilon_j^b\rangle \ \ ,\ \ R_{ij}=\langle\epsilon_i^o \epsilon_j^o\rangle
\end{align}
と与えられるとする。
解析誤差 $\bf{\epsilon}^a$ は、真の値 $\bf{x}^t$ を使い
\begin{align}
\bf{\epsilon}^a &= \bf{x}^a -\bf{x}^t \\
&= \left(\bf{x}^b - \bf{x}^t \right) +\bf{W}\left(\bf{y}-\bf{H}\bf{x}^b\right) \\
&= \left(\bf{x}^b - \bf{x}^t \right) +\bf{W}\left(\left(\bf{y}-\bf{H}\bf{x}^t \right)-\bf{H}\left(\bf{x}^b - \bf{x}^t \right)\right) \\
&= \bf{\epsilon}^b+\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right)
\end{align}
と書ける。
次に、解析誤差共分散行列 $\bf{P}^a$ を求めると、条件１と条件２を使い
\begin{align}
\bf{P}^a&= \left\langle \bf{\epsilon}^a \left(\bf{\epsilon}^a \right)^T \right\rangle \\
&= \left\langle \left[\bf{\epsilon}^b+\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right) \right] \left[\bf{\epsilon}^b+\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right) \right]^T \right\rangle \\
&=\left\langle \bf{\epsilon}^b \left(\bf{\epsilon}^b \right)^T + \bf{\epsilon}^b\left(\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right) \right)^T +\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right) \left(\bf{\epsilon}^b \right)^T +\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right)\left(\bf{W}\left(\bf{\epsilon}^o-\bf{H}\bf{\epsilon}^b\right) \right)^T \right\rangle \\
& = \left\langle \bf{\epsilon}^b \left(\bf{\epsilon}^b \right)^T\right\rangle - \left\langle \bf{\epsilon}^b \left(\bf{\epsilon}^b \right)^T\right\rangle \bf{H}^T\bf{W}^T - \bf{W}\bf{H}\left\langle \bf{\epsilon}^b \left(\bf{\epsilon}^b \right)^T\right\rangle \\
& + \bf{W}\left\langle \bf{\epsilon}^o \left(\bf{\epsilon}^o\right)^T\right\rangle \bf{W}^T + \bf{W}\bf{H}\left\langle \bf{\epsilon}^b \left(\bf{\epsilon}^b \right)^T\right\rangle \bf{H}^T\bf{W}^T \\
&=\bf{B} - \bf{B} \bf{H}^T\bf{W}^T - \bf{W}\bf{H}\bf{B} \\
& + \bf{W}\bf{R}\bf{W}^T + \bf{W}\bf{H}\bf{B} \bf{H}^T\bf{W}^T
\end{align}
と書ける。また$\bf{P}^a=(\bf{P}^a)^T$であり対称行列である。また解析誤差分散は、解析誤差共分散行列の対角成分 $\mathrm{Tr} \bf{P}^a = P^a_{ii}$である。
最後に、解析誤差分散$P^a_{ii}$の最小となる$\bf{W}$を求める。$P^a_{ii}$を$W_{pq}$で微分して
\begin{align}
\frac{\partial P^a_{ii}}{\partial W_{pq}}&= -B_{ik}H_{mk}\delta_{ip}\delta_{mq}+\delta_{ip}\delta_{kq}R_{km}W_{im} +W_{ik}R_{km} \delta_{ip}\delta_{qm} - \delta_{ip}\delta_{qk}H_{km}B_{mi} +\delta_{ip}\delta_{kq}H_{km}B_{mn}H_{ln}W_{il}+W_{ik}H_{km}B_{mn}H_{ln}\delta_{lq}\delta_{ip} \\
&= -2B_{pm}H_{qm}+2W_{pk}R_{kq}+2W_{pk}H_{km}B_{mn}H_{qm}\\
&=0
\end{align}
$\bf{W}$は、
\begin{align}
\bf{W} = BH^T \left(R + HBH^T \right)^{-1}
\end{align}
となる。最終的に、推論値$x^a$は
\begin{align}
\bf{x}^a&=\bf{x}^b +W(y-Hx^b)\\
&= \bf{x}^b+BH^T \left(R+HBH^T \right)^{-1} \left(y-Hx^b \right)
\end{align}
と求まる。
後々のために$\bf{W}$を
\begin{align}
\bf{W} &= \left(\bf{B}^{-1}+H^TR^{-1}H \right)^{-1}\left(\bf{B}^{-1}+H^TR^{-1}H \right) \bf{B}H^T \left(R + HBH^T \right)^{-1}\\
&=\left(\bf{B}^{-1}+H^TR^{-1}H \right)^{-1}\left(\bf{H}^T+H^TR^{-1}HBH^T \right) \left(\bf{R} + HBH^T \right)^{-1}\\
&=\left(\bf{B}^{-1}+H^TR^{-1}H \right)^{-1}\bf{H}^TR^{-1} \left(\bf{R} + HBH^T \right) \left(\bf{R} + HBH^T \right)^{-1}\\
&=\left(\bf{B}^{-1}+H^TR^{-1}H \right)^{-1}\bf{H}^TR^{-1}
\end{align}
と変形しておく。すると$x^a$は、
\begin{align}
\bf x^a = x^b + \left(B^{-1} - H^T R^{-1}H \right)^{-1} H^T R^{-1}(y-Hx^b)
\end{align}
である。
３次元変分法
予測値$x^b$と観測値$y$は、ガウス分布に従うとする。
\begin{align}
P_b\left(\bf x^b|x\right)&=\frac{1}{\left(\sqrt{2\pi } \right)^n \left|\bf B \right|^{1/2} } \exp\left\{-\frac{1}{2} \bf \left(x-x^b \right)^T B^{-1} \left(x-x^b \right) \right\} \\
P_o\left(\bf y|x\right)&=\frac{1}{\left(\sqrt{2\pi } \right)^m \left|\bf R \right|^{1/2} } \exp\left\{-\frac{1}{2} \bf \left(H(x)-y \right)^T R^{-1} \left(H(x)-y \right) \right\}
\end{align}
$n,m$は、格子点・観測点の数である。$x$は初期値であり、この２つのガウス分布を使い、尤度を
\begin{align}
\mathcal{L}\left(\bf x|x^b,y \right) = P( \bf x^b,y|x ) = P_b\left( \bf x^b|x \right)P_o\left( \bf y|x \right)
\end{align}
とする。尤度を使い最適な初期値$x$を求めることを考える。つまり、予測値と観測値の確率分布が最大となる初期値$x$を最尤推定により求める。
最尤推定は、最大となる尤度を求めることであるが、便利上、対数を取りマイナスをつけた量を誤差関数$J(x)$として
\begin{align}
J(x)&=-\log\mathcal{L}\left(\bf x|x^b,y \right) \\
&= \frac{1}{2} \bf \left(x-x^b \right)^T B^{-1} \left(x-x^b \right) + \frac{1}{2} \bf \left(H(x)-y \right)^T R^{-1} \left(H(x)-y \right) +\ldots
\end{align}
この誤差関数$J(x)$が最小値なる初期値$x$を求める。また$\ldots$は$J(x)$の最小化に関係ない部分である。
誤差関数$J(x)$を微分すると、$\bf B,\bf R$の対称性から
\begin{align}
\bf \nabla J(x) &= \bf B^{-1} \left(x-x^b \right) + \bar{H}^{T} R^{-1} \left(H(x)-y \right)
\end{align}
と書ける。また$\bf \bar{H}$は、観測演算子$H(x)$の微分であり
\bar{H}=\nabla H(x)=
\begin{pmatrix}
\frac{\partial H_1}{\partial x_1} & \frac{\partial H_1}{\partial x_2}& \cdots & \frac{\partial H_1}{\partial x_n} \\\
\frac{\partial H_2}{\partial x_1} & \frac{\partial H_2}{\partial x_2}& \cdots & \frac{\partial H_2}{\partial x_n} \\\
\vdots & \vdots & \ddots & \vdots\\\
\frac{\partial H_m}{\partial x_1} & \frac{\partial H_m}{\partial x_2}& \cdots & \frac{\partial H_m}{\partial x_n}
\end{pmatrix}
となる。数値計算では、勾配法を使い
\begin{align}
x \leftarrow x- \gamma \nabla J(x)
\end{align}
最適な初期値$x$を求める。
観測演算子が線形の場合、つまり、$\bf H(x)=Hx$を考える。$\bf \nabla J(x)=0$となる $x$ を $x^a$ として計算すれば、
\begin{align}
\bf \nabla J(x) &= \bf B^{-1} \left(x^a-x^b \right) + H^T R^{-1} \left(Hx^a-Hx^b-(y-Hx^b) \right) \\
&=0 \\
&\Rightarrow \bf \left(B^{-1} - H^T R^{-1}H \right) x^a = \left(B^{-1} - H^T R^{-1}H \right) x^b + H^T R^{-1}(y-Hx^b) \\
&\Rightarrow \bf x^a = x^b + \left(B^{-1} - H^T R^{-1}H \right)^{-1} H^T R^{-1}(y-Hx^b)
\end{align}
が求まる。この式は、最適内挿法から得られた式と同じである。
最後に
データ同化の基本的な事項だと思われる、最適内挿法と3次元変分法について説明した。3次元変分法において、観測演算子が線形ならば、最適内挿法に対応している。3次元変分法では、観測演算子が非線形な形でも計算することができるが、最適内挿法は線形でなければならない。最適内挿法の非線形への拡張は、拡張カルマンフィルターやアンサンブルカルマンフィルターとなるだろう。
実際に使う場面は、時間発展する偏微分方程式に適用される。機会があれば、3次元変分法を時間発展の系でも使える4次元変分法について紹介したい。3次元変分法と同様に、誤差関数を計算し勾配を求める。勾配の求める方法をアジョイント方といい、深層学習で用いる誤差逆伝播法と似た計算をする。
参考文献
記事は、以下の書籍を参考にしています。
データ同化 観測・実験とモデルを融合するイノベーション
淡路敏之・蒲地政文・池田元美・石川洋一 編著