はじめに
明治大学総合数理学部 Advent Calender 2017の18日目の記事です。今回は偏微分方程式と深層学習を繋ぐ話をしていきます。
目標
偏微分方程式の新しいSimulation方法であるDeep Galerkin Methodについて説明する。
Galerkin法について
そもそもGalerkin Methodとは有限要素法(FEM : Finite Element Method)のアルゴリズムの一つである。そうなると有限要素法とはなんだとなるので、以下にざっくりではあるが説明していく。
有限要素法とは偏微分方程式をSimulationする手法の一つである。その前に他のSimulation方法を紹介すると、数値解析の講義などでまず先に習う微分方程式の手法は常微分ならEuler法(Wikipedia)であったり、偏微分方程式なら陽解法や陰解法(Wikipedia)などといった差分法があげられる。これら手法は時間を等間隔に分割して微分方程式を近似的に解く手法である。
しかし、これら差分法とは違い、有限要素法ではまず空間を任意の大きさの三角形のメッシュで分割する。そして、その三角形を「うまく」貼り合わせることで、与えられた偏微分方程式を近似していくのが有限要素法である。この貼り合わせ方を与えるのがGalerkin Methodである。
有限要素法により得られる近似はある意味で元の偏微分方程式の最良の近似となっていることがわかっている。少し数学の言葉を使う。有限要素法で得た近似解は各メッシュ上で区分線形な関数を線形結合することによって得るものであるが、これら区分線形関数たちの張る空間上では、この近似解が元の解に「一番」近い解となっている。
(もっと立ち入りたい場合は関数解析と有限要素法を学んでください。キーワードとしてはSobolev空間、射影定理などなどです。)
Deep Galerkin Methodについて
Deep Galerkin Methodとは、深層学習を用いた偏微分方程式のSimulation手法である。そのアイディアは至極直感的で、Deep Neural Networkの関数の表現力を使って、偏微分方程式の解を近似しようというものである。具体的に紹介する。
まずは次のような偏微分方程式を考える。
\begin{eqnarray}
\begin{cases}
\dfrac{\partial u}{\partial t} + \mathcal{L} u(x, t) = 0, & \quad (x, t) \in \Omega \times [0, T] \\
u(x, 0) = u_0(x), & \quad x \in \Omega\\
u(x, t) = g(x, t) & \quad x \in \partial \Omega
\end{cases}
\end{eqnarray}
ここで、$\mathcal{L}$は楕円形微分作用素(Laplacianとかを考えてもらえれば良い)、$\Omega$は$\mathbb{R}^n$の有界領域であり、$\partial \Omega$は$\Omega$の境界である。
この偏微分方程式を近似するように、Deep Neural Networkを表す$f(x, t; \theta)$を最適化していこう。ここで使うDeep Neural Networkはinput層が$x \in \mathbb{R}^n$の$n$個分、$t$の1個とバイアス分の合わせて$n + 2$個分、hidden層が$h$の任意個、outputが$f(x, t; \theta)$の1個分である。
最適化にはまず評価関数を用意する必要があるが、次のように与える。$$J(\theta) = \left|\left|\frac{\partial f}{\partial t} + \mathcal{L}f\right|\right|^2_{\Omega \times [0, T], \gamma_1}+ \left|\left|f(x, t; \theta) - g(x)\right|\right|^2_{\partial \Omega \times [0, T], \gamma_2} + \left|\left|f(x, 0, \theta) - u_0(x)\right|\right|^2_{\Omega, \gamma_3}$$ここで、$\gamma_i$は何らかの確率密度関数であり、$||\cdot||$は$\gamma_i$分の重み付きの$L_2$ノルムすなわち
||\cdot||^2_{A, \gamma_i} = \int_A |f(y)|^2 \gamma_i(y)dy
により与えられる。つまり、元の偏微分方程式の各式と($L_2$の意味で)なるべく近くなるようなDeep Neural Networkの表現$f$を$J(\theta)$を最小化($f$真の解なら0となる)するような$\theta$を見つけることで実現しようという話である。
しかし、このままでは2つほど問題がある。一つ目は実際に$||\cdot||$を求めるには積分が必要であるという点。これは実装時に大変な問題になる。もう一つはどうやったら$\frac{\partial f}{\partial t}$と$\mathcal{L}f$を求められるのかという点。こちらはDeep Neural Networkの活性化関数として微分可能な関数、例えばSigmoid関数を使用していれば実際に求められる。
1つ目の問題
積分をするのは現実的はない。ならば$J(\theta)$を推定してしまえばよいのである。その手段は以下のとおりである。最適化のアルゴリズムも兼ねて説明する。
- $(x_n, t_n) \in \Omega \times [0, T], \ (\xi_n, \tau_n) \in \partial \Omega, w_n \in \Omega$をそれぞれ$\gamma_1, \gamma_2, \gamma_3$に従うように取ってくる。
- $s_n = [(x_n, t_n), (\xi_n, \tau_n), w_n]$として、$$J(\theta_n, s_n) = \left(\frac{\partial f}{\partial t}(x_n, t_n; \theta_n) + \mathcal{L}f(x_n, y_n; \theta_n)\right)^2 + (f(\xi_n, \tau_n; \theta_n) - g(\xi_n))^2 + (f(w_n, 0; \theta_n) - u_0(w_n))^2$$を求める。
- $\theta_n$と$J(\theta, s)$を何らかの最適化手法を用いて最小化する。
2で定義した$J(\theta_n, s_n)$では関数同士を比べる$L_2$ノルムではなく一つの点$s_n$における差を元にした評価関数となっている。$J(\theta_n, s_n)$は適当に定義したのではなく、$J(\theta)$の不偏推定量になっていることが、$s_n$のサンプリング手法からわかる。従って、積分をすることなく、最適化を行うことができる。
また、サンプリング手法を見てもらうと、ランダムに点を取っていることが分かるので、このDeep Galerkin Methodはメッシュに依存しないメッシュフリーな手法であることが分かる。
ところで、3番で何らかの最適化手法と述べたが弊Advent Calender 17日目の『深層学習の最適化アルゴリズム集』に詳しく載っている。
2つ目の課題
$\frac{\partial f}{\partial t}$と$\mathcal{L}f$をどうやって求めるのかについてであるが、これは単純に計算すればよいだけである。もちろんコンピュータに任せて自動微分などをしてもらえば求められるが、手計算でも求められるので実際に求めると良い。一応、$f$の一階微分及び、$f$のLaplacianは以下の設定のもとで求めておいたので、答え合わせに使ってほしい。
- 活性化関数 : Sigmoid
- input層 :
- $x_0$ : バイアス
- $x_1 \cdots, x_n$ : 空間変数
- $x_{n + 1} = t$ : 時間変数
- hidden層 : 計h層
- $m_1, \cdots, m_h$ : 各層のノード数(バイアス以外)
- $y_{0}^{[1]}, \cdots y_{0}^{[h]}$ : バイアス
- $y_{1}^{[k]}, \cdots, y_{m_k}^{[k]}$ : $k$層目のノード
- output層
- $b = y_{0}^{[h + 1]}$ : 出力活性
- $f = \sigma(b)$ : 出力($\sigma$は適当な$C^2$級関数)
- 重み
- $W^{[k]} = \left(w_{i_{k - 1}, {i_k}}^{[k]}\right) $ : ($k - 1$)層目から$k$層目への重み
$z^{[k]} = w_{i_{k - 1}, {i_k}}^{[k]} \cdot y_{i_k}^{[k]} \cdot (1 - y_{i_k}^{[k]})$とすると、
\frac{\partial f}{\partial x_i} = \frac{\partial \sigma}{\partial b}\left(\sum_{i_1, \cdots, i_{h + 1} = 0}^{m_1, \cdots, m_h, 0}w_{i, i_1}^{[1]}\left(\prod_{k = 1}^{h}z_k\right)\right)
\Delta f = \sum_{i = 1}^n \left[\left\{\frac{\partial^2 \sigma}{\partial b^2}\left(\sum_{i_1, \cdots, i_{h + 1} = 0}^{m_1, \cdots, m_h, 0}w_{i, i_1}^{[1]}\left(\prod_{k = 1}^{h}z_k\right)\right)^2\right\} + \left\{\frac{\partial \sigma}{\partial b}\left(\sum_{i_1, \cdots, i_{h + 1} = 0}^{m_1, \cdots, m_h, 0}w_{i, i_1}^{[1]}\left[\sum_{p = 1}^h \left\{\left[z_{p + 1}\left(\sum_{i_1, \cdots, i_p = 0}^{m_1, \cdots, m_p}w_{i, i_1}^{[1]}\left(\prod_{l = 2}^p z_l\right)\right)(1 - 2y_{i_p}^{[p]})\right]\left[\prod_{\substack{k = 1 \\ k \ne p}}^h z_{k + 1}\right]\right\}\right]\right)\right\}\right]
というとんでもない式が得られる。一般の楕円形作用素についてやろうものなら、さらにひっちゃかめっちゃかになる。
問題点
上で厳密に$\frac{\partial f}{\partial t}$と$\Delta f$の厳密な表示を与えたが、これをコンピュータに計算させようものなら大変時間がかかってしまう。しかも、これらが評価関数に入っている以上、毎ループこれらの値を計算せねばならない。これら式の計算時間に大きく関与しているのはhidden層の数であるので、上手い具合のパラメータを見つけることが必要とされる。
結論
Deep Galerkin Methodを紹介した。