📂機械学習

論文レビュー: 物理情報基盤ニューラルネットワーク(PINN)

概要

• レファレンスと数式の番号や表記法は、論文をそのまま踏襲する。

Physics-informed neural networks (PINN^{[ピン]と読む})は、数値的に解くために設計された微分方程式の人工ニューラルネットワークであり、2018年Journal of Computational Physicsに発表された論文Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equationsで紹介されました。この論文の著者は、応用数学、機械工学のM. Raissi, P. Perdikaris, G.E. Karniadakisです。

この論文で述べられている物理情報^{physics information}は、壮大に見えるかもしれませんが、実際には与えられた微分方程式自体を意味すると考えても良いでしょう。つまり、‘微分方程式を人工ニューラルネットワークで解く際、与えられた微分方程式を利用します’と言っているのと同じです。機械学習の論文を読むときは、このように見栄えの良い名前に惑わされないよう注意が必要です。

微分方程式の数値的解法においてPINNが注目される理由は、損失関数に関するアイデアがシンプルで理解しやすく、実装も簡単だからでしょう。実際に論文の例では非常にシンプルなDNNが紹介されています。

一般的に言われるPINNはSection 3.1で紹介されるモデルを指します。

実装

• PyTorchパイトチ

0. 抄録

著者はPINNを’与えられた非線形偏微分方程式を満たしながら、教師あり学習問題を解くために訓練された人工ニューラルネットワーク’と紹介しています。この論文で主に扱う2つの問題は、‘data-driven solution and data-driven discovery of partial differential equations’です。性能評価のために、流体力学、量子力学、拡散方程式などの問題を解いてみました。

1. 序論

最近の機械学習とデータ分析の進歩は、画像認識^{image recognition}、認知科学^{cognitive science}、ゲノム学^genomicsなどの科学分野で革新的な結果をもたらしていますが、複雑な物理的、生物学的、工学的システムに対しては（データ収集コストが高いため）少ない情報で望ましい結果を導き出す必要がある困難があります。このような*小さなデータ領域^{small data regime}*では、DNN、CNN、RNNなどの先進技術の収束性が保証されていません。

[4-6]で、データ効率が良く（=少ないデータで）、物理情報を学習できる（=微分方程式を解くことができる）方法についての研

究が進んでいます。非線形問題への拡張は、この論文の著者であるRaissiの後続研究[8,9]で提案されました。

2. 問題設定

人工ニューラルネットワークで表される関数は、入力値（偏微分方程式でのソリューション$u$の座標$x, t$を指す）とパラメータによって関数値が決定されるが、これら2種類の変数に対して微分を行うために自動微分^{automatic differentiation}を活用する。

このようなニューラルネットワークは、観測されたデータを支配する物理法則に起因する任意の対称性、不変性、または保存原理を尊重するように制約されている。これは、一般的な時間依存かつ非線形の偏微分方程式によってモデル化される。

この論文でこの文章が難しいと感じられるかもしれないが、私の考えでは簡単に言えば提案された人工ニューラルネットワークであるPINNが、与えられた微分方程式を満たす必要があるということだ。後述するが、微分方程式を満たす必要があるという条件を損失関数として使用するためである。

この論文の目的は、数理物理学におけるディープラーニングを進化させる新しいパラダイムのモデリングと計算パラダイムを提示することである。そのために、前述したように、この論文では主に2つの問題を扱う。一つは偏微分方程式のデータ駆動ソリューション^{data-driven solution}であり、もう一つは偏微分方程式のデータ駆動発見^{data-driven discovery}である。使用された全てのコードとデータセットはhttps://github.com/maziarraissi/PINNsで確認できる。この論文では、$L1$、$L2$、ドロップアウトなどの正則化なしに、ハイパーボリックタンジェントを活性化関数として用いたシンプルなMLPが使用されている。各例では、ニューラルネットワークの構造、オプティマイザー、学習率などが具体的に紹介される。

この論文では、以下のようなパラメータ化された非線形偏微分方程式の一般的な形^{parameterized and nonlinear partial differential equations of the general form}を扱う。

$$\begin{equation} u_{t} + \mathcal{N}[u; \lambda] = 0,\quad x \in \Omega,\quad t \in [0,T] \end{equation}$$

ここで、$u=u(t,x)$は(1)を満たす隠れた^{(=与えられていない=知られていない)}関数、つまり(1)のソリューションであり、$\mathcal{N}[\cdot; \lambda]$は$\lambda$でパラメータ化された非線形演算子（NonlinearのNに由来する）であり、$\Omega \subset \mathbb{R}^{D}$である。多くの数理物理学の問題^{problems in mathematical physics}は上記のような形で表される。例えば、1次

元粘性バーガース方程式を見てみよう。

$$u_{t} + uu_{x} = \nu u_{xx}$$

これは(1)で$\mathcal{N}[u; \lambda] = \lambda_{1} uu_{x} - \lambda_{2}u_{xx}$、$\lambda = (\lambda_{1}, \lambda_{2})$の場合である。与えられた方程式(1)に対して、扱うべき2つの問題はそれぞれ以下の通りである。

• 偏微分方程式のデータ駆動ソリューション: 固定された$\lambda$に対して、システムのソリューション$u(t,x)$は何か？
• 偏微分方程式のデータ駆動発見: 観測されたデータを最もよく表現するパラメータ$\lambda$は何か？

3. 偏微分方程式のデータ駆動ソリューション

セクション3では、以下の形式の偏微分方程式からデータに基づいたソリューションを見つける問題について扱う。

$$\begin{equation} u_{t} + \mathcal{N}[u] = 0,\quad x \in \Omega,\quad t \in [0,T] \end{equation}$$

つまり、$(1)$でパラメータ $\lambda$ が固定されている状況である。セクション3.1とセクション3.2ではそれぞれ連続時間モデルと離散時間モデルを扱う。方程式を見つける問題はセクション4で扱う。ここで言う’データ’の意味は以下で詳しく説明する。

3.1. 連続時間モデル

$(t,x) \in \mathbb{R} \times \mathbb{R}$ とすると、$u : \mathbb{R}^{2} \to \mathbb{R}$ である。これを人工ニューラルネットワークで近似するが、次のように実装されるシンプルなMLPを使用する。Juliaでは、

using Flux

u = Chain(
    Dense(2, 10, relu),
    Dense(10, 10, relu),
    Dense(10, 1)
    )

PyTorchでは、

import torch
import torch.nn as nn
import torch.nn.functional as F

layers = [2, 10, 10, 1]

class network(nn.Module):
    def __init__(self):
        super(network, self).__init__()
        layer_list = [nn.Linear(layers[i], layers[i+1]) for i in range(len(layers)-1)]
        self.linears = nn.ModuleList(layer_list)
        
    def forward(self, tx):
        u = tx

        for i in range(len(layers)-2):
            u = self.linears[i](u)
            u = F.relu(u)

        u = self.linears[-1](u)
        
        return u

u = network()

これで $u$ は入力ノードが $2$ つ、出力ノードが $1$ つの人工ニューラルネットワークとなる。$(2)$ の左辺を次のような関数 $f = f(t,x; u)$ として定義しよう。

$$\begin{equation} f := u_{t} + \mathcal{N}[u] \end{equation}$$

ここで $u$ は人工ニューラルネットワークであるため、$f$ も隠れ層のパラメータを持つ一種の人工ニューラルネットワークである。上記のような $f$ を物理情報に基づいたニューラルネットワーク^{physics-informed neural network, PINN}と呼ぶ。言い換えれば、与えられた偏微分方程式そのものである。$f$ に含まれる微分は自動微分で実装され、$u$ と同じパラメータを共有する。人工ニューラルネットワーク $u$ が $(2)$ のソリューションを適切に近似していれば、$f$ の関数値はどこでも $0$ であるべきだ。ここから、$f \to 0$ になるように人工ニューラルネットワークを学習させることが推測できる。

$(t_{u}^{i}, x_{u}^{i})$ を初期値、境界値が定義された領域の点とする。

$$(t_{u}^{i}, x_{u}^{i}) \in( \Omega \times \left\{ 0 \right\}) \cup (\partial \Omega \times [0, T])$$

$u_{\ast}$ を実際のソリューションとすると、初期条件と境界条件が与えられたということは、次のような値が与えられたということと同じである。

$$\left\{ t_{u}^{i}, x_{u}^{i}, u^{i} \right\}_{i=1}^{N_{u}},\quad u^{i} = u_{\ast}(t_{u}^{i}, x_{u}^{i})$$

理論上はこれらの値を無限に持つことになるが、数値的な問題では有限の点のみを扱えるので、$N_{u}$ 個を持っているとする。人工ニューラルネットワーク $u$ は $(t_{u}^{i}, x_{u}^{i})$ を入力として受け取り、$u^{i}$ を出力する必要があるので、これらがそれぞれ入力と対応するラベルとなる。

$$\text{input} = (t_{u}^{i}, x_{u}^{i}),\qquad \text{label} = u^{i}$$

これがPINNで学習する‘データ’である。それでは、損失関数を次のように設定できる。

$$MSE_{u} = \dfrac{1}{N_{u}} \sum\limits_{i=1}^{N_{u}} \left| u(t_{u}^{i},x_{u}^{i}) - u^{i} \right|^{2}$$

また、$f$は適切な点集合（理論的には解 $u_{\ast}$ が定義されるすべての点で満たされるべきだが、数値的には有限の点しか扱えない）$\left\{ t_{f}^{i}, x_{f}^{i} \right\}_{i=1}^{N_{f}}$で$(2)$を満たさなければならない。これらの適切な点を論文ではコロケーションポイント^{collocation points}と呼ぶ。コロケーションポイントに対して以下の損失関数を設定する。

$$MSE_{f} = \dfrac{1}{N_{f}}\sum\limits_{i=1}^{N_{f}} \left| f(t_{f}^{i}, x_{f}^{i}) \right|^{2}$$

つまり、$MSE_{f}$が$0$に近づくことは、物理的情報（偏微分方程式）を満たすことを意味する。したがって、人工ニューラルネットワーク $u$ を訓練するための最終的な損失関数は以下の通りである。

$$MSE = MSE_{u} + MSE_{f}$$

論文では、$MSE_{f}$を使用することで物理的情報を制約として設けることは、[15, 16]で初めて研究されたが、PINN論文ではこれを現代的な計算ツールで検討し、より困難なダイナミックシステムに適用したと説明されている。

物理情報に基づく機械学習^{physics-informed machine learning}という用語自体は、Wangの乱流モデリング^{turbulence modeling}に関する研究[17]で初めて使用されたとされる。しかし、PINN以前の研究では、サポートベクターマシン、ランダムフォレスト、FNNなどの機械学習アルゴリズムが単に使用されていたと説明されている。PINNがこれらと区別される点は、一般的に機械学習に使用されるパラメータに対する微分だけでなく、解の座標 $x, t$ に関する微分も考慮している点である。つまり、パラメータ $w$ を持つ人工ニューラルネットワークで近似された解を $u(t,x; w)$ とするとき、以前に提案された方法は偏微分 $u_{w}$ のみを利用したが、PINNは $u_{t}$ や $u_{x}$ などを利用して解を求める。このようなアプローチにより、少量のデータでも解をうまく見つけることができると説明されている。

この手続きがグローバル最小値に収束するという理論的な保証はないにもかかわらず、与えられた偏微分方程式が適切に定義されており、その解が一意であり、十分に表現力のあるニューラルネットワークアーキテクチャと十分な数のコロケーションポイント $N_{f}$ が与えられている場合、我々の方法は良好な

予測精度^{good prediction accuracy}を達成することが経験的に確認されていると論文には述べられている。

3.1.1. 例（シュレーディンガー方程式）

この例では、周期的な境界条件と複素数値を取る解に対して、提案された方法がうまく機能するかを重点的に確認する。例として、以下の初期条件と境界条件が与えられるシュレーディンガー方程式を扱う。

$$\begin{align*} ih_{t} + 0.5h_{xx} + \left| h \right|^{2}h &= 0,\quad x\in [-5, 5], t\in[0, \pi/2], \\ h(0,x) &= 2\operatorname{sech} (x), \\ h(t,-5) &= h(t,5), \\ h_{x}(t,-5) &= h_{x}(t,5) \end{align*}$$

問題の解 $h_{\ast}(t,x)$ は $h_{\ast} : [0, \pi/2] \times [-5, 5] \to \mathbb{C}$ として複素関数値を持つ。しかし、関数の出力が複素数になるように人工ニューラルネットワークを定義するのではなく、実部を担当する $u(t,x)$ と虚部を担当する $v(t,x)$ の2次元ベクトルが出力されるように定義する。簡単に言えば、入力と出力のノードがそれぞれ2つのMLPとして定義することである。

$$h(t,x) = \begin{bmatrix} u(t,x) \\[0.5em] v(t,x) \end{bmatrix}$$

この問題におけるPINN $f$ は以下の通りである。

$$f := ih_{t} + 0.5h_{xx} + \left| h \right|^{2} h$$

$h(t,x)$ と $f(t,x)$ のパラメータは、初期値に対する損失 $MSE_{0}$、境界値に対する損失 $MSE_{b}$、物理情報に対する損失 $MSE_{f}$ を最小化するように学習される。

$$MSE = MSE_{0} + MSE_{b} + MSE_{f}$$

$$\begin{align*} \text{where } MSE_{0} &= \dfrac{1}{N_{0}}\sum_{i=1}^{N_{0}} \left| h(0, x_{0}^{i}) - h_{0}^{i} \right|^{2} \quad (h_{0}^{i} = 2\operatorname{sech} (x_{0}^{i})) \\ MSE_{b} &= \dfrac{1}{N_{b}}\sum_{i=1}^{N_{b}} \left( \left| h(t_{b}^{i}, -5) - h(t_{b}^{i}, 5) \right|^{2} + \left| h_{x}(t_{b}^{i},-5) - h_{x}(t_{b}^{i},5) \right|^{2} \right) \\ MSE_{f} &= \dfrac{1}{N_{f}} \sum\limits_{i=1}^{N_{f}} \left| f(t_{f}^{i}, x_{f}^{i}) \right|^{2} \end{align*}$$

• 論文には $MSE_{b}$ の式にタイプミスがあるので注意すること。

ここで、$\left\{ x_{0}^{i}, h_{0}^{i} \right\}_{i=1}^{N_{0}}$ は初期値データ、$\left\{ t_{b}^{i} \right\}_{i=1}^{N_{b}}$ は境界でのコロケーションポイント、$\left\{ t_{f}^{i}, x_{f}^{i} \right\}_{i=1}^{N_{f}}$ は $f$ に対するコロケーションポイントである。

データセットを生成するために、従来のスペクトルメソッド^{spectral methods}を使用した。$h(0,x)$ での初期値データの数は $N_{0} = 50$、境界値データの数は $N_{b} = 50$ とし、ランダムに選んだ。また、$f$ のコロケーションポイントの数は $N_{f} = 20,000$ である。人工ニューラルネットワークは、100個のノードを持つ線形層を5層、層間の活性化関数としてハイパーボリックタンジェント $\tanh$ を使用して構築した。

Figure 1.

Figure 1では、上の図は予測された解 $\left| h(t, x) \right|$ のヒートマップを示している。下の図は、時間がそれぞれ $t = 0.59, 0.79, 0.98$ のときの予測された解と実際の解がどれだけ一致しているかを示している。相対的 $L_{2}$ ノルム^{relative $L_{2}$-norm} は $0.00197 = 1.97 \cdot 10^{-3}$ で、予測された解が正確な解と比較して約 $0.02$ の差異があることを意味している。したがって、PINNは少ない初期値データでシュレーディンガー方程式の非線形挙動を正確に捉えることができる。

現在扱っている連続時間モデルは、初期値が少なくてもうまく機能するが、コロケーションポイントの数 $N_{f}$ が十分に多くなければならないという潜在的な制約がある。これは空間の次元が2以下の場合はあまり問題にならないが、高次元の場合、必要なコロケーションポイントの数が指数関数的に増加する可能性があるため、問題になる可能性がある。そのため、次のセクションでは、多くのコロケーションポイントを必要としないようにするために、古典的なルンゲ・クッタ時間ステップスキームを活用した、より構造化されたニューラルネットワークを提案する。

3.2. 離散時間モデル

セクション3.1では、解を連続時間に対して近似した。この場合、人工ニューラルネットワークは全体の領域に対して同時に学習され、任意の点 $(x,t)$ に対して出力がある。このセクションでは、セクション3.1とは異なり、離散時間について扱う。つまり $t_{n}$ の値を知っているとき、$t_{n+1}$ の値を人工ニューラルネットワークで近似する方法について説明する。$(2)$ に $q$ ステージのルンゲ・クッタ法を適用すると、以下のようになる。

$$u(t_{n+1}, x) = u(t_{n}, x) - \Delta t \sum_{j=1}^{q} b_{j}\mathcal{N}\left[ u(t_{n}+c_{j} \Delta t, x) \right]$$

ここで $u^{n}(x) = u(t_{n}, x)$, $u^{n+c_{j}} = u(t_{n} + c_{j}\Delta t, x)$ と表記すると、

$$\begin{equation} \begin{aligned} u^{n+1} &= u^{n} - \Delta t \sum_{j=1}^{q} b_{j}\mathcal{N}\left[ u^{n+c_{j}}\right] \\ \text{where } u^{n+c_{j}} &= u^{n} - \Delta t \sum_{i=1}^{q} a_{j,i}\mathcal{N}\left[ u^{n+c_{i}}\right] \quad j=1,\dots,q \end{aligned}\tag{7} \end{equation}$$

上記の $q+1$ 個の式で、右辺の $\sum$ 項を全て左辺に移行する。そして、左辺を $u_{i}^{n}$ のように表記する。

$$\begin{equation} \begin{aligned} u_{q+1}^{n} &:= u^{n+1} + \Delta t \sum_{j=1}^{q} b_{j}\mathcal{N}\left[ u^{n+c_{j}}\right] = u^{n} \\ \\ u_{1}^{n} &:= u^{n+c_{1}} + \Delta t \sum_{i=1}^{q} a_{1,i}\mathcal{N}\left[ u^{n+c_{i}}\right] = u^{n} \\ u_{2}^{n} &:= u^{n+c_{2}} + \Delta t \sum_{i=1}^{q} a_{2,i}\mathcal{N}\left[ u^{n+c_{i}}\right] = u^{n} \\ &\vdots \\ u_{q}^{n} &:= u^{n+c_{q}} + \Delta t \sum_{i=1}^{q} a_{q,i}\mathcal{N}\left[ u^{n+c_{i}}\right] = u^{n} \end{aligned}\tag{9} \end{equation}$$

すると、これらの全ての値が $u^{n}$ に等しくなることがわかる。

$$u^{n} = u_{1}^{n} = u_{2}^{n} = \cdots = u_{q+1}^{n} \tag{8}$$

したがって、セクション3.2で言及されている物理情報とは、与えられた初期条件・境界条件と $(8)$ を指す。今度は $u(t_{n+1}, x)$ を求めるために二つの人工ニューラルネットワークを定義する。セクション3.1で使用した人工ニューラルネットワークは、正確な解 $u_{\ast}$ に収束することを期待する $u$ と、$u$ が満たすべき微分方程式 $f$ だったが、ここでは少し異なる。まず、人工ニューラルネットワーク $U$ を次の関数として定義する。

$$U : \mathbb{R} \to \mathbb{R}^{q+1}$$

つまり、入力層のノードが $1$ つ、出力層のノードが $q+1$ つのニューラルネットワークである。このニューラルネットワークの出力を以下の値とする。

$$U(x) = \begin{bmatrix} u^{n+c_{1}}(x) \\[0.5em] u^{n+c_{2}}(x) \\ \vdots \\[0.5em] u^{n+c_{q}}(x) \\[0.5em] u^{n+1}(x) \end{bmatrix} \tag{10}$$

このニューラルネットワークは、添付されたコード内の PhysicsInformedNN クラスで定義されている neural_net に該当する。

つまり、以下の学習プロセスで $U$ の出力の最後の成分が $u(t_{n+1}, x)$ に収束することを期待している。二番目のニューラルネットワークは、$U$ の出力と $(7)$ の定義を利用して、次のように定義される関数である。

3.2.1. 例（アレン・カーン方程式）

離散時間モデルにおける例として、以下の初期条件と周期的な境界条件が与えられるアレン・カーン方程式を取り上げる。

$$\begin{equation} \begin{aligned} &u_{t} - 0.0001u_{xx} + 5 u^{3} - 5u = 0,\qquad x\in [-1, 1], t\in[0, 1], \\ &u(0,x) = x^{2} \cos (\pi x), \\ &u(t,-1) = u(t,1), \\ &u_{x}(t,-1) = u_{x}(t,1) \end{aligned}\tag{12} \end{equation}$$

この例における $(9)$ に含まれる非線形演算子は以下の通りである。

$$\mathcal{N}[u^{n+c_{j}}] = -0.0001u_{xx}^{n+c_{j}} + 5(u^{n+c_{j}})^{3} - 5u^{n+c_{j}}$$

タイムステップ $t^{n}$ での $u$ の値を $u^{n,i}$ と表記する。

$$u^{n,i} = u^{n}(x^{n,i}) = u(t^{n}, x^{n,i}),\qquad i=1,\dots,N_{n}$$

問題は $u^{n}$ が与えられたときに $u^{n+1}$ を計算することであり、$\left\{ x^{n,i}, u^{n,i} \right\}_{i=1}^{N_{n}}$ は与えられたデータセットである。$(8)$ により、このデータセットに対して以下が成り立つ。

$$u^{n,i} = u_{1}^{n}(x^{n,i}) = \cdots = u_{q+1}^{n}(x^{n,i})$$

したがって、以下のような二乗誤差の合計^{sum of squared error (SSE)}を損失関数とする。

• ここではなぜ $MSE$ ではなく $SSE$ を使用するのかは明確ではない。連続時間モデルでは $MSE$ を使用していたが、離散時間モデルでは $SSE$ を使用しており、何らかの実験的な理由があると思われる。

$$SSE_{n} = \sum\limits_{j=1}^{q+1} \sum\limits_{i=1}^{N_{n}} \left| u_{j}^{n} (x^{n,i}) - u^{n,i} \right|^{2}$$

各 $u_{j}^{n}$ は $(9)$ によって計算され、この計算にはニューラルネットワーク $U$ の出力が使用される。このロスは、添付されたコード内の PhysicsInformedNN クラスで定義されている net_U0 に対応する。そして $U$ の出力は $(12)$ の境界条件を満たさなければならないため、以下のような損失関数を設定する。

$$\begin{align*} SSE_{b} &= \sum\limits_{i=1}^{q} \left| u^{n+c_{i}}(-1) - u^{n+c_{i}}(1) \right|^{2} + \left| u^{n+1}(-1) - u^{n+1}(1) \right|^{2} \\ &\quad+ \sum\limits_{i=1}^{q} \left| u_{x}^{n+c_{i}}(-1) - u_{x}^{n+c_{i}}(1) \right|^{2} + \left| u_{x}^{n+1}(-1) - u_{x}^{n+1}(1) \right|^{2} \ \end{align*}$$

これらの合計が最終的なロスである。

$$SSE = SSE_{n} + SSE_{b}$$

Figure 2.

Fig. 2では、上の図が正確な解のヒートマップを示している。下の図では、$t=0.1$ での $u$ を知っているときに $t=0.9$ での値を予測した結果を示している。下の左側の図では、青い線が正確な解であり、$\color{red}\mathsf{X}$ がデータとして使用された点を示している。下の右側の図では、青い線が正確な解であり、赤い線が予測された解である。

暗黙のルンゲ・クッタ法 (IRK) では $u^{n+c_{j}}$ を計算するためにすべての $j$ に対する連立方程式を解く必要があるため、$q$ が大きくなると計算コストが大幅に増加するが、この論文で提案されている方法では $q$ が大きくなってもそれに伴う追加コストは非常に少ないと説明されている。また、$q$ が小さい場合、IRK ではタイムステップ $\Delta t$ が大きいと正確な予測ができないが、PINN の場合は $\Delta t$ が大きくても正確に予測できると説明されている。

4. 偏微分方程式のデータ駆動発見

この章では、観測データがある場合に、偏微分方程式 $(1)$ のパラメータ $\lambda$ を見つける問題について扱う。詳細は以下の例を通じて説明する。

4.1. 連続時間モデル

$f$ を以下のように $(1)$ の左辺として定義しよう。

$$f = u_{t} + \mathcal{N}[u; \lambda]$$

セクション3の $(3)$ と異なる点は、$\lambda$ が固定された定数ではなく、学習が必要な未知のパラメータになったことである。

4.1.1. 例（ナヴィエ–ストークス方程式）

セクション4.1.1では、非圧縮性流体の実際のデータに関する例として、ナヴィエ–ストークス方程式によって表されるケースを紹介する。以下の2次元ナヴィエ–ストークス方程式を考える。

$$\begin{equation} \begin{aligned} u_{t} + \lambda_{1}(uu_{x} + vu_{y}) &= -p_{x} + \lambda_{2}(u_{xx} + u_{yy}) \\ v_{t} + \lambda_{1}(uv_{x} + vv_{y}) &= -p_{y} + \lambda_{2}(v_{xx} + v_{yy}) \end{aligned} \tag{15} \end{equation}$$

ここで、$u(t,x,y)$ は流体の速度ベクトルの $x$ 成分、$v(t,x,y)$ は $y$ 成分である。また、$p(t,x,y)$ は圧力、$\lambda = (\lambda_{1}, \lambda_{2})$ は未知のパラメータである。ナヴィエ–ストークス方程式の解は発散が $0$ となる条件を満たすため、以下が成立する。

$$\begin{equation} u_{x} + v_{y} = 0 \tag{17} \end{equation}$$

ある潜在関数 $\psi (t, x, y)$ に対して、以下のように仮定する。

$$u = \psi_{y},\quad v = -\psi_{x}$$

つまり、流体の速度ベクトルを $\begin{bmatrix} \psi_{y} & -\psi_{x}\end{bmatrix}$ とすると、$u_{x} + v_{y} = \psi_{yx} - \psi_{xy} = 0$ であるため、自然に $(17)$ を満たす。$u$ と $v$ を個別に求めるのではなく、$\psi$ を人工ニューラルネットワークで近似し、その偏微分によって $u, v$ を得る。実際の速度ベクトル場に対して、以下のように測定された情報があるとする。

$$\left\{ t^{i}, x^{i}, y^{i}, u^{i}, v^{i} \right\}_{i=1}^{N}$$

これに基づいて損失関数を以下のようにする。ここで、$u = \phi_{y}$, $v = -\psi_{x}$ であることを思い出す。

$$\dfrac{1}{N} \sum\limits_{i=1}^{N} \left( \left| u(t^{i}, x^{i}, y^{i}) - u^{i} \right|^{2} + \left| v(t^{i}, x^{i}, y^{i}) - v^{i} \right|^{2} \right)$$

そして、$(15)$ の右辺を左辺に定理し、それぞれを $f$ と $g$ として定義する。

$$\begin{equation} \begin{aligned} f &:= u_{t} + \lambda_{1}(uu_{x} + vu_{y}) + p_{x} - \lambda_{2}(u_{xx} + u_{yy}) \\ g &:= v_{t} + \lambda_{1}(uv_{x} + vv_{y}) + p_{y} - \lambda_{2}(v_{xx} + v_{yy}) \end{aligned}\tag{18} \end{equation}$$

すると $f, g$ の値は $\psi$ によって以下のように表される。（$p$ もニューラルネットワークで近似する）

$$\begin{align*} f &= \phi_{yt} + \lambda_{1}(\psi_{y} \psi_{yx} - \psi_{x}\psi_{yy}) + p_{x} -\lambda_{2}(\psi_{yxx} + \psi_{yyy}) \\ g &= -\phi_{xt} + \lambda_{1}(-\psi_{y} \psi_{xx} + \psi_{x}\psi_{xy}) + p_{y} + \lambda_{2}(\psi_{xxx} + \psi_{xyy}) \\ \end{align*}$$

損失関数に $f(t^{i}, x^{i}, y^{i}) = 0 = g(t^{i}, x^{i}, y^{i})$ という情報を加え、最終的に以下のようにする。

$$\begin{aligned} MSE &:= \dfrac{1}{N} \sum\limits_{i=1}^{N} \left( \left| u(t^{i}, x^{i}, y^{i}) - u^{i} \right|^{2} + \left| v(t^{i}, x^{i}, y^{i}) - v^{i} \right|^{2} \right) \\ &\qquad + \dfrac{1}{N} \sum\limits_{i=1}^{N} \left( \left| f(t^{i}, x^{i}, y^{i}) \right|^{2} + \left| g(t^{i}, x^{i}, y^{i}) \right|^{2} \right) \end{aligned} \tag{19}$$

次に、入力ノードが $3$ つ、出力ノードが $2$ つの人工ニューラルネットワークを定義する。この出力を $\begin{bmatrix} \psi (t, x, y) & p(t, x, y) \end{bmatrix}$ とする。すると、上記の損失関数を計算することができる。

ノイズがある場合とない場合のデータについて実験を行い、どちらの場合も高い精度で $\lambda_{1}, \lambda_{2}$ を予測できたことが示されている。また、圧力 $p$ に関するデータが与えられていない場合でも、人工ニューラルネットワークがパラメータと共に $p$ もかなり正確に近似できることが示された。具体的な実験設定、結果、参照解の求め方については、論文に詳しく記載されている。

5. 結論

この論文では、与えられたデータが満たす物理法則をエンコードする能力があり、偏微分方程式で説明できる新しい種類のニューラルネットワーク構造である物理情報に基づいたニューラルネットワークを紹介した。この結果から、物理モデルに対してディープラーニングが学習できることがわかった。これは多くの物理的シミュレーションに応用可能である。

しかし、著者は提案された方法が偏微分方程式を解くための既存の方法、例えば有限要素法^{finite element method}、スペクトル方法^{spectral methods}などを置き換えるものだと考えるべきではないと述べている。実際にセクション3.2.ではルンゲ・クッタ法をPINNに適用している。

PINNを実装するために、どれくらいの深さのニューラルネットワークが必要か、どれくらいのデータが必要かなどのハイパーパラメータに関する問題についても、著者は解決策を提案しようとした。しかし、ある方程式で効果的な設定が他の方程式ではそうではないことが観察されたと述べている。