『日経Robotics デジタル版（電子版）』のサービス開始を記念して、特別に誰でも閲覧できるようにしています。
有限要素法は2次元または3次元空間中の微分方程式の近似解を求める代表的な手法であり構造力学、流体、熱伝導、電磁場解析など広く使われている。例えば、領域$\Omega \subset \mathbb{R}^d$上で定義される時刻$t$、位置$\mathbf{x}$を入力とする関数$u(t,\mathbf{x})$が偏微分方程式$\partial_t u = F(t, \mathbf{x}, u, \partial_{\mathbf{x}}u,\partial^2_{\mathbf{x}}u, \ldots)$を満たす際の解を求めることができる。
Technical University of Munich教授のStephan Günnemann氏のグループが提案したFinite Element Network（FEN）1）は有限要素法の考えをグラフニューラルネットワーク（NN）に導入して、少数の観測データから、未知のダイナミクス（上記の$F$）を推定し、将来の値を予測する。高い予測性能を達成しつつ、人間にも解釈可能な結果を与え、既存の知識（例えばダイナミクスは対流項を含む等）を導入することができる。今回はこのFENについて紹介する。
有限要素法の基礎
有限要素法は空間を有限数のメッシュに分割し、偏微分方程式を解く。メッシュの頂点集合を$\mathcal{X} = \{\mathbf{x}^{(i)}\in \mathbb{R}^d\}_{i=1}^N$、各メッシュセル（以下セル）を$\mathcal{T} = \{\Delta^{(j)} | |\Delta^{(j)}|=d+1\}_{j=1}^{N_{\mathcal{T}}}$とする。例えば2次元であれば三角形に分割される。推定対象の関数$u$はスカラー場とし、ベクトル場の場合は複数のスカラー場を予測する問題とする。
有限要素法では関数$u$を$N$個の基底からなる線形結合$\tilde{\mathcal{U}} = \text{span} \{\phi^{(1)}, \ldots, \phi^{(N)} \}$で近似する。今回は基底にP1基底とよばれる区分線形関数（hat関数）を使うことにする。これは$j$番目の基底関数$\phi^{(j)}$は入力が$\mathbf{x}^{(j)}$の時$1$、$\mathbf{x}^{(j)}$と隣接する頂点との間で$0$に向かって線形に小さくなっていき、その他の領域では$0$となっているような関数である。この基底関数を使って$u$は次のように近似される。
\[u (\mathbf{x}) = \sum_{j=1}^N c_j \phi^{(j)}(\mathbf{x})\]
基底関数は時間不変であり係数$c_j$のみが時間に依存することに注意されたい。
Galerkin法はダイナミクス$F$を最もうまく近似する$u \in \tilde{\mathcal{U}}$を求めるため、近似誤差$R(u) = \partial_t u - F(t,\mathbf{x}, u, ...,)$と許容解空間$\tilde{\mathcal{U}}$が直交するという制約を利用する（真の解から許容解空間に垂線を下ろす）。この制約は内積を$\langle u, v \rangle_{\Omega} = \int_{\Omega} u(\mathbf{x}) v(\mathbf{x}) d\mathbf{x}$と定義した時、次のように表される。
\[\langle R(u), v \rangle_{\Omega}= 0 \quad \forall v \in \tilde{\mathcal{U}} \Leftrightarrow \langle R(u), \phi^{(i)} \rangle_{\Omega} = 0 \quad \forall i = 1, \ldots, N\]
この制約条件を変形すると次の線形方程式が得られる。
\[\mathbf{A} \partial_t \mathbf{c} = \mathbf{m}\]
ここで、$A_{i,j} = \langle \phi^{(i)}, \phi^{(j)}\rangle_{\Omega}$は質量行列とよばれ、$\mathbf{c}$は基底係数ベクトルであり、$m_i= \langle F(t, \mathbf{x}, u, \ldots),\phi^{(i)}\rangle_{\Omega}$は$F$の基底毎の影響をまとめたものである。この方程式を解くことで$\partial_t \mathbf{c}$、そこから$\partial_t u$が求められ、ODE（常微分方程式）を解くことで$u$が求まる（Method of Linesとよばれる）。
観測から場のダイナミクスを推定し将来を予測する
FENはこの有限要素法の解き方を参考にしてダイナミクス$F$を推定し、ODEソルバを使って将来時刻の特徴量を予測する。$d$次元中の$N$個の点$\mathcal{X} = \{\mathbf{x}^{(i)} \in \mathbb{R}^d\}_{i=1}^N$が与えられ、それぞれの点で$m$個の特徴量$\mathbf{y}^{(t_0, i)}\in \mathbb{R}^m$を測定したとする。
はじめに空間をドロネー三角形分割を使ってメッシュに分割する。また、先ほどと同様にP1基底を使って関数を近似する。
\[u_0(\mathbf{x})_k = \sum_{i=1}^N y_k^{(t_0)} \phi^{(i)}(\mathbf{x})\]
この時、P1基底関数を使っているため（基底は対象とする頂点上のみで$1$をとる）、頂点に対応する観測値と係数が一致する。
\[u_0(\mathbf{x}^{(i)}) = \mathbf{y}^{(t_0, i)} \]
そして、有限要素法の議論と同様にして次の線形方程式が得られる。$\mathbf{A}$の意味は同じであり、$\mathbf{M}$は$m$個の特徴量があることに対応し、特徴毎の基底の影響を列毎に並べている（後の$M_{i,j}$の式を参照）。
\[\mathbf{A} \partial_t \mathbf{Y}^{(t)} = \mathbf{M}\]
この式は$\mathbf{A}$の逆行列が解ければ$\partial_t \mathbf{Y}^{(t)} =\mathbf{A}^{-1} \mathbf{M}$と求まるが、計算量が大きい。そこで代わりに$\tilde{\mathbf{A}}_{ii} = \sum_{j=1}^N \langle \phi^{(i)},\phi^{(j)} \rangle _{\Omega}$を対角成分に持つ対角行列で$\mathbf{A}$を近似し、この逆行列を適用する。次に$\mathbf{M}$は頂点$i$に隣接するセル集合$\mathcal{T}_i$のセル（$CH(\Delta)$は$\Delta$の凸包）毎の影響に分解できる。
\[M_{ik} = \langle F(t, \mathbf{x}, u, \ldots)_k, \phi^{(i)}\rangle_{\Omega} = \sum_{\Delta \in \mathcal{T}_i} \langle F(t,\mathbf{x}, u, \ldots)_k, \phi^{(i)} \rangle_{CH(\Delta)} \]
さらにこの右辺は次のように変形できる。
\[\langle F(t, \mathbf{x}, u, \ldots)_k, \phi^{(i)}\rangle_{CH(\Delta)} = F^{(i)}_{\Delta, k} \cdot \langle 1, \phi^{(i)}\rangle_{CH(\Delta)} \]
このスカラー値$F^{(i)}_{\Delta, k}$はメッシュ毎に定数であり、ニューラルネットワークを使って予測する。
\[f_{\Delta, \theta}^{(t, i)} := f(t, \mu_{\Delta},\mathbf{x}_{\Delta}, \mathbf{y}_{\Delta}^{(t)}; \theta)^{(i)} \approx F^{(i)}_{\Delta}\]
ただし、$\mu_{\Delta}$はセルの中心、$\mathbf{x}_{\Delta}$はセルの頂点の中心からの相対座標集合、$\mathbf{y}_{\Delta}^{(t)}$はセルの頂点の特徴集合である。この$f$はセル毎のダイナミクスを表す。また$\langle1, \phi^{(i)} \rangle_{CH(\Delta)}$は、メッシュ分割に対し決まるので、一度計算しておき、キャッシュしておけばよい。
グラフNNでダイナミクスを予測する
この$\partial_t \mathbf{Y}^{(t)}$を予測する過程は次のようなグラフNNで表すことができる。このグラフは頂点集合$\mathcal{X}$とハイパー枝$\mathcal{T}$から構成される。グラフNNは頂点$v$毎に特徴ベクトル$\mathbf{h}_v$を持ち、ステップ毎に頂点を次のように更新する。
\[\mathbf{h}_v^{(p)} = f_{upd} ( \mathbf{h}_{v}^{(p-1)}, \sum f_{msg}(\Delta)_{\mathbf{x}^{(i)}} ) \] \[f_{msg}(\Delta)_{\mathbf{x}^{(i)}} = f_{\Delta}^{(t, i)} \cdot\langle 1, \phi^{(i)} \rangle_{CH(\Delta)} \] \[f_{upd}(\mathbf{h}^{(i)}, \mathbf{m}^{(i)} ) = \tilde{A}_{ii}^{-1}\mathbf{m}^{(i)} \]
ただし、$\mathbf{m}^{(i)} = \sum_{\Delta \in \mathcal{T}} f_{msg}(\Delta)_{\mathbf{x}^{(i)}}$は頂点$\mathbf{x}^{(i)}$で集約された特徴ベクトルである。
このグラフNNの1ステップが有限要素法による$\partial_t \mathbf{Y}^{(t)} = \mathbf{A}^{-1}\mathbf{M}$を求める操作に対応し、これを$P$ステップ適用した場合は$P$個の$\partial_t\mathbf{Y}^{(t)}$の各時刻の微分が得られる。このダイナミクス$f$の関数形で帰納バイアスを導入できる。例えば$f$が時刻や位置に依存しない場合は、それらを入力として使わないことで表現できる。この研究では対流があることを考慮した関数形を設計し、実験を行った。学習は各観測点の各時刻の測定値とFENが出力が一致するように学習する。なお、ODEソルバ上も誤差逆伝搬できるため、全体をEnd-to-Endで学習できる。
実験
実験では実際の観測で得られた海水の表面温度と水流、実際の煙流、シミュレーター上でシリンダー周辺の流体データを使って学習し、将来を予測できるかを評価した。FENは他の最高精度を達成したグラフNNと同程度の精度を達成し、特に観測点以外の点での予測などにおいては従来手法を超える性能を達成した。また、ダイナミクスに対流を考慮したモデルはさらに予測精度を改善できたとともにダイナミクスの表現として対流とそれ以外の効果などを明示的に分解することができた。
まとめ
FENのような事前知識をモデルに組み込む試みは物理シミュレーションを中心に広がっている。こうした事前知識を組み込むことで汎化性能を改善できるだけでなく、人間にとっても解釈しやすいモデルを作ることができる。今回はダイナミクスの推定問題であったが、今後の発展としてシミュレーションの高速化、高精度化も考えられる。
Preferred Networks 代表取締役 最高研究責任者
2006年にPreferred Infrastructureを共同創業。2010年、東京大学大学院博士課程修了。博士（情報理工学）。未踏ソフト創造事業スーパークリエータ認定。東京大学総長賞。
