Physics-Informedニューラルネットワーク wi2 ニューラル演算子
レビュー論文[*1]、[*2]を基に、PINNs†について整理した。PINNsとは、 Physics-informed Neural Networksの略語である。
† SEO対策的として、PINNからPINNsにした。
Ⅰ-2 Case1 PGNNを使った渋滞予測
Ⅰ-3 PINNsのベストプラクティス
Ⅰ-4 Case2 疫学×PINNs
Ⅰ-5 PINNsのベンチマーク
Ⅰ-6 PINNs学習において、高精度かつ高速なオプティマイザを提案する論文
Ⅰ-7 Case3 人工衛星の軌道推定
Ⅰ-8 コロケーション点の構成を扱った論文2編
Ⅰ-9 Case4 アルツハイマー病の早期診断
補足Ⅰ-1 PINNsは有限要素法に勝てない？
(0) 前捌き
PINNsの前半部分、Physics-informedの正式な和訳はないようであるが、｢物理学に基づいた｣、｢理論に基づいた｣、｢物理情報付き｣といった訳語が一般に用いられている。以下では訳語を用いることなく、表記をPINNsで統一する。
歴史的な話をすると、事前知識を学習プロセスに組み込みという概念は、少なくとも1994年には表出している。明確な形でPINNsが提示されたのは2017年であるが、その隆盛は、Raissi の論文[*3]がarXivに投稿された 2019年からと目されている。
議論の拡散を避けるため、多様で複雑な物理系を理解するために様々な努力が重ねられてきたこと、シミュレーションもその一環であることは、所与とする。
(1) 物理からPINNsを見る
PI＝物理＋NNと考えられるから、まず物理からPINNsを見よう。物理からニューラルネットワーク(NN)を見たとき、物理がNNを欲する場面は、例えば、次のような場面であろう：㈠物理法則が部分的に知られている(保存則は知られているが、構成方程式は知られていない)＋データが不十分という状況下で、物理系を理解したい。㈡確率的偏微分方程式で記述される物理系を理解したい。㈢厳密には、非整数階(fractional)微積分方程式・偏微分方程式で記述される｢乱流、粘弾性体、その他の異常輸送過程など、長距離の時空間相互作用を伴う問題｣を理解したい(ちなみに、非整数階常微分方程式系で疫学モデルを構築したケースは、こちら)。
機械学習･深層学習の文脈からは、㈠のケースが分かりやすい。複雑な物理系で、教師データ(有用な実験データ、測定データ)が十分に存在することは前提とできない(もちろん、シミュレーション･データを教師データとするという発想はありだが、それはサロゲート･モデリングのアプローチである)。
そこで、大量のデータとマシンパワー頼みで解を導出するというNNのセントラル･ドグマから脱出し、支配方程式を介して物理則を満たすという制約をつけて、NNを構築するという新アプローチが、PINNsだとも言える。物理則を満たすという制約は、支配方程式そのものを損失関数として用いることで、学習プロセスに導入される。
なお、敵対的生成ネットワーク(GAN)は、高次元の確率的偏微分方程式を解くための強力なアプローチである。
(2) ニューラルネットワークからPINNsを見る
次に、PINNsをNNから見よう。機械学習･深層学習寄りの言葉を使うと、PINNsの構図がより腹落ちし易いかもしれない。PINNsは、以下のようにも表現可能だろう：支配方程式の解に関する体系的に構造化された事前知識を｢ドメイン知識｣として用いることで、精度向上を目指したニューラルネットワーク。
方法論的(アルゴリズム的)には、PINNsは『支配方程式を解く問題を、損失関数(PINNsでは、右辺が0の支配方程式と同値)を最小化する最適化問題に変換している』と表現できる。ここまで来ると、イメージが湧いてくるだろう。もちろん、正則化等の議論は、PINNsでもあるが、ここでは、そこまで立ち入らない。
細部を理解するために、PINNs をモジュール分解する。PINNsは、1⃣近似、2⃣損失(関数)、および3⃣最適化の 3 つからなると考えることができる。1⃣～3⃣については、(4)を参照。
なお、Physics-informed機械学習という用語もあるが、その意味は自明であろう。また、PINNsの変形アプローチには、PCNN、CPINNs、(hp-)VPINNsなどがある。
(3) シミュレーションからPINNsを見る
物理の支配方程式は、非線形偏微分方程式であることが多い。放物型や双曲型ならまだしも、楕円型だと厄介である。ナヴィエ･ストークス方程式のように、一般解の存在が数学的に証明されていない方程式もある。それは、さておき、物理における支配方程式は、解析解が得られないケースがほとんどであり、解を導出するために数値計算が選択される。数値計算には、有限差分法/有限要素法/有限体積法といった決定論的手法と、モンテカルロ法のような確率論的手法がある。シミュレーションからPINNsを見るとは、PINNsと数値計算の優劣を比較することを意味している。
PINNsの優位性は、 ①逆問題が得意、②次元の呪いから解放、③堅牢性、である。④(部分的に)柔軟、⑤(部分的に)高速、を上げても良いかも知れない。⑥外挿が可能、も上げられる。
①実課題への適用を想定した場合の有用性という意味で、逆問題が得意という特徴は、特筆に値する。逆問題は、常に解けるとは限らない。逆問題が解ける条件(適切性条件)を満たしている場合、その逆問題はwell-posedと呼ばれる。一方、満たしていない場合はill-posedと呼ばれる。PINNsは、逆問題が適切性条件を満たさない場合でも、逆問題の解決に有効で効率的、という性質を持つ。また、逆問題を解くために必要な追加作業が少ない、という性質も有する。
GANを始めとする生成モデルをイメージすると、PINNsは逆問題が得意、の納得感が増すかもしれない。
②次元の呪いとは、問題の次元が増加すると、計算コストが指数関数的に増加することを言う。コルモゴロフ型偏微分方程式 (熱方程式またはブラック-ショールズ方程式)では、残差が｢問題の次元に依存しない値で抑えられる｣ことが証明されている[*4]。一方、例えば有限要素法では、問題の次元とともに指数関数的に増加する。
③ここで言う堅牢性とは、｢観測不可能な潜在的な非線形状態変数を記述する能力｣を指している。平たく言うと、ある変数の測定値が欠落していても、PINNsはその変数を正確に予測することができる。また、初期条件や境界条件が指定されていない場合や、(支配方程式の)パラメータの一部が未知である場合など、従来の数値計算手法が通用しないような場面でも対応可能である。
④グリッドベースの有限差分法では不要だが、有限要素法や有限体積法では不可欠なメッシュ生成が不要なため、不規則領域や移動領域の問題を容易に扱うことができる。
⑤モンテカルロ法に比べれば、PINNsは高速である。またPINNsでは、計算コストの高いメッシュ生成が不要という意味で高速、と言っても良いかもしれない。並列化が容易な点も高速性には寄与する。ただ、順問題でグリッドベースの有限差分法に比べて速いか、というと、そんなことはない。むしろ勾配降下法を用いているため、解くという箇所に限定すれば、"非常に遅い"。
⑥PINNsは、境界値問題において空間外挿を行うことができる[*5]。
一方、劣位性としては、以下が上げられる。
㊀NN自体の解明が十分ではない。そのため、例えば、なぜPINNsが次元の呪いから解放されるのか、理由がわかっていない(ディープなNNは、"滑らかな方向"を見つけることができるため、次元の呪いを回避できるようである)。さらに言えば、NNは、勾配に基づく最適化によって偏微分方程式の解を求めることができるのか？という根本的な疑問に対する答えは得られていない。
㊁PINNsの場合、損失関数が非凸関数であるため、解の一意性が保証されない。またPINNsは勾配ベースの最適化であるから、大域的な最小値が見つかるとは限らない(NNの幅を”十分広く”とれば、大域的最適解が勾配法で求まることは証明されている)。
㊂PINNsは固定メッシュやグリッドを必要としないため、トータルでは計算コストが低いとされる。有限要素法や有限体積法では、メッシュ生成の計算コストが高い(純粋に解く、という意味では、速くない)。しかしPINNsは、学習段階で境界条件を確立する必要があるため、境界条件が変化した場合は、新しいネットワークを作成する必要がある。
㊃PINNsに限らずNNに共通であるが、層が深くなるにつれて、勾配消失に対して脆弱になる。
(4)PINNsのお作法
1⃣ 近似
① アーキテクチャ
PINNsで用いられる主要なNNのアーキテクチャは、フィードフォワードNN(FFNN)、畳み込みNN(CNN)、再帰型NN(RNN)である。
② 方法論[*6]
損失関数L(θ)を最小化するNNのパラメータ･セットθを見つける。順問題の場合は、L(θ)＝wf×Lf(θ)＋wb×Lb(θ)。逆問題の場合は、L(θ)＝wf×Lf(θ)＋wb×Lb(θ)＋wi×Li(θ)、である。wf、wb、wiは、単なる重み†。Lf、Lb、Liについては、❷①を参照。
† 単なると述べた意図は、物理的な視点からは｢単なる｣、ということであった。機械学習･深層学習の視点からは、"単なる"は、語弊がある。この重みを適当に調整することで、｢複雑形状を有する問題｣や｢マルチスケール問題｣で精度が向上する(グランドトルゥースにより近づく)。☞Ⅰ-5を参照。
③ 正則化
あまり議論されていないものの、チホノフ正則化(≒リッジ回帰)が、一般的に使われているようである。逆問題との平仄を考慮しているのかもしれない。
④ 活性化関数
深層学習の学習プロセスで一般に使われる活性化関数は、ReLU、シグモイド、双曲線正接(tanh)である。一般的な深層学習の文脈では、Swish のような調整可能な活性化関数を学習することが推奨されているものの、PINNsにおいては、無限に微分可能なtanhを使用する傾向がある。
【参 考】
NVIDIAのサイト[*42]には、｢self-scalable tanh(Stan)活性化関数は、PINNsにおいて、より優れた収束特性を示し、精度を向上させることが知られている｣とある→22年9月にアップデートされた、NVIDIA Modulus v22.09では、Stan活性化関数がサポートされている。[*43]
論文[*44]によれば、(PINNsの活性化関数として採用する場合)tanhには、以下2つの問題がある。
㊀ PDE残差†計算に由来する勾配消失問題。この勾配消失は、活性化関数の飽和が原因である。そして、この問題は微分階数が高いほど、問題の(方程式の)次元が大きいほど、シビアになる。
†PDE(偏微分方程式Partial Differential Equation)残差とは『ネットワーク(PINNs)からの出力が、正解でないために支配方程式を満たさない故に、発生する誤差(残差)』を指している。
㊁ 重みの大きさを変化させる必要があるため、学習が不安定になる、という問題。PINNsの出力は、PDEと初期/境界条件によって大きさがわからないため、正規化することができない。このことと、tanhの出力が-1～1であることのミスマッチが学習を不安定にさせる。つまり、最終層で重みを、急激に大きくしなければならない(可能性がある)ため、学習が不安定になる(ことがある)。
2⃣ 損失関数
損失関数は、一般に、平均自乗誤差(MSE)の形式で表現される。
① 順問題
損失関数は、順問題では、㊀と㊁との2つに分解される・・・㊀Lf：支配方程式の不一致に関する損失関数。つまり、厳密解(uexact)では、右辺が0となる支配方程式において、uexact≠uであるが故に、0とはならない損失である。ここでuは、NNの出力＝近似解である。㊁Lb：初期条件及び境界条件の不一致に関する損失関数。つまり、uexact≠uであるが故に、初期条件及び境界条件の値がズレることで生じる損失である。
② 逆問題
逆問題の場合は、逆問題が満たす式が加わる・・・㊂Li：逆問題が満たす式の不一致に関する損失関数。uexact≠uであるが故に、0とはならない損失である。逆問題の場合は、Ldataの必要性が増す。
③ Ldataについて
順問題の損失関数は、LfとLbで構成される。ただ、PINNsにはNNの側面もある†1。したがって厳密には、学習データを与える時空間座標{ti,xi}におけるu(ti,xi)と、学習データの値uiが異なることによる損失(学習データによる損失Ldata)も考慮する必要がある。[*1]と[*2]では考慮されているが、[*6(再)]では考慮されていない(というより、現れないことの方が多い)。
一見、Ldataは、Lfでカバーできるようにも見える(すなわち、Ldataを考慮しなくても構わないように思える)。ただ、前述の通り、PINNsには解の一意性が保証されない。そのため必ずしも、Lfで、Ldataをカバーできるとは限らない(と理解している)。故に、dataを考慮する意味はある。そして、サンプリングが上手く出来ていないとか、コロケーション点が余りに少ないといった状況だと、Ldataを考慮しない影響が大きくなるだろう(と理解している)†2。
また、Ldataは正則化項の役割も果たす、と思われる(PINNsでは、複数の損失関数があり、それぞれが、正則化項のような役割を果たしているかもしれない)。
流体力学におけるPINNsのレビュー論文[*84]には、｢損失関数の定義は、問題に依存するため、問題の種類によって一部の項が消える場合がある。既知のパラメータと初期/境界条件を使用して、流体力学の順問題を解く場合、Ldataは、必ずしも必要ではない｣とある。
†1 NN＝近似器という意図。
†2 実際のところ、コロケーション点の追加や再サンプリング等を通じて、(相対L2)誤差を減らす取り組みが成される。つまり、いずれにせよ、誤差を減らす取り組みの中で、｢コロケーション点が少なすぎるorサンプリングが不十分｣という問題は解消される。従って、順問題においてLdataは、実質的に｢不要｣と考えて良いのだろう。
3⃣⃣ 最適化
ほとんどの PINNs 文献で損失関数は、勾配ベースの最適化スキームであるAdam(Adaptive Momentum)を使用したミニバッチ･サンプリングと、制限付きブロイデン･フレッチャー･ゴールドファーブ･シャンノ(Broyden-Fletcher-Goldfarb-Shanno)法(L-BFGS法)を使用している。
ただし、学習の開始時点からL-BFGSを使うと、上手く学習できないことがあるらしい。そのような場合、最初にAdamを使って粗い最適化を実施し、適当なタイミングでL-BFGSに切り替えるという2段階アプローチが採用される。さらにいうと、切り替えタイミングも重要である。つまり、｢Adam＋L-BFGS｣2段階アプローチだと試行錯誤が必要になる。計算コストは、L-BFGS＞Adamである(ので、Adamで事足りれば、幸い)。
PINNs用の新しい最適化手法｢エネルギー自然勾配法｣については、こちらを参照。
4⃣ 自動微分
自動微分は、PINNsにおいて地味に重要である。PINNsの損失関数は、支配方程式であり、そのほとんどは偏微分方程式であるから、微分計算が欠かせない。PINNsは、微分可能な関数を使ってNNを構築しているので、NN内で表現する近似解uの微分を厳密に実行できる。詳しくは例えば、[*53]を参照。
TensorFlowやPyTorchを含むほとんどの機械学習ライブラリは自動微分を提供する。例えば、TensorFlowにおいて∂u/∂t は、tf.gradients(u, t) として計算できる。2階微分は tf.gradients を2回適用して計算する。高階微分の計算をより効率的に行う機械学習ソフトウェアライブラリも存在する。
(5)実応用例
応用分野として、｢流体力学、分子動力学、量子化学、材料科学、地球物理学｣といった分野が頭出しされている。イメージを掴むため、具体例をあげる。
(なお、特出ししたトピックスは、以下の通り：渋滞予測、感染予測、衛星の軌道予測。)
1⃣ 磁場閉じ込め核融合プラズマ
入力データ＝プラズマの観測データ(電子密度及び温度)、物理モデル＝磁化プラズマの二流体方程式であるブラジンスキー方程式、出力＝プラズマ乱流場。→プラズマ診断やモデル検証を効率化できる。
2⃣ 血行動態
入力データ＝頭部MRI画像データ、物理モデル＝動脈中の血流動態(流量・圧力)を表現するための力学モデルであるウィンドケッセル･モデル、出力＝三次元頭蓋内動脈瘤の血行動態。→事前に学習されたモデルは、新しい患者の状態にすばやく変更できる。
3⃣ 地球物理学
①理研は、地震に伴う地殻変動のモデリングにPINNsを活用した(22年12月1日[*7])。ここでは、PINNsを｢物理法則を組み込んだ深層学習｣と和訳している[*8]。
入力データ＝地球内部の位置データ、物理モデル＝岩石内部に働く力と岩石の変形の関係を記述する弾性体力学、出力＝地殻の変位量。→複雑な形状を持つ場合、地下構造の性質が急激に変化する場合においても、正しく解析できることを実証した。PINNsの特性が活かされている。なお、地震時に生じる断層両側の動きのズレを表現するため、必要に応じて、入力(位置データ)を変換している。
②国立研究開発法人海洋研究開発機構は、ベイズ推定を援用するPINNs(ベイジアンPINNs)の枠組みを使って、地震波速度(～地下構造)の｢不確実性定量化｣を行った、と発表した[*35](23年10月11日)。なお、本件で、PINNsは｢物理情報に基づく深層学習｣と訳されている。
不確実性定量化(Uncertainty Quantification:UQ)[*36]は、やや広い概念であるが、数値計算の文脈では、変数に誤差(例えば、離散化誤差)が含まれていることを所与として、その誤差が、数値計算結果に与える影響を定量化することを言う。
より具体的に、数値計算の文脈で、不確実性定量化を述べると、次のようになる：誤差を何らかの確率分布Aに従う変数とすれば、数値計算結果も何らかの確率分布Bに従う変数考えることができる。Aからサンプリングされた誤差と、Bからサンプリングされた計算結果の間との間に、統計的な関係を探ることが、(数値計算の文脈における)｢不確実性定量化｣である。
さらに、ベイズ推定を使った数値計算の文脈で言うと、先の確率分布Bが、事後分布になる。つまり、ベイズ推定を使うと、自然に｢不確実性定量化｣が行える。なお、ベイジアンPINNsの文脈では、速度とか外力と言った、物理量が｢パラメータ｣になり、いわゆるθはハイパーパラメータである。
アーキテクチャ的には、ベイジアンPINNs[*37]は、㊀ベイジアン･ニューラルネットワークと、㊁PINNsとが、融合したニューラルネットネットワークである。㊀で、適当な｢パラメータの事前分布｣を設定し、㊁でハイパーパラメータ(～尤度)を最適化計算で推定する。その結果を、㊀に送り、｢パラメータの事後分布｣の推論を行う。本研究では、関数空間粒子ベース変分推論(スタイン変分勾配降下法、SVGD[*38])を使って、事後分布を変分推論[*39]している。
4⃣ デジタルツイン
①以下の理由から、PINNsは、デジタルツインの構築に適する：物理ベースの計算モデルは、計算コストが高い前処理と校正手順(メッシュ生成や初期条件及び境界条件の校正など)に大きく依存している。このため、リアルタイムの意思決定環境での使用が妨げられている。さらに、多くの複雑な自然界の物理モデルは、高々保存則として、部分的に知られているだけであるため、適切な構成則が仮定されていない限り、閉じた方程式系を提供することはない。
②大阪大学は、PINNsを使うことで、るつぼ内のシリコン融液の流速分布及び温度分布を、正確に再現できたと発表(応用物理学会秋季学術講演会、22年9月14日)[*9]。PINNsを現実の系に適用し、流動と温度におけるその物理的整合性を確認した世界で最初の研究事例、と主張している。なお、本研究ではPINNsを｢物理法則の支配方程式を学習したニューラルネットワーク｣と和訳している。
5⃣ その他
①産業用アプリケーション：潤滑油劣化プロセスや腐食疲労による亀裂の成長とベアリングの疲労を表現するために活用されている。また、常微分方程式に簡略化された動揺方程式(スイング方程式)を解くことによって、電力システム アプリケーションに対処するためにも活用されている。スイング方程式は、発電機の同期ダイナミクスを記述する方程式である。
②強化学習と深層学習を組み合わせた｢深層強化学習｣を利用するためのツールとして、PINNs を使用できる。
③DeepONet[Ⅳ【2】を参照]とPINNsによって符号化された物理を組み合わせることで、電気対流や極超音速などにおいて、外挿を用いたリアルタイム高精度予測を達成することが可能となっている。極超音速は2022年、防衛分野において突然パワーワードとして浮上した、あの極超音速である。
④不飽和浸透問題：不飽和浸透問題における逆問題(パラメータ解析)にPINNsを適用していた例[*78]。入力データ＝不飽和浸透挙動の計測データ、物理モデル＝Richards式、出力＝不飽和浸透特性パラメータ。オプティマイザーはL-BFGS、活性化関数はtanh。ユニット数20、隠れ層7。※細かく言うと、PINNsというよりPGNNか?
【参考】
別のレビュー論文[*18]によると、2019～2022年半ばにかけて公開されたPINNsに関する論文の国別割合(やや定義が曖昧であるが･･･)は、米国が29%、中国が25%、ドイツが10%。日本はその他に含まれていて詳細不明。
(6)ソフトウェア
1⃣ ソルバー
NVIDIA Modulus(旧SimNet)、DeepXDE(米ブラウン大で開発→米ペンシルバニア大でメンテ、アパッチ･ライセンス)、NeuroDiffEq(ハーバード大応用計算科学研究所で開発･メンテ、MITライセンス)、NeuralPDE(MITライセンス)、が主要なソルバーである。PyDEns(言語:Python、バックエンド:TensorFlow、アパッチ･ライセンス)は、2020年に開発停止したようである。Modulus 以外は、オープンソース。NeuralPDE以外は、言語が全てPython。NeuralPDEはJulia。バックエンドは、ModulusとDeepXDE がTensorFlow。NeuroDiffEqはPyTorchで、NeuralPDEはJulia。
①DeepXDEは非整数階偏微分方程式にも対応している。また、複雑な領域形状をサポートしており、 ユーザーコードをコンパクトに保ち、 数学的定式化に近い形にすることが可能。さらに、DeepXDEは、すべてのコンポーネントが疎結合であるため、構造化されており、高度に設定可能。
なお、[*6]のfirst authorであるLuがDeepXDEの開発者である。
②NeuroDiffEq に対しては『対応する損失を追加するのではなく、明示的な境界制約を強制するため、ライブラリがサポートしていない任意の境界には不十分であるように見える』という記述がある。
③NVIDIA Modulus は、空間領域構成法および平面充填操作による構造物定義を提供し、幅広いジオメトリをパラメーター化できる。損失関数は、通常行われる総和(いわゆる∑)ではなく、積分を使用している。さらに、学習率アニーリングを使用して、最適化スキームにおける損失関数の重みパラメーターを微調整する。さらに、他の多くのパッケージとは異なり、強形式及び弱形式の 偏微分方程式を処理できる。
2⃣ ラッパー
SciAnn(言語:Python、バックエンド:TensorFlow、MITライセンス)、ADCME(言語:Julia、バックエンド:TensorFlow、MITライセンス)、GPyTorch(言語:Python、バックエンド:PyTorch、MITライセンス)、Neural Tangents(言語:Python、バックエンド:JAX、アパッチ･ライセンス)が主要なラッパーである。
①SciANN は、バッチ最適化や転移学習のためのモデルの再利用など、Keras の機能の多くを継承している。 またニューラルネットワークの構築を抽象化するように設計されているため、複雑な関数形式を設定する柔軟性を提供する。
②ADCMEは、逆問題を解くための計算フレームワークを提供するライブラリである。TensorFlow の機能を拡張および強化することによって開発された。
③GPyTorchは、ガウス過程の計算を行うライブラリである。自動微分機能は装備されている。ベイジアンPINNsを実施する場合、不確実推定を行う場合には、有効かもしれない。
④Neural Tangentsは、Google AIが開発したライブラリです。無限幅のニューラルネットワークを高速かつ簡単に構築及び学習させることを目的としている。学習モデルとして、ガウス過程とニューラル･タンジェント･カーネル(NTK)(割と正式な日本訳は、神経接核であるが、これでは意味が全く分からない!)が実装されている。NTKは、NNをカーネル法の枠組みで理解するアイデアである(カーネル法は、第1層目を固定した横幅無限の2層ニューラルネットワークと理解できる)。
(7)課題と対応策
0⃣ 前振り
Ⅰ-3ベストプラクティスとⅠ-5ベンチマークにおいて、まとまった形で課題と対応策を議論している。Ⅰ-3では、(主に、課題を解消するべく考案された複合モデルを提示している[*34]とは異なり)｢PINNsの学習効率と精度を向上させるベストプラクティス｣を紹介している。ペンシルベニア大学の研究者による論文[*40](23年8月16日@arXiv)である。Ⅰ-5では、まず課題と対応策を、それぞれカテゴライズする。対応策では選択肢を示しつつ、最適な策を選んでいる。その上で、様々な物理系(支配方程式)に対応策を適用し、その有効性を議論している。
下記1⃣では、(PINNsに限定されないNN全体に及ぶような)もう少し緩く、散発的･五月雨式に課題を上げる。2⃣では(1⃣と重複するが)課題を上げ、対応策を軽くあてている。
1⃣ 五月雨式故に網羅的なPINNsの課題
まず、シミュレーションにおけるPINNsの劣位性(課題)を、簡潔に再掲する。
㊀ なぜ、PINNsが次元の呪いから解放されるのか、理由がわかっていない。
㊁ 損失関数が非凸関数であるため、解の一意性が保証されない。
㊂ 境界条件が変化した場合は、新しいネットワークを作成する必要がある。
㊃ 層が深くなるにつれて、勾配消失に対して脆弱になる。
加えて、以下の課題が指摘されている。
❶ 観測データを情報の点源と見なし、情報を外部に伝播する前にそれらの点で残差を最小化する必要がある逆問題については、難がある。
❷ マルチスケール問題では、高周波成分の学習に苦戦し、しばしば学習に失敗する。
❸ PINNsの学習における学習率の影響は、十分に調査されていない。
❹ 点毎に定義される損失関数を用いたNNは、高次元の問題では成功するが、非平滑な伝導性・透過性を持つ拡散方程式のような低次元の特殊なケースでは失敗することもある。
❺ 効果的なアーキテクチャの設計は、経験的に行われており、非常に時間がかかる。
2⃣ 五月雨式故に網羅的なPINNsの課題とその対応策：其の弐
(PGNN及び)PINNs、PeNN([*A-7]参照)、ニューラル演算子(DeepONet、フーリエ･ニューラル演算子)を含む広範なレビュー論文[*34]では、PINNsの制限(課題)として、以下をあげている(当然、㊀～㊃及び❶～❺と重複する)。なお、[*34]は流体力学と固体力学を、主要アプリケーションとして意識している。
⓵ PINNsは深いネットワークを使用するため、勾配消失の問題に直面する可能性がある。⇔㊃
👉 勾配消失に関しては、活性化関数をtanhからStanに替えるという選択がある(論文[*44]参照)←[*34]に書いてあるわけではない。
👉 他の解決策例として、(HybridNet、高密度畳み込みエンコーダ/デコーダ･ネットワーク、自己回帰エンコーダ/デコーダ･モデル、TF-Net、DiscretizationNet、及びPhyGeoNetなどの)畳み込みフィルタを使用した、離散学習技術の採用。
⓶ 実際の、3次元問題では法外に遅くなる可能性がある。
👉 ニューラル演算子(あるいはPhysics-Informedニューラル演算子)を使うという選択がある←[*34]に書いてあるわけではない。
⓷ 低次元の時空間パラメータ化に制限を課す。～❹
⓸ PINNs の収束を示す｢確実な｣理論的証明はない。⇔㊁
⓹ 重みを持つ損失関数のパーツにおける、重みを最適に選択するためのガイドラインはない。→損失関数の異なる項が学習中に互いに競合する可能性があり、この競合によって学習プロセスの安定性が低下する可能性がある。⇔❺
⓺ 高周波または、マルチスケール構造によって支配される問題の、非線形偏微分方程式を解くことができないことがよくある。⇔❷
👉 '解決策'として、Ⅰ-3【1】(3)フーリエ特徴量埋め込み、を参照←[*34]に書いてあるわけではない。
☛ [*45]は、フーリエ･ニューラル演算子(FNO)が『高周波数モードでエネルギー減衰が遅い、非線形性の高い偏微分方程式で生じる複雑な演算子を近似できる』と主張。
👉 他の解決策例として、スイス連邦工科大学チューリッヒ校の研究者による、有限基底PINNs(FBPINNs)がある[*71]。
⓻ 任意の新しいインスタンスについて、PINNs は新しいニューラルネットワークを学習する必要がある。⇔㊂
👉 ニューラル演算子(あるいはPhysics-Informedニューラル演算子)を使うという選択がある←[*34]に書いてあるわけではない。
⓼ 異種媒体、たとえば複数の材料で構成される複合スラブにおける逆問題の解決策を学習する際に困難に遭遇する。～❶
👉 解決策例として、分散PINNs(DPINNs)。
⓽ 汎化性能に問題がある。←[*34]はPeNN推しであり、PeNNは汎化性能が高いと主張。
住友電工システムソリューションと京都大学は、正確で説明可能な渋滞予測を可能とする新しい深層学習モデルQTNet を発表した(23年8月7日)。高度道路交通システム(ITS)において「説明可能性が低いブラックボックス予測は相応しくない」というモチベーションの下で開発し、東京都における1時間先の渋滞長予測で誤差40m以下を達成とアピールした[*23]。QTNetの成果をまとめた論文(以下、本論文[*24])は、8月6日に、AI分野(機械学習とデータマイニング)の国際会議「The 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining」にて発表された。
【1】本論文の主張
データ駆動型ニューラルネットワークと待ち行列理論に基づくドメイン知識を組み合わせたQTNetは、実データセットを用いた｢渋滞待ち行列の長さ(以下、渋滞長)｣予測実験で、従来手法に比べて、高い予測精度を示したと主張。深刻な渋滞においても、従来手法に比べて精度が高い。
【2】事前整理
(0) 先行研究整理[*24],[*25]
道路における自動車の渋滞予測(あるいは人流予測)という多変量時系列予測の精度を上げるには、❶流れデータの時間的ダイナミクスに加えて、❷データ間の空間依存性を考慮した予測モデルを構築する必要がある。❶単独に関しては、まず自己回帰和分移動平均モデル(ARIMAモデル)やベクトル自己回帰モデル(VARモデル)が採用された。しかし、ARIMAやVARは時系列データに定常性を仮定していたため適合性が低く、代わりに再帰型ニューラルネットワーク(LTSMやGRU)が用いられることとなった。❷単独に関しては、交通網をグラフで表して、その構造をニューラルネットワークの入力に用いるグラフ畳み込みニューラルネットワーク(GCN)が使用されている。
❶＋❷として、時空間グラフ ニューラル ネットワーク(STGNN)がある。STGNNは、交通データの時空間依存性を考慮することにより、交通予測の大幅な改善を達成した(例えばSTGNNに関するサーベイ論文として[*26]がある)。ただし、空間グラフを(機械学習モデルを使用する場合のような、自動で作成するケースと異なり)人の手を使って作成すると、交通データに固有の空間依存性を常に表しているとは限らず、予測パフォーマンスが低下する。この問題は、自己適応(型)隣接行列を使用することで解決できる。具体的には、交通データからグラフの頂点(ノード)埋め込みを推論することで、ノードの隣接行列を学習して、交通データに固有の適切な空間依存性を抽出する。
ここまでで、交通データの時間依存性と空間依存性を考慮する枠組みを用意し、さらに交通データに固有の空間依存性を取り込む工夫を凝らした。具体的な名称を上げれば、グラフWaveNetというSTGNNを採用する(詳細は後述)。残る問題は、渋滞長である。渋滞長は、信号制御のための交通需要の計算や、ドライバーへの道路状況の直感的な伝達に使用されるため、交通管理には不可欠である。しかし、これまで予測の対象になっていなかった。
(1) 物理学に基づくニューラルネットワーク(Physics-Guided Neural Network:PGNN)
本論文はPINNsとは、異なるPGNNを使っている。以下、PINNsとPGNNとの差異について説明する。なお、本論文では、物理学に基づく深層学習(PGDL)と記述されているが、PGNNと区別しない。
1⃣ 概要[*27]
PINNsは、支配方程式が存在する物理現象を扱う枠組みと考えられる。つまり、解析対象である物理系に対する知識が十分であることを暗黙の仮定としている。解析対象である物理系に対する知識が不十分であっても、その不十分な知識(部分的知識)を活用するという枠組みが、2021年に提案された、物理学に基づく機械学習(PGML)である。部分的知識とは、｢単純化されたモデルや経験則から得ることができる｣とされている。PGMLアプローチの基本的な考え方は、｢ニューラルネットワーク内の層の1つに、部分的知識を注入し、学習プロセスを導くこと｣である。
[*27]では、モデル空間をニューラルネットワークに限定したPGMLとして、PGNNを提案している。このPGNNの基本的な考え方は、主成分回帰(PCR)を”一般化する”ことらしい。PCRでは、回帰モデルへの入力として、高次元の説明変数を使う代わりに、低次元の潜在変数を使用する。こうすることで、回帰変数の複雑さを削減し、回帰分析をロバストにできる。[*27]によると、PCRには2つの大きな問題がある：①主成分分析を用いて計算された潜在変数は、説明変数の線形結合にしかなりえない、②回帰タスクが、潜在変数の計算から切り離されている。
PGNNアプローチでは、❶潜在変数は、追加の特徴(部分的知識)で補うことができる、❷潜在変数の計算と回帰の両方が、ニューラルネットワーク内で組み合わされる。
2⃣ 損失関数
PINNsは、支配方程式が存在する物理現象を扱う枠組みであり、支配方程式を損失関数として採用する。PINNsにおける支配方程式は、PGNNでは｢部分的知識｣に相当する。部分的知識を支配方程式寄りの言葉に直すと、”物理学に基づいて導出される関係式”と表現できるだろう。したがって、PGNNでは、この”物理学に基づいて導出される関係式”を損失関数として採用することになるが、一つ問題がある。PINNsでは、支配方程式を損失関数とすることで、学習モデルの出力が(少なくとも、ある程度)正しいことが保証されるのに対して、PGNNでは、物理学に基づいて導出した関係式を損失関数としたとしても、学習モデルの出力が正しいことは保証されない(はずである)。そこで、PGNNでは、残差項を導入する。つまり、物理学に基づいて導出した関係式＋残差項がPGNNの損失関数とする。👉この残差を(あくまでPINNsの枠組みとして、)深層ニューラルネットワークを補完的に使用して最小化する、という論文[*73]がある。
本論文における｢物理学に基づいて導出した関係式｣は、後述する。
3⃣ まとめ
PINNsとPGNNは、思想としては同じと考えられる。つまり、物理学に基づく関係式(ドメイン知識)をニューラルネットワークに持ち込み、学習器の出力をイイ感じに制御することで、精度を上げるという大枠は同じである。
【参考】
[*34]によれば、PGNNの"制限事項"は以下の通り。
① PGNNは物理法則に制約されていない。このため、学習データが、基礎となる物理化学的属性の全範囲をカバーしていない場合、外挿機能を提供しない。科学分野では通常、学習データは疎らである。
② PGNNは、複雑な非線形問題や、(層流から乱流まで、のように)物理化学的属性の範囲が非常に広い場合では、内挿機能(すら)も欠如する。
③ PGNNは、(一般に)逆問題に対処できない。
④ PGNNは、解像度不変ではない。つまり、低解像度で学習したり、高解像度で直接推論したりすることはできない。ただし、これは、PINNsでも同様。一方PeNN(PENN)は、解像度不変。
⑤ PGNNは、データセット全体にわたる入出力の相互依存関係を学習する。このようなプロセスでは、異なる入力と出力のペア間の関数の依存関係のわずかな変動をノイズとみなして、平均的なソリューションを生成する可能性がある。つまり、汎化性能が低い。
⑥ PGNNは、学習データセットが多様な場合、基礎となるプロセスを学習するのに苦労する可能性がある。
(2) グラフWaveNet
本論文のモデルQTNetは、グラフWaveNetの枠組みを採用しており、その改良版と考えられる。時間的畳み込みの手法として、Dilated畳み込みと因果畳み込みが使用されている。
グラフWaveNetの元となっているWaveNetは、音声という時系列データを処理するために、Googleが開発した(畳み込み)ニューラルネットワークであり、合成音声も生成できる。WaveNetは、Dilated畳み込みとCausal(因果)畳み込みが組み合わされているところに、特徴がある。Dilated畳み込みとは、｢層が深くなるにつれて、畳み込みの間隔を開ける｣という工夫(手法)で、計算を効率化する目的で行われる。因果畳み込みとは、｢過去のタイムステップだけを用いて、畳み込みを行う｣という手法で、受容野を広くできる[*28]。ちなみに、WaveNetは、洪水予測に使うことが検討されているようである[*29]。
グラフWaveNetは、グラフ構造を入力として用いるWaveNetである。そして、グラフWaveNet最大の特徴は『自己適応隣接行列(self-adaptive adjacency matrix)』をネットワークに組み込み、隣接行列自体を学習パラメータとすることで、地点間の空間的関係性を捉える、という仕組みにある[*30]。道路セグメント間の空間依存関係を抽出するために、グラフWaveNetは、まず行方向および列方向の正規化隣接行列を使用して、グラフ畳み込み演算を適用する。さらに、自己適応隣接行列 𝑨= Softmax(ReLU(𝑬1𝑬T2)) を利用して、交通量データから、固有の空間依存関係を発見する(𝑬は、対象とするノードを含む、エッジを表すベクトルのようである)。
(3) 物理学に基づいて導出した関係式：砂時計モデルと渋滞長
本論文では、｢砂時計モデル｣を仮定した上で、待ち行列理論を適用して、移動速度と交通量(交通流率)を変数として、渋滞長(待ち行列長さ)を表現する関係式を導いている。これが、本論文における｢物理学に基づいて導出した関係式｣である([*24]の著者が知る限り、このような取り組みは、これまでないそうである)。
本論文では、学習モデルからに出力における説明可能性･解釈可能性の重要性を強調している。説明可能性･解釈可能性は、次のように担保されると主張している☛ 渋滞長、移動速度、交通量を含む QTNet の出力は、待ち行列理論から導かれた｢物理学に基づいて導出した関係式｣を満たすように制約される。つまり、予測された速度と交通量は、渋滞長を正確に予測するだけでなく、予測変数の関係は互いに一致している。したがって、QTNet は完全なブラック ボックスではなく、｢物理学に基づいて導出した関係式｣に基づいて予測の説明を提供することができる。
ネットワーク･アーキテクチャ的に説明すれば『中間層が｢移動速度と交通量｣という説明可能･解釈可能な物理量を出力し、それらを使って｢渋滞長｣を出力しているから、説明可能･解釈可能が担保される』と表現できるかもしれない。
【3】シミュレーションのセットアップ
(1) データ
東京で運用されている交通管制システムによって収集された、道路セグメントごとの｢渋滞長、移動速度、交通量｣で構成されるデータが使用された。期間は、2020年 10月から2021年7月。データは 50 秒ごとに記録されているが、重大なノイズが含まれているので、適切な平滑化を行って、5 分ごとにダウンサンプリングされている。
セグメントは、都内の幹線道路を主要な交差点で分割した1,098の道路セグメントで構成されており、セグメント内に存在する渋滞は、各セグメントの下流端の交差点を起点とする渋滞として解釈される。移動速度は、セグメントを通過する車両の平均速度から計算され、交通量はセグメント内に存在するいくつかの車両感知器からの測定値の平均値(交通量の単位は、単位長さ当たりの車両台数で、交通流率とも呼ばれる)。
入力データは、過去 1 時間の｢渋滞長、移動速度、交通量｣で、1 時間先の渋滞長さが、出力である。
(2) 損失関数
QTNet の全体的な損失関数は、4つの項の和として表される：𝜆1Lq＋𝜆2Lnq＋𝜆3Ls＋𝜆4Lf。Lqは渋滞長に関する損失関数で、渋滞長に基づく重みが掛けられている。その理由は、㊀渋滞予測は、短い渋滞よりも、長い渋滞を正確に予測することが重要であり、㊁渋滞データの大部分はゼロだから、である。
Lnqも渋滞長に関する損失関数である。道路セグメントの渋滞長は、そのセグメントの長さによって上限が決まる。このため、同じ渋滞長であっても、道路セグメントの長さに応じて、異なる影響を与える可能性がある。そのため、セグメント長の逆数で重み付けした、別の損失関数Lnqが導入された。
LsとLfは、それぞれ移動速度と交通量に関する損失関数である。速度は、観測値が小さい場合の誤差を強調するために、観測値の逆数で重み付けされる。これは、値が小さいほど交通渋滞が深刻であることを示すためである。
𝜆i(i=1～4)は、渋滞長、移動速度、交通量のさまざまなスケールのバランスを取るための、重み(ハイパーパラメータ)である。ちなみに予備実験から、𝜆1=0.02、𝜆2= 4.0、𝜆3= 5.0、𝜆4= 0.007と設定されている。
(3) その他
❶ハードウェアとソフトウェア
Python 3.7.4 と TensorFlow 2.4.1 を使用してモデルを実装し、単一の NVIDIA Tesla T4 を備えたサーバーで実行した。
❷ニューラルネットネットワークの設定
初期学習率は0.001、エポック数は100に設定。オプティマイザはAdamを使用。10 エポック以内に改善が観察されない場合に学習を停止する、早期停止(Early Stopping)を採用した。
【4】シミュレーション結果
(1) 比較した学習モデル：ベースライン
HA(Historical Average：移動平均法)、VAR(自己回帰モデル)、ARIMA(自己回帰和分移動平均モデル)、LSTM(Long Short Term Memory)、DCRNN[*31]、AGCRN[*32]、グラフWaveNet、MegaCRN[*33]と、QTNetとを比較した。
(2) 評価指標
渋滞長の評価指標は、RMSE(自乗平均平方根誤差)とMAE(平均絶対誤差)である。
渋滞長は、道路セグメントと時間の大部分でゼロであるため、観測全体のパフォーマンスを評価する場合、常にゼロを予測する方法が好まれる可能性が高い。しかし、これでは深刻な渋滞を正確に予測できない。この事態を避けるため、次の4パターンで評価した。
㊀ すべてのデータ: テストデータ内のすべての観測値を使用。
㊁ ゼロ以外のデータ: 待ち行列の長さゼロを削除。
㊂ 最悪 10%: 待ち行列の長さが 200 m を超える観測。
㊃ 最悪 5%: 待ち行列の長さが 300 m を超える観測。
(3) 結果
1⃣ 渋滞長に関する総論
㊀ QTNetがRMSE、MAEともに最良。次点は、RMSEがMegaCRN。MAEがAGCRN。
㊁ RMSEはHAが最良(※)で、次点がQTNet。MAEはQTNetが最良で、次点はMegaCRN。
㊂ QTNetがRMSE、MAEともに最良。次点は共に、グラフWaveNet。
㊃ QTNetがRMSE、MAEともに最良。次点は共に、グラフWaveNet。
※ 本論文では、この理由を｢渋滞長の予測における、周期特徴の使用の重要性を示唆していると考えられる｣と整理している。ここで言う周期特徴とは、｢ラッシュアワーや週末の小旅行など、日ごと、週ごとに、渋滞の原因となる事象が繰り返し発生する｣ことを指している。
2⃣ QTNetの改善度合
ベースラインからのQTNet の改善は、渋滞長が長いほど大きい。
【5】考察など
(0) ITSという言葉を久しぶりに目にした･･･。
(1) 工学系の問題設定では、構成方程式と呼ばれる式が、しばしば登場する。教科書的に説明すると、｢物理則から導かれる支配方程式の数より、未知数の数が多い｣場合に、対象となる物理系を解析するために｢物理系の物理的性質から導出される関係式｣を構成方程式、と呼ぶ。PGNNは構成方程式を損失関数として採用するPINNsだと考えれば、肚落ちし易いかもしれない。物理系の解析において、支配方程式に加えて構成方程式を使って解析する場合もあるから、その場合はPINNs＋PGNNというハイブリッドである、という捉え方でよいのだろう。
(2) パターン思考的にL1正則、L2正則を使うのではなく、対象系の物理的性質から関係式を導いて、機械学習･深層学習の精度を上げるという手法が登場したのは、2021年と最近のようである。今後の進展が楽しみな分野だと思う。
(3) トヨタが自信をもって発表(23年2月9日)したMegaCRNよりも良い結果が得られたと主張している。計算コストの比較なども知りたいところである。
ペンシルベニア大学[*41]の研究者は、PINNsの学習効率と精度を向上させるベストプラクティスを紹介した論文[*40](以下、本論文)を発表した(23年8月16日@arXiv)。本論文は、｢PINNsにはいくつかの問題が指摘されている。しかし、それは、ネットワーク･アーキテクチャに起因するものではなく、運用上の問題である｣と主張する。
本論文の特徴は、㈠施策の中で、具体的な数値を提示している、㈡アブレーション研究をしっかり行っている、ことにあると思われる。
【1】学習効率と精度を向上させる具体的な施策
具体策は、PDE の無次元化(1)、適切なネットワーク アーキテクチャの選択(2)～(4)、および適切な学習アルゴリズム(5)～(7)の採用という3部構成である。
(1) 偏微分方程式(PDE)の無次元化
PDE の無次元化は、㊀ 一般的な初期化スキーム(例えば、Glorotスキーム)の前提を満たすため、㊁入力変数と出力変数とのスケール不均衡を緩和するため、㊂、オプティマイザーに余計な仕事をさせないため、に重要な役割を果たす。
※ NVIDIA Modulus v22.09では、PDEの無次元化が推奨されている[*42]。
(2) 幅、深さ等 ー適切なネットワーク アーキテクチャの選択1/3
① 幅と深さがそれぞれ 128 ～ 512、3 ～ 6 の範囲のネットワークを採用することを勧めている。
② 活性化関数は、tanh関数(双曲線正接関数)を勧めている。
③ 初期化スキームは、Glorot スキームを勧めている。
さらに、ランダム･フーリエ特徴量埋め込みに用いる、スケール因子σの具体的な大きさも示している。実務的には、σ ∈ [1, 10] とすることを勧めている。
(4) ランダム重み分解 ー適切なネットワーク アーキテクチャの選択3/3
ネットワーク内の各ニューロンに関連付けられた重み(ウェイト)から成る｢重み行列W｣を、指数化したスケール係数ベクトルを対角成分とする対角行列diag(exp(s))と、行列Vの積に分解する。つまり、W＝diag(exp(s))×V、と分解する。ここで、Wを構成する｢重みベクトル｣は、スカラー係数(ニューラルネットワークの文脈で言うと、学習可能パラメータ)と、ベクトルvの積に分解されている。スカラー係数から成るベクトルが、先のスケール係数ベクトルである。行列Vはベクトルvから成る。また、sは多変量正規分布からサンプリングされる。多変量正規分布の、平均(ベクトル)は0.5若しくは1.0、分散共分散行列(∑＝σ×恒等行列)のσは0.1を勧めている。
本論文は、重み行列を分解することで、損失の収束とモデルの精度が一貫して確実に向上する、と主張している。
(5) 時間的因果関係の尊重 ー適切な学習アルゴリズムの採用1/3
”時間依存”PDE を解くときに PINNs が時間的因果関係に違反する可能性があるため、誤った解に収束しやすいことが示されている。これは、PINNs がすべての PDE残差を同時に最小化する傾向がある一方で、初期の正しい解を取得する前であっても、後の時点で PDE残差を最小化する方向に望ましくない偏りがあるためである。ここで、PDE残差とは『ネットワークからの出力が、(未だ)正解でないために支配方程式を満たさないために、発生する誤差(残差)』を指している。
対応策として、時間軸に沿った重み｢時間的重み｣を導入する。最適化プロセス内で欠落している因果構造を強制するために、まず時間領域を分割する。次に、時間領域の i 番目セグメント内の損失を導入する。この損失は、支配方程式から(PDE残差から)の損失である。時間領域全体の損失は、{i番目の時間的重み×i 番目セグメント内の損失}を iについての総和したものである(数学記号で表せば、∑i)。
※ NVIDIA Modulus v22.09では、｢causual weighting スキーム｣によって、物理的因果関係に違反する連続時間PINNsが内包するバイアスに対処している。
学習中の損失のバランスを自動的に調整するスキームとして、1⃣勾配ベースの重み付けスキームと 2⃣NTK (Neural Tangent Kernel神経正接核)ベースの重み付けスキームを提案･提示している。NTK の固有値は損失関数の収束率を特徴付ける。固有値が高いほど、収束速度が速いことを意味する。
２つのスキームのパフォーマンスは似ているが、 1⃣で更新された重み付けは、2⃣に比べて安定していない。一方、2⃣のスキームではより高い計算コストが必要である。結論として、通常は勾配ベースのスキームを最初の選択肢として採用すること勧めている。
(8) その他
① Adam オプティマイザーは、大幅なチューニングを行わなくても、一貫して優れたパフォーマンスを発揮する。
② 重み減衰の使用は、勧めない(特に順方向問題に対して)。結果の予測精度が低下する傾向があるためである。[☞つまり、AdamWは推奨されない？]
③ 初期学習率0.001と指数関数的減衰を組み合わせると、通常は良好な結果が得られる[例えば、本論文では、事例/アブレーション研究で、(ほとんど)2,000 減衰ステップごとに、減衰率 0.9 の指数関数的減衰を行っている]。
④ フルバッチ勾配降下法を使用して PINNs を学習すると、過学習する可能性がある。したがって、ランダムサンプリングの使用を強く勧める。
⑤ 周期的な境界条件をハード制約として厳密に課すと、学習の収束と予測精度が大幅に向上する。
⑥ MLPのアーキテクチャを修正する(詳細割愛)。
【2】為参考：アブレーション研究
(1) 前説
本論文では、包括的で完全に再現可能なアブレーション研究が、豊富な対象に対して提示されている。対象は、❶1次元アレン･カーン方程式[反応拡散方程式：PINNsの代表的な適用例として知られている]、❷１次元移流方程式[輸送現象を記述する方程式]、❸ストークス流れ[厳密解が存在する遅い流れ]、❹蔵本･シバシンスキー方程式[化学乱流を記述する方程式]、❺蓋駆動キャビティ流れ[レイノルズ数は3,200]、❻(トーラス上の)非圧縮性ナヴィエ･ストークス方程式[流体の基礎支配方程式]、❻(円柱周りの)非圧縮性ナヴィエ･ストークス方程式、である。評価項目は、 L2ノルムによるエラーと実行時間である。❻と❼は、それほど芳しい成果が得られていないので、❶～❺について、本論文で提案した施策の有効性を評価している。なお、❺は、カリキュラム学習を実施している。まず、層流(レイノルズ数100、400、1,000)から学習を進めていき、最終的にレイノルズ数3,200の乱流を学習している。❻と❼のレイノルズ数は100。
共通のセットアップは、以下の通り：4 つの隠れ層、各隠れ層あたり 256 個のニューロンのMLP を採用。初期化スキームは、Glorot スキーム。 Adamオプティマイザーを使用。(初期)学習率 10−3 から開始し、減衰率 0.9 の指数関数的減衰を行う。
異なるセットアップは以下の通り：活性化関数は、❸ストークス流れを除いてtanh(ストークス流れは、GeLU)。損失再重み付けスキームは、❶を除いて、勾配ベース(つまり、❶は神経正接核(TNK)ベース)。学習率の減衰ステップは、❶と❺を除いて2,000毎(❶は5,000、❺は10,000)。コロケーション･ポイント(計算点)は、❹を除いて、バッチサイズ8,192の計算ドメインから均一に(ランダムに)サンプリングした(❹は、4,096)。
(2) アブレーション研究の具体的な結果
(1)の❶～❺を対象に、L2ノルムエラーを評価指標として、最も効いた施策を、まとめた。
❶→ (3)フーリエ特徴量埋め込み
❷→ (6) 損失再重み付けスキーム
❸→ (4)ランダム重み分解
❹→ (6) 損失再重み付けスキーム
❺→ (3)フーリエ特徴量埋め込み
(注：2年以上前の古い論文・・・)米ブラウン大学及び中国･上海大学他[*65]は、｢予測精度が高い疫学モデルを選択するフレームワークの提供を目的とする｣研究論文(以下、本論文[*66])を発表した(2021年11月@nature computational science)。本論文は、二つの捉え方ができる。
パラメータを、｢感染からの経過時間に依存しない｣一定のスカラー値に固定した疫学モデルは、予測精度が低いであろうことは、容易に予想できる(実際、新型コロナウィルスの事例で散々、批判に晒された)。固定パラメータではなく、経過時間に依存するパラメータ(時間依存パラメータ)として捉える場合、予測精度は向上することが期待される。連立一階常微分方程式系で記述される疫学モデルの時間依存パラメータを推定することは、いわゆる逆問題に相当する。逆問題解析は、残差(ノルム)の最小値を探索する最適化問題に帰着する。本論文は、この最適化問題を、PINNs(Physics-Informedニューラルネットワーク)を使って解いた論文という捉え方ができる。PINNsを使った理由は、｢逆問題が得意｣であり、｢外挿｣を可能とするからである。
二つ目は、連立非整数階常微分方程式系で疫学モデルを構築し、予測精度を議論している論文という捉え方である(もちろんPINNsは使用する)。
【1】本論文の結論
結論は、以下の通りである。
(1) PINNsを利用することで、少なくとも短期(2週間)予測においては、時間依存パラメータを外挿できることがわかった。
(2) 整数階微分モデルは、非整数階微分モデルよりも堅牢である。これは、次の意味である：利用可能なデータが不十分な場合†、非整数階数モデルによる予測は、比較的大きな不確実性を示す。利用できるデータが増えれば、不確実性は減少する。
(3) 長期(3か月)予測の精度が最も高かったのは、9つの疫学クラスのダイナミクスを考慮するモデルであった。
† パラメータの識別可能性が欠如している場合。観測される変数分布から、パラメータが一意に定まるとき、｢パラメータは識別可能である｣と言われる。
【2】事前整理
(1) 疫学モデル
0⃣ 概要
疫学(epidemiology)とは、時空間的な病気の広がりと、それに関連する要因を調べることであり、感染者が宿主集団内にどのように侵入し、流行を引き起こすのかを個体群動態としてとらえる分野である[*67]。疫学モデルは、特に長期にわたるパンデミックの場合、疾患の経過全体にわたる疾患の症状、病原体の濃度および挙動の変化による｢時間依存パラメーター｣に対応する必要がある。
1⃣ 疫学9クラスと3パラメータ
本論文に現れる 疫学9クラス(なお補助的に、累積I及び累積Hが考慮されている)と3パラメータについて説明する。
まず、9つの疫学クラスを定義する。
㊀Susceptible(感受性と訳される)：感染していない(ので、健常者とも呼ばれる)、未感染。しかし、免疫を持っていないため、今後、感染する可能性がある者。㊁Expose(暴露と訳される)：感染したが、潜伏期間中であるために発症していない者。感染源とはならない。㊂Pre-symptomatic(発症前と訳される)：感染したものの、まだ発症していない者(今後、発症することが推測されている)。発症していないが、感染源となる。㊃Quarantined(社会的隔離と訳される)：感染して発症しているが、隔離されているため、感染源とはならない(という想定)。㊄Infectious(感染と訳される)：感染者。感染して発症し、感染源ともなる。㊅Asymptomatic(無症状と訳される)：感染したが、発症しなておらず、今後も発症しないと推測される者。感染源となる 。㊆Hospitalized(入院と訳される)：感染して発症しているが、隔離されているため、感染源とはならない(という想定。院内感染の可能性はある？)㊇Disease-related death(訳語は、不明)：感染を事由として亡くなった者。㊈Recovered/removed(回復/隔離と訳される)：回復は、感染した後に免疫を獲得して、症状が収まった者。感染Iから、一定の回復率(一般に、γと表記される)で回復する。今後、感染しないし、感染源ともならない。
次に、3つのパラメータを説明する。
⓵市中感染率βI(t)、⓶入院者の中で感染事由で死亡した患者の割合q(t)及び、⓷感染･発症した者の中で入院した者の割合p(t)。補足すると、βI(t)は、接触が発生した場合の一人当たりの感染率である。
2⃣ SIRモデル及び、その拡張
もっともシンプルな疫学モデルとして、SIRモデルがある。SIRは、上記1⃣で言うと㊀、㊄、㊈である。SIRモデルの前提は、㈠一度感染したら、(免疫を獲得し)二度と感染しない。㈡感染期間中のみ、他者に感染させる能力を有する、㈢感染者は一定の割合で回復する、である。㈡を精緻にすると、1⃣における㊁、㊅が現れ、SEIRモデルやSEIJRモデルとなる。
また、SIRモデルは、感染者の割合(及び回復者の割合)を予測するモデルと考えられるが、入院患者の割合や死亡者ぼ割合を予測したければ、SIHDR(あるいはSEIHDR、SEIJHDR)モデルを考慮することになる。
本論文では、7つのモデル(厳密に言うと9つ)を扱っている。詳細は、【3】(3)を参照。
(2) 非整数階微分[*68],[*69]
1⃣ 概要
本論文で用いられているfractional derivativesあるいはfractional orderのfractionalに対する訳語は、[*68]にならって、非整数階をあてた(fractional derivativesは、非整数階導関数)。非整数階微積分法の研究は、ライプニッツの時代から(つまり創世記から)始まっていたらしいが、フラクタル解析の必要性から研究が進展したようである。フラクタル構造体に対して解析的な扱いを行うためには、非整数のフラクタル次元に基づいた、非整数階微積分が必要となる。なお、非整数階微積分法の研究主に、旧ソ連や東欧圏で行われてきたようである。
非整数階微分の定義は、｢有限差分を基礎にした定義、積分を基礎にした定義、それ以外の定義｣の3つに大別されるそうである[*69]。本論文では、 積分を基礎にした定義に該当する｢Caputoの非整数階微分｣ が使用されている。加えて本論文では、非整数階微分の階数が、定数ではなく変数(パラメータ化)されている(ただし、パラメータ化するという拡張は、本論文のオリジナルではない)。
2⃣ 履歴の記述
整数階の微積分は、点における値で計算可能であるため、過去の履歴は無関係である。一方、非整数階の微積分は、過去からの履歴を積み上げた上に計算される。これを逆に見ると、非整数階微積分は、履歴現象を(より良く)記述できる、と捉えることができる。履歴現象で思い浮かぶのは、物性物理であれば、ヒステリシス。金融(金融工学)で言えば、経路依存のデリバティブ(ルックバックやアジアン、バミューダン等)が思い浮かぶ。
履歴→記憶効果とやや概念拡張すると、非整数階微(積)分の範囲が広がる。本論文では、新型コロナウィルスの拡散過程において、記憶効果が存在する(べき乗則スケーリングが確認されている)として、非整数階微分演算子を使った常微分方程式系で構成した疫学モデルも評価検討している。
3⃣ 蛇足
PINNsでは、各種ライブラリが装備する自動微分を使うことが標準的である。しかし、自動微分は、非整数階微分には対応していない(ので、使えない)。計算コストが高いので、削減する工夫も行われている。
【3】セットアップ
(0) アーキテクチャ
1⃣ モデル･アーキテクチャ
1階常微分方程式系で支配される”物理系(拡散過程)”に対して、ニューラルネットワークを使った代理モデルを構築しているという見方ができる。また、ガウス過程が使われることが多いブラックボックス最適化問題を、ニューラルネットワークを使って解決しているという見方もできる。
2⃣ ネットワーク･アーキテクチャ
PINNsの定式化では、活性化関数がtanhで、層毎に32 個のニューロンを持つ隠れ層10層のニューラルネットワーク(NN)を使用する。時間依存パラメーターは、活性化関数がtanhで層毎に20 個のニューロンを持つ隠れ層5層の個別NNによってパラメーター化される。
(1) データセット
0⃣ 共通
データセットは、 7 日間の移動平均ウィンドウを適用することで前処理され、感染発生の報告における平日と週末の変動を平滑化した。
1⃣ ミシガン州のデータ
ミシガン州データセットは、2020年4月15日以降の毎日の感染症例 Inew、現在の入院症例 H、および毎日の死亡症例 Dnewで構成される。
2⃣ ロードアイランド州のデータ
ロードアイランド州のデータ セットは、2020年3月4日以降の毎日の感染症例Inew、毎日の入院症例 Hnew、現在の入院症例 H、および毎日の死亡症例 Dnewで構成される。
3⃣ ニューヨーク市のデータ
ニューヨーク市のデータセットは、(2020年4月1日以降の？)累積感染者数Ic、累積入院者数Hc、累積死亡者数Dc、で構成される。開始日の記載は明確ではない。
4⃣ イタリアのデータ
イタリアのデータセットは、(2021年3月1日以降の？)毎日の感染者数Inew、回復者数Rnew、死亡者数Dnewで構成される。開始日の記載は明確でない。
(2) パラメータ
βI：市中伝染率
p：H(入院)クラスの疾病関連死者の割合
q：入院患者の割合
r：感染率
a：回復率
b：死亡率
d：ムダ時間モデルのムダ時間(日数)
ε：軽症感染から重症感染への移行割合 0.75
δ：無症状者の割合 0.6
χ：未症状者から発症者への移行割合 0.55
γ1：隔離割合 0.0365
1/α：平均潜伏期間 5.2日
1/α1：潜伏期間 2.9日
1/α2：未症状期間 2.3日
1/γα：無症状者の平均感染期間 6日
1/γ：症状者の平均感染期間 6日
1/ΦR：症状からの回復までに入院している期間 7.5日
1/Φd：死亡するまでに入院している期間 15日
1/dI：隔離までに感染している平均期間 2.9日
1/dH：平均入院期間 6.9日
1/dQ：平均隔離期間 10日
(3) 本論文で検討されたモデル
1⃣ 整数階微分モデルI1～I3
整数階微分モデルは、モデルに含まれるパラメータを、観測データを再現するように、適合させる。
パラメータは、(2)を参照。なお、正確には、非整数階モデルのF1とF3の階数を整数化したモデルI4とI5も整数階微分モデルに含まれるが、割愛。
❶I1：整数階SEIJHDRモデル。モデルのダイナミクスは、SEIJDHRと累積Iと累積Hを含めた9変数に対する9組の1階常微分方程式系で表現される。実効再生産数†は、βI×((1－δ)ε/γα＋δ/γ)である。変数記号については、(2)パラメータを参照。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→E→I→H→D、S→E→I→H→R、S→E→I→R、S→E→J→R、S→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,J,S)、E’←(I,J,S,E)、I’←(I,E)、J’←(J,E)、D’←(H)、H’←(I,H)、R’←(I,J,H,S)、IC’←(E)、HC’←(I)。ここで、例えば、S’←(I,J,S)は、Sの導関数(時間に関する1階微分)が、I,J,Sという３つの変数を使った多項式で表現できることを表している(プライムは、ドットの代わり)。
❷I2：整数階SEPIJDHRモデル。ダイナミクスは、SEPIJDHR＋累積I＋累積Hの10変数に対する10組の1階常微分方程式系で表現される。実効再生産数は、βI×((1－δ)ε/γα＋δ/γ＋χ/α2)。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→E→P→I→H→D、S→E→P→I→H→R、S→E→P→I→R、S→E→P→J→R、S→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,J,S,P)、E’←(I,J,S,E,P)、P’←(P,E)、I’←(I,P)、J’←(J,P)、D’←(H)、H’←(I,H)、R’←(I,J,H,S)、IC’←(P)、HC’←(I)。
❸ I3：整数階SEPIJDHQRモデル。ダイナミクスは、SEPIJDHR＋累積I＋累積Hの11変数に対する11組の1階常微分方程式系で表現される。実効再生産数は、βI×((1－δ)ε/γα＋δ/dI＋ε/α2)。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→E→P→I→H→D、S→E→P→I→H→R、S→E→P→I→Q→R、S→E→P→J→R、S→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,J,S,P)、E’←(I,J,S,E,P)、P’←(P,E)、I’←(I,P)、J’←(J,P)、D’←(H)、H’←(I,H)、Q’←(I,Q)、R’←(J,H,S,Q)、IC’←(P)、HC’←(I)。
† ｢ある時刻tにおける、1人の感染者による二次感染者数｣と定義される。
2⃣ 非整数階微分モデルF1～F3
非整数階微分モデルではパラメータではなく、非整数階微分の階数を、観測データを再現するように、適合させる。SIHDRのそれぞれに対して、階数κi(i=1～5)を適合させる。
❹F1：非整数階SIRモデル。SIR＋累積Iの4変数に対する4組の非整数階常微分方程式系でモデルのダイナミクスが表現される。実効再生産数は、βI/γ1。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→I→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,S)、I’←(I,S)、R’←(I)、IC’←(I,S)。❹～❻では、プライム記号’は、時間に関する非整数階微分を表すこととする。
❺F2：非整数階SIDRモデル。SIDR＋累積Iの5変数に対する5組の非整数階常微分方程式系で表現される。実効再生産数は、r/b。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→I→D、S→I→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,S)、I’←(I,S)、D’←(I)、R’←(I)、IC’←(I,S)。
❻F3：非整数階SIHDRモデル。SIHDR＋累積I＋累積Hの7変数に対する7組の非整数階常微分方程式系で表現される。実効再生産数は、βI/γ。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→I→H→D、S→I→H→R、S→I→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,S)、I’←(I,S)、H’←(I,H)、D’←(H)、R’←(I,H)、IC’←(I,S)、HC’←(I)。
3⃣ 時間遅延モデルT1
整数階微分モデルに、新型コロナウィルスの拡散過程における記憶効果を考慮したモデル。
❼T1：ムダ時間SIJHDRモデル。I1からEを取り除き、Sにムダ時間を導入(S(t)をS(t-d)とする。dが｢ムダ時間｣)したモデル。実効再生産数は、βI×((1－δ)ε/γα＋δ/γ)で、I1と同じ。
各クラスの関係をグラフ構造(非巡回有向グラフ)的に表現すると、S→I→H→D、S→I→H→R、S→I→R、S→J→Rとなる。当該モデルの常微分方程式系は、以下のように表現できる：S’←(I,J,S)、I’←(I,J,S)、J’←(I,J,S)、H’←(I,H)、D’←(H)、R’←(I,J,H)、IC’←(I,J,S)、HC’←(I)。❼におけるプライム記号は、時間に関する1階微分を表す。
【4】検証結果[*70]
(1) 短期予測における外挿性の確認
ニューヨーク市のデータセットと、整数階モデルI1～I3を使って、短期予測の精度を比較した。正確には、P1～P3という３つのアプローチで精度を比較した。
P1は、I1～I3を使う。推定したパラメーターβI(t)、p(t)、及びq(t)は、学習最終日(本稿では仮にTとする)における値で固定する。すなわち、 βI(T)、p(T)、及びq(T)とする。これらのパラメータを使った常微分方程式で、I(市中感染率)、H(入院割合)、D(死亡割合)に対して2週間後の将来値を計算し、実データと比較する。
P2は、I3のみを使う。I3で推定したパラメータを外挿する。外挿したパラメータを使った常微分方程式で、I、H、Dの将来値を計算し、実データと比較する。
P3は、代理モデルとしてのPINNsで将来予測値を計算する。
P2の精度が最も高かったため、短期(約2週間)予測では、パラメータの外挿が可能と結論している。
(2) 整数階モデルの堅牢性
ニューヨーク市のデータセットを使って、非整数階モデルF1と、F1を整数階化したI4モデル(実は、単なるSIRモデル？)を比較している。データが欠けている状態で出力する｢S、I、R｣を、説明の都合上本稿では、｢S0、I0、R0｣としよう。データを徐々に増やしていく状態で出力する｢S、I、R｣を｢Sj、Ij、Rj｣としよう。j=1,2,3及び0とできる。
添え字jを0から増やしていっても、I4モデルの出力は大きく変動しない。一方、
F1の出力は、変動が大きい。このことを持って、整数階モデルの方が、データにまつわる問題があっても(データが欠落している状態であっても)、堅牢である、と結論している。
(3) 長期予測
ニューヨーク市のデータセットを使い、(1)におけるアプローチP1を採用して、長期(約3か月)予測を行った。I、H、Dの将来値と実データを比較した結果、I3モデルの精度が最も高かった。
【5】考察
(1) ｢逆問題解析に適する｣、｢外挿可能｣というPINNsの2大特徴が発揮された研究といえるのだろう。
ただ残念ながら、短期予測ではパラメータ外挿が可能であった一方、長期予測では難しかったのだろう。固定パラメータで実施している。
(2) グラフ構造(非巡回有向グラフ)を入力とする、グラフ･ニューラルネットワークを適用すると、どうような結果(予測精度)がでるだろうか。
(3) なお、ブラウン大学はMITと共同で、ニューラル演算子を適用して、パンデミックに関する予測(感染者割合の予測)を行った研究成果を発表している(Ⅳ-2)。ベイズ的実験計画法＋ニューラル演算子で、極端現象の定量化を行うという研究であるが、少ないサンプル数で、良いモデルが得られた、と結論している。ちなみにモデルは、SIRモデルである。
(4) 非整数階モデルも、残念ながら活躍できていないが、そもそも記憶効果は、どれほど重要なのであろうか。記憶効果の重要性を明らかにするために、ムダ時間モデルの深堀りを行っても良かったのではないだろうか。思ったほど、記憶効果が重要でないのであれば、そもそも非整数階モデルは必要とならないはず。
中国･清華大学と米ペンシルベニア大学の研究者は、PINNs(Physics Informedニューラルネットワーク)の課題とその解決策について、まとめた論文[*74](以下、本論文)を発表した(2023年10月5日@arXiv)。
｢広範な偏微分方程式(PDE)にわたるPINNs の包括的な比較が不足している｣という問題意識の下、行われた研究であり、ベンチマークツール”PINNacle”という形で、実装化されている。
【1】全体構成
(1) 以下4項目への対応を、PINNsの課題として、特出しする：
1⃣複雑形状を有する問題、2⃣マルチスケール問題、3⃣高次元問題、4⃣非線形問題
(2) 課題解決策を、6つに分別する：
⓵重みの変更/再サンプリング、⓶オプティマイザー、⓷変分法に基づく損失関数、⓸正則化を導入した損失関数、⓹ドメイン分解、⓺適応型活性化関数
※ ⓵は1⃣と2⃣を、⓶は2⃣を、⓹と⓺は1⃣と2⃣を、を解決もしくは、経験的な改善を示したりするように設計されている可能性がある。
(3) 7つの物理系と逆問題(支配方程式22個)を、検証対象として選択する。
(4) 課題解決策(2)を実装したPINNsを選択し(3)に適用することで、解決策の有効性を検証する。
【2】選択された物理系(支配方程式)
(1) バーガース方程式・・・2個 (※バーカースに表記を統一)
①1次元バーガース方程式、②2次元バーガース方程式
(2) ポアソン方程式・・・4個
③2次元ポアソン方程式、④不規則形状を持つ2次元ポアソン･ボルツマン方程式、⑤複雑形状を持つ3次元ポアソン方程式、⑥多数の部分領域を持つ2次元ポアソン方程式
(3) 熱方程式・・・4個
⑦変数係数の2次元熱方程式、⑧マルチスケールの2次元熱方程式、⑨複雑形状を持つ2次元熱方程式、⑩時間発展が長時間に及ぶ2次元熱方程式
(4) ナヴィエ･ストークス(NS)方程式・・・3個
⑪2次元蓋駆動キャビティー流れ(レイノルズ数Re＝100)、⑫2次元バックステップ流れ(Re＝100)、⑬時間発展が長時間に及ぶ2次元NS方程式
※ NS方程式は非線形偏微分方程式であるが、Re＝100程度では、非線形性は弱い。
(5) 波動方程式・・・3個
⑭1次元波動方程式、⑮不均一媒体中の2次元波動方程式、⑯時間発展が長時間に及ぶ2次元マルチスケール波動方程式
(6) カオス的挙動を示す偏微分方程式・・・2個
⑰2次元反応拡散方程式(グレイ=スコット･モデル、GS方程式)、⑱蔵本シバシンスキー方程式(KS方程式)
(7) 高次元偏微分方程式・・・2個
⑲5次元ポアソン方程式、⑳5次元熱方程式
(8) 逆問題・・・2個
㉑2次元ポアソン方程式の逆問題、㉒2次元熱問題の逆問題
‖補足‖4つの課題に対する振り分け
1⃣＝④、⑤、⑨、⑫、⑮。2⃣＝⑥、⑧、⑯、⑰、⑱。3⃣＝⑲、⑳。4⃣＝⑬、⑰、⑱、㉑、㉒。
【3】選択されたPINNs
(1) 解決策を実装したPINNs
⓵重みの変更†1/再サンプリング†2
3種類のPINNsを選択した：勾配ノルムを使用して重みを変更する❹PINN-LRA、神経正接核を使用して重みを変更する❺PINN-NTK。そして、残差ベースで再サンプリングを実施する❻PINN-RAR。RARとは、残差ベースの適応的リファインメント(Residual Adaptive Refinement)の意味である。学習中にコロケーション点の分布をサンプリングし直して、適応的な改善を行う。
†1 ｢PDE損失、境界損失、データ損失｣という3項目から構成される、PINNsの損失関数における各項目の重みを、(適応的に)変更することを指している。一般に、ニューラルネットワークが支配方程式を正確に近似できないことから生じる損失をPDE(偏微分方程式)損失と呼ぶ。初期条件/境界条件とコロケーション点のズレにより発生する損失は、境界損失と呼ばれる。
†2 コロケーション点を変えてサンプリングすることを指す。結果として、重みの変更につながる。
⓶オプティマイザー(最適化器)
新しいオプティマイザーとして、❼MultiAdamを採用している。MultiAdamは、勾配の運動量を利用して損失項をパラメータ毎にバランスさせる(スケール不変の)オプティマイザーである[*75]。
⓷変分法に基づく損失関数
❽hp-VPINNが選択されている。hpはh-リファインメント＋p-リファインメントを意味している。VPINNのVは、変分法(Variational method)のVである。h-リファインメント、p-リファインメントは有限要素法(FEM)における用語で、h-リファインメントはメッシュ間隔を狭めることで、FEMの計算精度を上げるというアプローチを指している。p-リファインメントは高次の補間関数を使用することで、FEMの計算精度を上げるというアプローチを指している。hp-VPINNはFEM(正確には、ペトロフ･ガラーキン法を採用したFEM)にインスパイアされている。
⓸正則化を導入した損失関数
❾gPINNが選択されている。gPINNはPDE残差の勾配情報を利用し、勾配を損失関数に埋め込んでいる。gPINNのgは勾配(gradient)のgである。
⓹ドメイン分解
ドメイン分解は、ドメイン全体を複数のサブドメインに分割し、これらのサブドメイン上でサブネットワークを学習する方法である。マルチスケールの問題を解決するのに役立つが、複数のサブネットワークがあると学習の難易度が上がる。XPINN、cPINN、および FBPINNの3つが代表例である。本論文では、｢学習を安定させるために、ドメイン固有の正規化を適用する最先端のドメイン分解である｣❿FBPINNが選択されている(FBPINNは[*71]も参照)。ちなみに、FBPINNは差分法(有限差分法)にインスパイアされている。
⓺適応型活性化関数
⓫LAAF(Locally Adaptive Activation Function)と⓬GAAF(Global Adaptive Activation Function)が選択された。
(2) 比較対象となったPINNs
以下の3つがベースラインとなった：(プレーンな)❶PINN、境界損失の重みを大きくして学習したPINNである❷PINN-w(→⓵に該当する)、最適化アルゴリズムにL-BFGS法を使った❸PINN-LBFGS(→⓶に該当する)。L-BFGSと法は、制限付き-ブロイデン･フレッチャー･ゴールドファーブ･シャンノ(Broyden-Fletcher-Goldfarb-Shanno)法である。
【4】PINNsのセットアップ
(1) モデル･アーキテクチャ
ニューラルネットワークのアーキテクチャは、5つの層を持つ多層パーセプトロンである。各層のニューロン幅は100である。エポック数は20,000。この数は、エポック回数が多すぎる場合に起こり得る過学習を回避できる水準である。
コロケーション点の数は、2次元問題では8,192、3次元問題では32,768。境界/初期点は、2次元問題では2,048、3次元問題では8,192。
(2) ハイパーパラメータ
Adamオプティマイザーのハイパーパラメータ(1次モーメント及び2次モーメント用の係数)は、(0.9,0.999)を使用。学習率は１×10-3に設定。減衰させずに固定。
(3) その他備考
PINN-wは、PDE損失の重みを1、境界損失を100としている。PINN-LRAは、重みの更新にα=0.1を設定するが、これは原論文で推奨されている値である。gPINNでは、勾配項に0.01の重み、それ以外には1の重みを割り当てる。
PINN-RARでは、残差が最大となる新しい点を2000エポックごとに学習セットに追加する。
1次元波動方程式⑭では、境界損失を｢ディリクレ境界損失と非ディリクレ損失｣に分け、10,000エポック学習させた。
【5】検証結果
(1) 結果の整理
モデル(❶～⓬)の性能は、L2相対誤差で計測した(実際は、L1相対誤差を含めて、いくつかの指標が算出されているが、モデルの比較検証はL2相対誤差で行われている)。相対誤差は、グランドトルゥースとPINNsによる近似解との差である。グランドトルゥースは、有限要素法ソルバーのCOMSOLを用いて得られた数値解析解である。ただし、GS方程式とKS方程式にはChebfunを使用している(Chebfunはチェビシェフ関数Chebyshev functionの意味)。
表記1⃣～4⃣と⓵～⓺の内容については【1】を、①～㉒は【2】を、❶～⓬は【3】を参照。また、以下では、⇒は、最も優れたモデルであることを意味している。→は2番目の優れたモデルであることを意味している。
⑩、⑬、⑯、⑱については｢優れたモデル該当なし｣という結果になっている。これは、全モデルのL2相対誤差が、およそ1.0×100であったことを意味している。
❶なし
❷ →②、→⑲(3⃣)
❸⇒①、⇒⑧(2⃣)、⇒⑳(3⃣)
⓵❹⇒②、⇒⑤(1⃣)、⇒⑦、⇒⑲(3⃣) →㉒(4⃣)
⓵❺⇒③、⇒④(1⃣)、⇒⑭ →⑦、→⑧(2⃣)
⓵❻ →⑨(1⃣)
⓶❼⇒⑥(2⃣) →③、→⑤(1⃣)、→⑭、→⑰(3⃣4⃣)、→⑳(3⃣)
⓷❽⇒㉑(4⃣)、⇒㉒(4⃣) →⑮(1⃣)
⓸❾ →㉑(4⃣)
⓹❿⇒⑰(3⃣4⃣) →④(1⃣)
⓺⓫⇒⑨(1⃣)、⇒⑪、⇒⑫(1⃣) →①、→⑥(2⃣)
⓺⓬⇒⑮(1⃣) →⑪
(2) 本論文における検証結果の評価
本論文では、以下のように総括している。
⓵(つまり❹～❻)は、1⃣複雑形状を有する問題(④、⑤、⑨)や、2⃣マルチスケール問題(⑥、⑧、⑯)を含む、一部のケースで性能が向上している。この理由は『学習プロセス中に、より困難なドメインに、より多くの重みを適応的に配置した結果と考えられる』。一方で『これらの方法は、ナヴィエ･ストークス方程式(⑪～⑬)、波動方程式(⑭～⑯)、および3⃣高次元(⑲、⑳)、時間発展が長時間に及ぶ場合(⑩、⑬、⑯)には依然として困難を伴う』。
MultiAdam(❼)に対する本論文の評価は『いくつかの単純なケースとGS方程式(⑰)を解決するが、他の方法を大幅に上回るパフォーマンスはない』である。
変分法に基づく損失関数(❽)の評価は『逆問題の解決において顕著な優位性を示す』が、『順問題は、標準PINN を超える明確な改善は見られない』。
FBPINN(❿)は『GS 方程式で最小の誤差を達成している』。
(3) 検証結果総括
㈠ 本論文では『⓵重み変更(❹～❻)あるいは、⓹ドメイン分解(❿)は、1⃣複雑形状を有する問題と2⃣マルチスケール問題とに対する性能を、向上させる可能性がある』とある。👉⓵は同意だが、⓹(❿)は、その結果が出ていないと思われる。
㈡ 変分法に基づく損失関数(⓷❽)は、逆問題で良い性能を発揮する。👉これは、その通りであろう。
㈢ 本論文では『4⃣非線形問題に適切に対処できる方法はほとんどない』とある。👉これは、KS方程式に対して芳しい成果が出せなかったことを反映しているのであろう。GS方程式に対しては、まずまずの成果が出ている。
【6】考察
(0) 本論文の内容は、包括的なベンチマークで、大変有用だと思われる。
(1) 22件中4件(⑩、⑬、⑯、⑱)で、全モデルが機能しなかった。これはPINNsにとって、やや厳しい結果と言えるだろう。
(2) とは言いながらも、1⃣複雑形状に対しては、｢損失関数の重みを調整する手法❹と❺｣が有効と言えるだろう(本論文の結果をそのまま鵜吞みにすれば、結果的に"重み調整"につながる"再サンプリング"❻は、わざわざ選択するまでもないという結論で良いだろう)。加えて、適応型活性化関数⓫と⓬も有効と言えるだろう(本論文は、⓫⓬と❿を間違えているのではないだろうか)。
(3) 2⃣マルチスケール問題に対しては、オプティマイザーが有効であるように思える(❸🖋と❼)。本論文は、❹～❻及び❿が2⃣に対しても有効という書きぶりであるが、これも、間違えではないだろうか(註:❿は元々、マルチスケール対策として提案されている)。なお、マルチスケール問題対策として、勾配ブースティングをPINNsに組み合わせるというアプローチもある[*76]。
(4) 最大の問題は、やはり非線形問題に弱いということであろうか。[*40]では、乱流(高Re数)でNS方程式を解く場合、低Re数から学習を始める｢カリキュラム学習｣を提案している(Ⅰ-3 PINNsのベストプラクティス【1】(7)を参照)。なお、[*77]では(他の論文を引用する形で)カリキュラム学習が、｢時間発展が長時間に及ぶ｣場合に有効である("症状"を軽減する)ことを述べている。
(5) ちなみに｢損失関数の重みを調整する手法｣は、Ⅰ-3 PINNsのベストプラクティス【1】(6)でも扱っている。そこでは、本稿で表記するところの❹と❺が登場している。[*40]によれば、❺は❹より安定しているが、計算コストが高いので、❹を推奨している。
(6) ｢❹あるいは❺による重みを調整法｣＋｢最適なオプティマイザー(とりあえずL-BFGS🖋を使用)｣で対応するというのが、現実解になるだろうか(ちなみに、｢損失関数の重みを調整する｣手法単体では、高次元問題における困難を解消できないが、オプティマイザを含めると、解消できていると思われる)。本論文では、❹or❺＋❸or❼という合わせ技は披露されていないが、今後登場するのだろう(期待込み)。
(7) 非線形問題は、そもそも難しい。従って、PINNsの弱点という見方をするのではなく、PINNsにとっても難しいという解釈で良いのではないだろうか。
(8) 本論文ではPINNsに対して、4つの課題を提示した。これ以外にも、㊀計算時間が長い、㊁最適化が保証されない(しばしば収束しない)、㊂解析対象は小さい領域に限定されている、等がある。
㊀の対策として、例えば[*77]は新しいオプティマイザー｢エネルギー自然勾配降下法(E-NGD)｣を提案している(→下記Ⅰ-6を参照)。例えば2次元ポアソン方程式を解くケースで、最適化に要する時間は、Adamで約1時間(原論文ママ)のところE-NGDでは43秒であった(ちなみに精度もE-NGDが高い(というより、[*77]では精度を主問題にしている))。1次元熱方程式ではAdam約2時間に対し、E-NGD2分48秒。いずれにしても、オプティマイザーは重要であろう(PINNsに限らず、NN全体に言えるだろうが･･･)。
㊂の対応策には、ドメイン分解～FBPINN(有限基底PINN)[*71]がある。
🖋 本論文では、L-BFGS単体で、良い結果が出ている。全くダメ(相対L2誤差が、100)なケースもある。残念ながら、ケース･バイ･ケースということだろう。{Adam、L-BFGS、Adam＋L-BFGS}を試行錯誤するしかない。
独マックス･プランク研究所及びノルウェーSimula Research Laboratoryの研究者は、｢PINNs学習において、高精度な近似解が得られる新しいオプティマイザ｣を提案した論文(以下、本論文[*77])を発表した(23年8月15日@arXiv)。結構、重要な論文だと、認識している。
PINNsに用いられている従来のオプティマイザは精度が低い、という問題意識がある。相対L2誤差で計測した場合、勾配降下法(GD)、確率的勾配降下法(SGD)、Adam、BFGS などで、｢10－3程度よりも優れた精度が達成されることはほとんどない｣と指摘している(註：本論文では、BFGSで10－4程度の精度が出ている)。
【1】本論文の主張
本論文で提案する｢エネルギー自然勾配法(E-NGD)｣は、
(1) 精度が高い。相対L2誤差で10－7程度を叩き出す、
(2) 収束までに要する反復回数が少ない。つまり、収束までに要する時間が短い(従来比で数十倍)。
【2】事前整理
(0) リーマン計量
ざっくり言うと、曲がった空間の内積がリーマン計量である。ニューラルネットワークのパラメータがなす｢パラメータ空間｣を考える場合、パラメータ空間をユークリッド空間に限定する必然性はない。というよりむしろ、ユークリッド空間に限定すべきではないだろう。ニューラルネットワークの近似能力の高さは、その柔軟性にある。パラメータ空間をユークリッド空間に限定した場合、その柔軟性が棄損するだろう。従って、パラメータ空間としては、曲がった空間(リーマン空間)を考える必要があるだろう。
リーマン空間では、リーマン計量という名の内積を導入することで、距離を導入することができる。距離を導入することで、勾配を議論することができる。☞続きは、(2)を参照。
(1) 弱形式化(変分形式化)
偏微分方程式の弱形式化(あるいは変分形式化)を考えるとき、オイラー･ラグランジュ方程式が偏微分方程式の弱形式化になるような、エネルギー汎関数を考えることが標準的である。弱形式化では(支配方程式である)偏微分方程式を積分するので、微分可能性の要件が一つ下がる。つまり、例えば2階微分可能でなければならない、という要件が1階微分可能で良いに下がる。そのため｢弱い｣と呼ばれる。もちろん、代わりに、積分可能性が要求される(一般的に言えば、p乗可積分であることが要求されるが、通常の偏微分方程式は2階偏微分なので、通常は2乗可積分で十分)。
ただし、本論文は、弱形式の偏微分方程式を解くわけではない。何をするかは、☛【3】(1)を参照。
(2) 自然勾配の概要
[*79]では、まず勾配について思索している：｢関数fが、最も急に変化する方向は、パラメータθ を同じ微小距離dθだけ各方向に動かしたとき、fの変化分 ∇f=f(θ+dθ)ーf(θ)が最大になる方向dθ のことである｣。そして、次に、真っすぐな空間と曲がった空間に対して、dθについて考えている。
真っすぐなユークリッド空間であれば(dθは)、|dθ|2=dθTIdθと書ける。ここで、Iは単位行列(あるいは、ここでは恒等行列と呼んだ方が座りが良い)。一方、曲がったリーマン空間の場合は、|dθ|2=dθTG(θ)dθと書かれる。ここで、Gは(適当な)正定値行列である。この場合のGはG(θ)と表現されている通り、パラメータ(パラメータ空間の各”点”)毎に異なる(正確には、リーマン空間の各点の接空間毎に定義される)。
ここまでのお膳立ての下で、dθを求めると、dθ＝G-1(θ)∇f(θ)≐∇̃f(θ)となる。∇̃＝G-1∇を、自然勾配と呼ぶ。曲がった空間にまで一般化した勾配ということができるであろう。G＝IとするとI-1＝Iなので、自然勾配は、(通常の)勾配と一致する。自然勾配が定義できたので、次は、自然勾配法に移ろう。
dθは、｢関数fが最も急に変化する方向＝∇fが最大になる方向｣であった。そして、最適化によるパラメータ探索法としての勾配降下法では、｢関数fが最も急に変化する方向＝パラメータを更新する方向＝dθ≒Δθ｣でパラメータを探索するのであった。勾配降下法は、｢損失関数L(θ)を用いてパラメータθを、θk+1ステップ目＝θkステップ目－η∇L(θ)と更新する｣という枠組みである。ここで(適用なkに対して)、Δθ＝θk+1ステップ目ーθkステップ目である。またηは、深層学習の文脈で一般に｢学習率｣とよばれるパラメータである。まとめると、自然勾配降下法は、｢θk+1ステップ目＝θkステップ目－η∇̃L(θ)と更新する｣枠組みである、と言える。
(3) 自然勾配法とニュートン法
自然勾配法、ニュートン法そして、本論文で提案されているエネルギー自然勾配法は、下記で述べるように、統一した見方が可能である。
自然勾配は、適用な正定値行列Gの逆行列と勾配を掛け合わせたものであった。このGとして、具体的に何を用いるかを決める必要がある。当然ながら、ニューラルネットワークのパラメータ最適化に自然勾配法を適用するケースを考える。ニューラルネットワークのパラメータ最適化では、天下り的であるが、Gとしてフィッシャー行列を採用する。これは、パラメータが適当な確率分布に従っていると考え、確率分布の構造を考える(距離の概念を導入する)ことで、パラメータ空間の幾何学的構造を(勾配計算に)反映させている、と捉えることができる。実際、カルバック･ライブラ情報量を”距離”と見做して、(2)と同じようにdθを求めると、めでたく(あるいは見事に)フィッシャー行列が現れる[*80]。
一方、ニュートン法は、Gとして損失関数のヘッセ行列を採用したと見做すことができる。ヘッセ行列は正定値とは限らない(ので、問題によって、ニュートン法は収束しないケースがある等、しばしば不具合を起こす)。
本論文で提案されているエネルギー自然勾配法では、Gとして｢エネルギーグラム行列｣を採用している。グラム行列であるため、正定値ではなく半正定値である。そのため、逆行列ではなく、擬似逆行列を使用する。
(4) 自然勾配は、特定条件下で極めて速い
[*79]では、以下のように述べられている：フィッシャー(情報)行列は、最適点ではヘッセ行列に等しい。従って、最適点の近くでは、自然勾配法はニュートン法と同等で、2次の収束をする。しかし自然勾配法の長所は(中略)、学習の過程でプラトーに捉らわれずに、ここ[プラトー]をうまく抜け出せることにある。(中略)ニュートン法やその変種ではこの特性は出せない。
さらに、[*79]では、自然勾配法が、(場合によっては)数千倍も速い理由を、解説している。端的言うと、パラメータ空間が大きく歪んだ構造をしている場合、自然勾配法は速い。大きく歪んでいるとは、パラメータ空間が特異点を含んでいる場合である。パラメータ空間に特異点が発生している状態を、ニューラルネットワークの文脈で表現すると、出力が同じニューロンが複数あるという状態である。自然勾配法は、パラメータ空間の構造を考慮に入れてパラメータ探索を行っている。このため、特異点に捉えられることなく、探索が行えるので、結果として速く収束する、と結論している。
【3】本論文のアイデア
エネルギー自然勾配法は、偏微分方程式の弱形式を考える。その場合、偏微分方程式の求解は、変分エネルギーの最小化問題に帰着する。つまり、変分エネルギーを最小化するような｢方向｣に、パラメータを探索すること、物理学の支配法則に沿っている(ので、うまくいく可能性が高い？)。自然勾配法の枠組みに、焼き直すと、変分エネルギー(正確には、p-ディリクレ問題のp-ディリクレ･エネルギー[*81]で、かつ本論文ではp=2が扱われている)の2階偏導関数で表わされる｢エネルギーグラム行列｣を、【2】(2)の適当な”正定値行列”と見做なす。
算術的に言えば、自然勾配法では、フィッシャー行列の逆行列を計算する必要がある(実際は、逆行列の計算は計算コストが高いので、適応的に求められる)。ニュートン法では、ヘッセ行列の逆行列を計算する必要がある。
エネルギー自然勾配法では、形式上は、エネルギー･グラム行列の逆行列を求める必要がある。しかし、エネルギー･グラム行列はグラム行列なので、半正定値行列である。つまり、逆行列が計算できるとは限らない。そこで、擬似逆行列を使う。ただし、(ややこしいことに)実際には擬似逆行列も計算せずに、最小二乗問題に変換して、この問題を解いている。なぜかというと、並列計算で(マシンパワーを活用して)高速化できるからである。
【4】比較結果
本論文では、ヒルベルト自然勾配法という手法も扱われているが割愛する(従来のオプティマイザと性能差がなく、議論の大勢に影響がないため)。
(1) セットアップ
1⃣ 検証のセットアップ
ポアソン方程式と熱方程式を対象とする。学習は、初期値を変えて10回行われている。以下で示す結果は、中央値を採用している。比較検証するオプティマイザは、勾配降下法(GD)、Adam、BFGSである。比較指標は、相対L2誤差である(グランドトルゥースは、解析解。解析的に解けるケースを選んでいる)。
2⃣ ニューラルネットワークのセットアップ
ハードウェアは、NVIDIA GeForce RTX™ 30シリーズの3080 Laptop GPUが使用された。ソフトウェア･ライブラリはJAXが使用された(微分はJAXの自動微分モジュールで実行された)。活性化関数はtanhを採用(PINNsとして標準的)。ネットワークの深さは2(浅層)、幅は64。
(2) 2次元ポアソン方程式を対象とした場合
(コロケーション点)サイズは、内部が900で境界が120である(サイズ的には小さい。つまり対象系は、複雑形状でもないし、マルチフィジクスでもないことを意味している)。
GDは、8.2×10-3、Adamは1.1×10-3、BFGSは4.4×10-4である。対して、E-NGDは2.4×10-7であった。E-NGDとBFGSの反復回数は500回である一方、GDとAdamは2×105であった。1回の学習に要する時間は、GDで1.8×10-2秒、Adamで3.7×10-2秒、E-NGDで8.6×10-2秒である。BFGSは、1.8秒を要した。
反復回数が異なるので、学習終了(収束する)までに要したトータル時間では、GDで1時間(3600秒＝1.8秒×2回×105×10-2)、Adamで2時間3分20秒(本論文では、1時間6分となっているが計算が合わない)、BFGSで15分(1.8秒×500回＝900秒＝15分)。E-NGDで43秒(8.6秒×500×10-2＝8.6×5＝43)である。
[*79]で、自然勾配法は数千倍速いと書かれている。上記例では、E-NGDの反復回数が、ほぼ103倍少ない(実際は、400倍少ない)。故に、過言とは言えないだろう。
(3) 1次元熱方程式を対象とした場合
GDは1.6×10-2、Adamは1.0×10-3、BFGSは1.4×10-4。E-NGDは、6.3×10-6である。反復回数は、E-NGDが2×103回で、他は2×105回である。1回の学習に要する時間は、GDで2.2×10-2秒、Adamで3.8×10-2秒、E-NGDで8.3×10-2秒である。BFGSは、4.3秒を要した。
トータル時間は、GDで1時間13分20秒(本論文では1時間12分となっているが、これは1回の学習に要する時間の丸めによるものであろう)、Adamで2時間6分40秒(本論文では2時間6分)、E-NGDで2分46秒(本論文では2分48秒)。BFGSで2時間23分20秒(本論文では、35分48秒。これは反復回数500回に相当するので、誤りと思われる)。
本例では、反復回数は102倍少ない((2)と併せて考えると、数百倍という表現が妥当かもしれない)。
(4) 結果のまとめ
まず、E-NGDの精度は高い。また、1回の学習に要する時間では、E-NGDは長いが反復回数が少ないので(つまりは効率的に探索できていると考えられるので)、トータル時間は、短い。
【5】考察
(0) 本論文では、ポアソン方程式と熱方程式という"エントリーモデル"を対象としている。
(1) PINNsのベンチマークを提案した論文[*74]には、相対L2誤差での比較が行われている。相対L2誤差の水準感を知るために、その結果を参考にしよう。確かに、ほとんどのケースで10-2程度である。例外的に、高次元問題では10-4～10-3の精度が出ている。通常の方法であれば、精度は10-4程度が限界という主張は、妥当なのかもしれない📖。それを鑑みると、10-7～10-6という精度は、特筆に値するのだろう。
(2) 自然勾配法は｢統計的な意味でイイ感じ｣のパラメータを探索する方法で、エネルギー自然勾配法(E-NGD)は｢物理的な意味でイイ感じ｣のパラメータを探索する方法といえるだろう。物理法則に従うことを拘束条件とするPINNsとは、いかにもE-NGDとの相性が良さそうに感じる。
(3) ただし、本論文の例でE-NGDが良い結果を残したのは、"たまたま"であろう。たまたまという意味は、E-NGDの効用が限定的であることを必ずしも意味しないが、問題(物理系)によってその性能は変わるはずである。この場合の物理系は、言うまでもなく、境界条件を込みで考える。【2】(2)で言うところの、適当な正定値行列G(θ)の良し悪しは、パラメータ空間の構造－つまり問題としている物理系に依存するはずだから、E-NGDと相性が良い問題と悪い問題はあるだろう。
(4) [*82]では、自然勾配法について、｢自然勾配法はその高い収束性†にもかかわらず、実用上用いられることは極めて稀である。その原因として、深層学習の設定において、損失関数の非凸性により局所解に陥りやすいことや、学習が不安定になる問題が挙げられる｣と書かれている。収束性で言うと、E-NGDは自然勾配法に"劣る"だろう。それは、エネルギーグラム行列が半正定値行列であるため、必ずしも逆行列が存在しないからである(このため、逆行列の代わりに、擬似逆行列が使われる)。ニュートン法のように収束しないケースがあるだろう。
(5) 本論文は、Ⅰ-5でカバーされていないPINNsの欠点＝計算時間が長い、への対応策という位置付けで取り上げた。Ⅰ-5で取り上げていない他の欠点として収束問題があった。結局、収束問題は解決しない。自然勾配法であれば収束しないは避けられるだろうが、遅い。E-NGDは速いが収束しないケースがある(はず)。実際、本論文でも言及されている通り、1次元熱方程式のケースであっても、1回の学習に要する時間の｢最悪値｣は、GDやAdamより10～100倍、BFGSより1000倍遅い。このケースでは、初期値が悪いという理由であるが、一般的には初期値だけで対処はできないだろう。
👉 ディリクレ問題であれば、初期値だけの問題だろうか？
👉 pディリクレ問題では、エネルギー汎関数が顕に示されるから、事前にE-NGDが有効か否かは、判断できるように思える。
(6) いずれにしても、選択肢が増えるというのは良いことであろう。E-NGDが有効な物理系の探索は、興味をそそられる。今後、そのような研究が行われるであろうか。
(7) [*83]によると、自然勾配法は｢収束は早いものの、汎化性能が低くなることが懸念されていた｣とある。パラメータ空間の幾何構造を反映させて(効率良く探索して)いるのであるから、なるほど、汎化性能は低くなる可能性はあるだろう。学習データと"本番"データで、パラメータ空間の幾何構造が異なっていたら、確かに影響を受けそうである。E-NGDを適用するケースでは、学習データも"本番"データも物理系としては同じはずなので、汎化性能は、それほど低下しないのではないだろうか。ただし、課題としては認識しておくべきであろう。
(8) 尚[*83]には、｢自然勾配法は、勾配降下法に比べて、強い正則化を用いる必要がある｣とも書かれている。E-NGDも同じかもしれないが、PINNsでは｢損失関数の重みの調整｣が有効に機能するから、それで同時にカバーできるかもしれない(と思うのはオプティミスト過ぎるだろうか)。
† 収束性とは、大域的収束性である。つまり、いかなる初期値からでも停留点(局所解)に到達することを主張しているに過ぎない(ので、このステートメントに矛盾はない)。また、いかなる初期値とは言っても、到達速度(収束率)の議論はあるので、初期値を選ぶ必要はある。
📖
追加の参照例1❚ 上記[*74]だけではn=1で心許ないので、他の例も示す。有限要素法とPINNsを比較した論文(詳細は、こちらを参照)では、1次元ポアソン方程式で、10-6が得られている。2次元及び3次元ポアソン方程式では、10-2～10-3である。オプティマイザは｢Adam＋L-BFGS｣である。
追加の参照例2❚ [*87]においても、バニラPINNsで｢良くて10-3｣は、概ね成立している(拡散方程式では、10-4も出ている)。ただし対象は、対流方程式で、オプティマイザはAdam及びAdam＋L-BFGS。反応方程式や反応拡散方程式では、良くて10-2、悪いと100である(物理パラメータによって、大きく変動する)。☛[*87]は、アルゴリズムを工夫して、相対L2誤差を低減するという論文。[*87]で提案されている手法では、反応方程式や反応拡散方程式でも、10-3を達成している。[*87]と[*88]に関しては、こちらを参照。
追加の参照例3❚ [*88]では、バニラPINNs✖1次元バーガース方程式で、10-4～10-3である。☛[*88]は、多種多様なケースで高性能な適応(型)サンプリング手法を構築した、という論文。
Ⅱ-1 アメリカンオプションの評価
[*11]は、米ダートマス大学数学科の学部生が書いた論文なので、アメリカンオプションの評価にPINNsを使ってみた、というだけの論文である。作成は2021年である。やってみただけの論文が成立するということは、アメリカンオプションの評価にPINNsを使うことは、ほとんどやられていないことを示すだろうから、紹介してみた、という趣旨である。
(0) 前捌き－PINNsへの期待
数々の批判はあるものの、ヨーロピアンオプションは曲がりなりにも、解析解が存在する。アメリカンオプションには解析解がないため、古くは、エクセルでも簡単に計算できるBAW(Barone-Adesi and Whaley)近似解(1987年)や、エクセルで作るのはやや面倒な格子モデルが使われてきた。PCの性能がアップしてからは、元来は、アメリカンオプションには不向きなモンテカルロ法(最小自乗モンテカルロ法：LSM法、2001年)も使われるようになった。
物理や工学でおなじみの有限差分法や有限要素法も、使われている。ちなみに有限差分法では、よく知られているBrennan-Schwartzアルゴリズムが提案されたのは1977年と古いが、その有効性が数学的に証明されたのは、1996年とのことである[*12]。
単純な(プレーンバニラ)アメリカンオプションの価値評価に、大掛かりな手法を取り入れる必要はない。しかし、行使条件を工夫した新株予約権の評価を行うような場合には、柔軟性を付加できるモンテカルロ･シミュレーションは有効である。しかし実のところ、新株予約権の評価は、エクセルでいうところの循環参照計算が必要であるため、計算負荷は(PCにとっては)相当高い。そういう意味で、PINNsが代替となる余地はあるだろう。
さらに、流行りのXVA(X Valuation Adjustment、Xは総称を意味する)計算を考えると、PINNsが活躍する余地は広がるのかもしれない。XVAについては、(節を改めて)後述する。
(1)セットアップ
❶入力値
資産価格(underlying)＝[0,500]のランダムな整数(PyTorchの乱数生成関数を使用) ※underlyingの上限が500ドル
行使価格＝100ドル
ボラティリティ＝40%
行使期間＝3年(つまり満期T＝3年)
時間t＝[0,300]のランダムな整数×100
時間増分Δt＝5×10-3
リスクフリーレート＝3%
配当＝0(従って、ペイオフはヨーロピアンオプションと等価であり、自由境界問題でもない)
❷モデル
Raissiの論文[*3(再)]によって提示された枠組みを、そのまま使用している。ただし、空間座標xが資産価格Sになっている。なお、論文にはパイソンで作成されたモデルのコードが付記されている。
層とニューロンの数＝10層、50ニューロン
活性化関数＝ReLU
重みの初期化＝He初期化(PytorchでKaimingを使用)
損失関数＝支配方程式による損失×100＋境界条件(S=0及びt＝T)による損失。平均自乗誤差で計算。支配方程式による損失を100倍オーバーウェイしている。
正則化＝なし
境界条件＝ディリクレ条件
最適化＝Adam
学習＝8×10-3の学習率で710エポック＋1×10-3の学習率で4700エポック
❸ベンチマーク
Brennan-Schwartzアルゴリズムの計算値を真値として、比較している。
なお、ニューラルネットワークを使ったオプション評価のレビュー論文[*13]を見ると、BAW法やLSM法がベンチマークとなっている。
(2)結果と考察
❶行使価格に近い資産価格で、モデルはオプション価値を正確に予測することが困難。
❷満期から遠い時刻の予測は、近い時刻の予測ほど正確でない。moneynessを予測することで、満期間近の予測を改善できるかもしれない。
❸層の数を10から20、30と増やしても、精度は大きく向上しない。
❹ニューロン数を増やしても、精度は大きく向上しない。
❺インスタンス数を増やしても、精度は大きく向上しない。
❻最大学習率は1.25×10-2から1.5×10-2の間である。学習率を10-2まで上げると、学習を実行できなくなることがある。
❼確率的勾配降下法を用いてモデルを学習させようとしたが不可能だった。これはTensorFlowを使ったRaissiらのアルゴリズム[*3(再)]でも、同様だった。
Ⅱ-2 逆問題：ボラティリティ･サーフェイスの予測
【1】前捌き
韓国の研究者が、PINNsを使ってインプライド･ボラティリティ･サーフェイスを予測した、という論文[*14](以下、本論文)を arXivで公開した(2023年1月2日)。正確には、インプライド･ボラティリティの予測タスクにおいて、PINNsベースの手法と他の手法を比較し、PINNsベースの手法が優秀だったという内容である。PINNsベースの手法は2つ。他の手法とは、深層学習ベースの手法×3、統計的手法×2である。
ステレオタイプな発想をすれば、PINNsは逆問題が"得意"だから、インプライド･ボラティリティの算出に適用するのは相応しい。ところが、実は、そうでもないという内容である。
【2】セットアップ
(1) データ
(インプライド)ボラティリティ･サーフェイスを予測するための、元データとして、2004/1/5～2021/8/13までの日次S&P500指数ヨーロピアンコールオプションを使用。具体的には、学習セットは、2004/1/5～2019/12/31。テストセットは、2020/1/1～2021/8/13。学習セットのうち、新しい日付の20％のデータを、検証セットとして選択している。
各日のボラティリティ･サーフェス･データを作成する際には、マネネス(moneyness)および満期日のペアごとに、取引されたオプション価格の対応するボラティリティ値を使用したマネネスとは、行使価格を原資産価格で除した値である(m=K/S)。
学習で用いるボラティリティ･サーフェスは、SciPy補間パッケージを用いて3次スプライン補間を行い生成している(3次スプライン曲線は、曲げエネルギーが最小であり、物理的に最も自然な曲線を描く)。スプライン補間する理由は、日次のボラティリティの生データは、同じマネネスと満期について格子状のフォーマットで提供されていないからである。
マネネスの範囲は0.9～1.1(つまりイン･ザ･マネー、アウト･オブ･ザ･マネーともに、ディープではない)。権利行使が可能な期間つまり満期の範囲は、0～1(年)である。それぞれの範囲を20分割して、ボラティリティをサンプリングしている。つまり、1日あたり20×20のマトリックスとして与えられる合計400個のボラティリティ値を使用している。｢最も流動性が高く、取引量の多いオプションの特性に基づいて｣このように選択している。
(2) モデル
モデルは、以下の7つである。なお、モデル❸～❼はベストなセッティング、比較されている。また、❸～❺及び❼は、過去10 日間のデータを用いて11 日目のデータを予測するように設定されている。学習率のスケジューリングは、エポック間の検証損失(MAPEで計算)の減少に基づく、学習率スケジューリングを採用している。
◆統計的手法
❶ベクトル自己回帰(VAR)、❷自己回帰和分移動平均モデル(ARIMA)
◆深層学習ベースの手法。一般的な概要説明を付記する(本論文での具体的なセットアップは割愛)。
❸Conv_LSTM(Convolutional Long-Short Term Memory)
Conv_LSTMは時空間予測のための再帰型ニューラルネットワーク(RNN)版で、入力から状態への遷移に畳み込み構造を持っている。Conv_LSTMのアーキテクチャは、LSTMにおける行列の乗算を畳み込み演算に置き換えている。Conv_LSTMは、グリッド内の特定のセルの将来の状態を、そのローカルな隣の入力と過去の状態に基づいて決定する。このアーキテクチャは、多くの時空間予測タスクに適用されている。
本論文でのセッティングは、エポック数が100。バッチサイズ32。初期学習率0.1%。LSTMカーネルサイズが、3×3。隠れ層64層。
❹SA-Conv_LSTM(Self Attention Conv_LSTM)
SA-ConvLSTM は Conv_LSTM の改良版であり、自己注意メモリモジュールの実装が追加されている。自己注意メモリモジュールは、(元々の出自では)画素間の大域的空間依存性を反映するために、各セルの隠れ特徴マップとメモリユニットに対して、自己注意演算を行う。得られた2つの特徴マップは集約される。出力される隠れマップは、この特徴量情報を直接利用する。メモリユニットは、ゲート機構に基づいて更新され、セル間の時間依存性を保持する。他のRNNベースのモデルと比較して、SA-ConvLSTMは、時空間依存性が顕著に改善されている(と言われている)。
本論文でのセッティングは、エポック数が100。バッチサイズ32。初期学習率0.1%。LSTMカーネルサイズが、3×3。隠れ層64層。
❺Conv_TF(Convolutional Transformer)
Conv_TFは、時空間予測を行うために、トランスフォーマー･アーキテクチャを変更し、畳み込み自己注意メカニズムを導入することで、時空間予測を行う。畳み込み自己注意機構は、SA-ConvLSTMと同様の自己注意操作を行うが、シーケンス内の全てのテンソル及び、生成されたクエリ毎のテンソルに適用する。オリジナルのトランスフォーマーと比較した他の顕著な違いは、特徴埋め込み、位置のエンコーディングの変更、および合成フィードフォワードネットワーク(SFFN)の使用である。
本論文でのセッティングは、エポック数が100。バッチサイズ16。初期学習率0.1%。隠れ層32層。アテンション･ヘッドは4。
◆PINNsベースの手法
❻プレーンなPINNs
プレーンなPINNsモデルは、コールオプション価格予測用ネットワークCDNN(S, τ, m, r)と、ボラティリティ関数用ネットワークσDNN(S, τ, m, r)の2つの深層ネットワーク(DNN)から構成されている。ニューラルネットの入力は、原資産価格S、満期までの期間τ、マネネスm、無リスク金利rである。空間座標はSに、mとrが加わっている。なお本論文では、CDNN(S, τ, m, r)で予測されたコールオプションの価格も、CDNNと表記されている。それは、インプライド･ボラティリティでも同様である。損失関数LPINNsは、｜σDNNーσbase|＋｜ブラックショールズ式(CDNN、σDNN)|である。ここで、σbaseは、補間されたボラティリティ･サーフェスから抽出された値である。LPINNsは、逆問題の損失関数＋順問題の損失関数、と解釈することができるだろう。
CとσのDNNは、10,000ノードの1つの隠れ層とソフトプラス活性化関数から構成されている。自動微分は、 PyTorchが提供するautogradを用いている。ノルムはL1ノルムのようである。オプティマイザーの記述はない。
本論文でのセッティングは、エポック数が2000(転移学習によって、それぞれ1,000エポックを2サイクル学習させた)。バッチサイズ256。初期学習率10%。
❼PI-Conv_TF(Pythics-Informed Convolutional Transformer)
PI-Conv_TFでインプライド･ボラティリティを予測し、その値をLPINNsに代入して、さらに学習を進めるというアーキテクチャ。損失関数は、プレーンなPINNsと同様に、ボラティリティのパートとブラックショールズ方程式のパートで構成されるが、後者にはハイパーパラメータが掛かっている。
本論文でのセッティングは、エポック数が100。バッチサイズ16。初期学習率0.1%。隠れ層32層。アテンション･ヘッドは4。また、ハイパーパラメータを0.1としている。
【3】比較結果
(1) 比較する指標
予測されたインプライド･ボラティリティを平均絶対パーセント誤差(MAPE)を用いて、比較する。VARについては、各日の入力ベクトルを、ボラティリティ･サーフェスのポイント毎の値を用いて作成している。ARIMA については、ボラティリティ・サーフェス上の同一点における過去のデータを入力とし、ボラティリティ・サーフェスの点順の値を予測している。
また、予測されたインプライド･ボラティリティをブラック･ショールズモデルに代入して得られるコールオプション価格も、平均絶対パーセント誤差(MAPE)を用いて、比較する。
(2) 平時の比較
❻プレーンなPINNsでの予測は、散々な結果である。❷ARIMAを除く全てのモデルに劣後している。ボラティリティに関しては、MAPE16.48%。一方、❸～❺は、5.79～8.42％である。PI-Conv_TFになると、4.92%と劇的に改善する。
定性的に評価すると、❶VAR、❷及び❻は、サーフェスの変化を捉えることができない。❸Conv_LSTMと❹SA-Conv_LSTMモデルは、境界での予測値が一様ではない。❺Conv_TFモデルは、堅牢な性能を示すが、アット・ザ・マネー付近のように、ボラティリティがマネネスに対して急変する場合、誤った予測をすることがある。
それに対して、❼PI-Conv_TFモデルは、ブラック・ショールズ方程式を、PINNsの損失関数として組み込んでいるため、ボラティリティ･サーフェイスの急激な変化を正確に予測し、この問題を軽減することができる。
(3) 有事の比較
本論文では、｢サブプライムローン危機とコロナ禍初期｣という歴史的な変動期において、類似のデータを学習していないときに、予測を行う際にうまく機能するかどうか？も検証している。
サブプライムローン危機に対して、学習期間を2004/1/5～2008/9/25。テスト期間を2008/9/26～2009/5/11としている。コロナ禍初期は、学習期間を2019/5/12～2020/3/4(本論文では、2009/5/12～2020/3/4となっているが、明らかに誤植である)。テスト期間を2020/3/5～2020/4/21としている。
その結果、平時で優秀であった❼は、サブプライムでMAPE65.23%(最下位❺66.51%についで悪い)、コロナで22.66%と、冴えない結果となった。サブプライムで最優秀は❸の20.37%、コロナでも最優秀は❸で16.94%である。この理由を「トランスフォーマー･アーキテクチャが、LSTMアーキテクチャよりも本質的に複雑であるため、variantデータを与えないと、過学習のリスクが高くなるためと思われる」と分析している。
【4】考察
(1)プレーンなPINNsは、かませ犬として使われている。様々な分野で揉まれている深層学習ベースの手法は強く、PINNsは未だこなれていない、という印象。そうは言いながらも、(損失関数を精緻化するとか、正則化を工夫するなどで)もう少しチューンすれば、PINNs自体のパフォーマンスは上がるのでは？とも思う。本論文ではそうせず、畳み込み＋注意機構という深層学習で実績のある手法と組み合わせている。
PINNsを使う意味があるのかと思うが、PINNsのおかげで、❼はアット・ザ・マネーのコールオプションの価格予測で優れている、とされている。ちなみに、この結果は、Ⅱ-1での結果と矛盾している。
(2)スイス連邦工科大学チューリッヒ校(ETHZ)数学専攻(応用数学)の修士論文[*15](21年6月)は、インプライド･ボラティリティの算出[*16]を含めて、デリバティブ周りにPINNsを適用している[*17]。この修論によれば、1次元(原資産1個)の場合、PINNsは、インプライド･ボラティリティを精度よく算出できる。しかし2次元以上(高次元)だと、精度が著しく落ちる。その理由を、｢高次元でインプライド･ボラティリティを求める逆問題は、ill-posedであるから｣と説明している。いずれにしても、PINNsは金融に限らず、こなれていないのだろう。
(3)インプライド･ボラティリティの算出は、PINNsにとって得意ではないのかもしれない。Ⅰ(7)❶｢観測データを情報の点源と見なし、情報を外部に伝播する前にそれらの点で残差を最小化する必要がある逆問題については、難がある｣に該当するだろうか？
掲題について、まず、量子化PINNsと(古典)PINNsを比較した論文[*19](以下、論文1)を紹介する。次に、量子化PINNsの精度に影響を与える、"要素"を調査した論文[*22](以下、論文2)を紹介する。
なお、quantum-based PINNsという英語を量子化PINNsと訳している。
Ⅲ-1 流体解析で、量子化PINNsがPINNsより良い結果を出した、という論文
Ⅲ-2 古典PINNsとの比較で、量子化PINNsを評価した論文
Ⅲ-3 Qase1 PINNsは、量子最適制御でも有効と主張する論文
スイスの量子技術スタートアップTerra Quantum(とスペシャリティ･ケミカルを製造販売する独エボニック･インダストリーズ)による論文1[*19](23年4月21日に、arXivにて公開)は、シミュレーション分野で、量子化PINNsが古典PINNsより優れた結果を出した、と主張する。該当シミュレーションは、流体解析シミュレーションである。
量子化PINNsは、NISQを前提とする量子古典ハイブリッドPINNsである(損失関数の定式化において、もっとカッチりした量子化PINNsもある)。量子化PINNsは、変分量子回路(パラメータ付き量子回路、アンザッツ)を古典PINNsに追加して、構築されている。
【0.5】参考
西マルチバース･コンピューティング、英オックスフォード･クォンタム･サーキッツ、米ムーディーズ・アナリティクスは、QPINNs(ママ。本稿でいうところの、量子古典ハイブリッドPINNs)を使った洪水予測モデルを開発すると発表した(23年10月31日)[*54]。洪水リスク評価と管理を、より正確かつ効率的にするために量子コンピューティングを利用する。代理(サロゲート)モデルの使用によって洪水予測を高速化する前提で、PINNsの使用→量子化によりさらに高速化、という発想。計算対象は、2次元浅水流方程式。発注者は、英｢環境･食糧･農村地域省｣。
3社は、イノベートUKが支援する量子技術促進基金(Quantum Catalyst Fund)コンペティション第1フェーズへの参加権を獲得した。第1フェーズの期間は3ヶ月(23年11月30日まで)。第2フェーズは最長15ヶ月で、2024年1月に開始。第2フェーズの承認は、第1フェーズの成功に基づく。第1フェーズでは、£0.1milの資金提供を受けた。第2フェーズに移行すれば、予算は£1.2mil。
【1】論文1の主張
論文1が、流体解析を対象に量子化PINNsと古典PINNsを比較した理由は、以下の通り：PINNs が流体解析で正確な解を提供するには、ニューラルネットワークが高い表現力 (多種多様な、場合によっては複雑な問題の解を学習する能力) を備えていることが重要である。高い表現力は、量子計算の強みである。論文1では、3次元Y 字型ミキサーに対して、3次元非圧縮性流体の定常流を仮定して、ナヴィエ･ストークス(NS)方程式(と連続の式)を、量子化PINNs及びPINNsで解いている。
量子化PINNsが古典PINNsより優れた結果を出したと主張するが、優れているとは、｢同程度の精度に到達する(収束する)までの時間が短い(エポック数が小さい)｣という意味である。
【2】PINNs及び量子化PINNs
(1) 損失関数
論文1の量子化PINNsはハイブリッドであり、パラメータ更新(最適化計算)は、(アーキテクチャとしては)古典コンピュータを使って勾配降下法にて行われる(実際は、全ての計算が量子シミュレータで実行された)。損失関数は、量子化PINNsと古典PINNsで同じであり、最小自乗誤差(MSE)にて、計算される(コロケーションポイントにわたって平均化する)。
全体の損失関数は、㊀＋㊁である：㊀PINNsにより得られる解は近似解であるが故に、右辺＝0と置いた支配方程式(連続の式とNS方程式)がゼロにならない、という損失関数。㊁初期条件･境界条件を与える点と計算を実行する点(コロケーションポイント)が異なることにより、 右辺＝0と置いた支配方程式がゼロならない、という損失関数。
(2) ニューラルネットワークの構成
古典PINNsの最初の層は(問題が 3次元であるため)3 つのニューロンで構成され 、次に n = 64 個のニューロンを含む5 個の隠れ層がある。最後から 2 番目の層には 16 個のニューロンがある。最後の全結合層で、3 つの成分を持つ速度ベクトルと、1つの成分を持つ圧力を出力する。活性化関数は、シグモイド関数を採用している。
量子化PINNsは、古典PINNs、変分量子回路(VQC)、及び最後の全結合層で構成されている。VOCは符号化層、変分層、測定層に分けられる。量子化PINNsでは、予備的なデータ処理は古典コンピューターで行われる。このデータは、符号化層の量子ゲートのパラメータに符号化される。変分層で学習が行われ、測定層で量子ビットが測定される。最終的に、出力として古典ビットのセットが生成される。
【3】セットアップ、学習及びシミュレーション結果
(0) 共通
すべての学習ハイパーパラメーター(学習率、スケジューラー･パラメーター、バッチ サイズ) は、量子モデルと古典モデルの間で共有されている。
境界条件は、すべりなし条件を設定。流入口は、速度プロファイルを固定(放物線状の速度プロファイル)。流出口は、圧力を固定。またパラメータ(物性値)は、密度 ρ = 1.0 kg/m3 及び動粘度 ν = 1.0 m2/sを与えている。
ちなみに、密度は、ほぼ空気(1.2kg/m3)と一緒。動粘度は、20℃における水の100万倍。粘度に直すと水の1000倍程度であり、イメージとしては、トニック･シャンプーといったところであるが、現実的に想像し難い流体である。
(1) 古典PINNs
量子化PINNsと比較するための学習は、以下のように実行された：Adamオプティマイザーを使用したミニバッチ勾配降下法を使用して、100 エポックの間、学習が行われた。オプティマイザーの1ステップが、1エポックに相当する。
なお、古典PINNsでは、転移学習の有効性を確認するための1100エポックの学習も行われているが割愛する[*20](ちなみに、1100 エポックの学習に要した時間は、単一の NVIDIA A100 GPUで5 分)。
(2) 量子化PINNs
古典PINNsと同様に、Adamオプティマイザーを使用したミニバッチ勾配降下法を使用して、100 エポックの間、学習が行われた。なお、すべての計算は、量子シミュレータ(QMware)で行われた。なお、量子コンピュータ（NISQマシン)を使った量子化PINNsも実行されたが、その結果は量子シミュレータの結果とは大きく異なったようである。使用したNISQマシンは、英オックスフォード･クオンタム･サーキッツのLucyと米リゲッティ･コンピューティングの AspenM-3(どちらも超伝導方式)。
導関数の計算(⋍量子回路層の逆伝播)には、随伴微分法(adjoint differentiation)を使用している。この方法は、量子シミュレーターで、変分量子回路の導関数を効率的に計算するのに役立つ、そうである。
(3) 結果
量子化PINNsの結果は、ゼロとなって欲しい損失関数が(100エポックの学習後)、0.03程度になっている。それに対して、古典PINNsの結果は、21%悪かったとしている(つまり、量子化PINNsの結果を0.03とすれば、0.036)。
古典PINNsの1100エポックの学習結果を見た限り、同程度の精度に達するまでは、ざっくり750エポック程度を要すると思われる。大きく改善されていると言うことが可能であろう。ちなみに、先の1100エポックの学習＝5分をあてはめると、750エポックの学習時間は3分24秒程度である。
【4】考察など
1⃣ 表現力
論文1も、量子(化PINNs)＞古典(PINNs)の根拠を、量子回路による量子計算の表現力にもとめている。これはコンセンサスを得ていると思われるものの、｢アンザッツの表現力が増しても、パラメータの探索能力は向上しない｣という論文[*21](23年6月20日)もある(こちらを参照)。量子機械学習は、何が出てくるか、まだまだわからない。量子＞古典の根拠を量子回路の表現力に求めることは妥当である一方、量子回路の表現力を(回路設計の工夫等で)継続的に高めさえすれば、量子>>古典になるわけではないのだろう。
2⃣ 適正性
古典PINNs及び量子化PINNsは、連続の式(つまり、質量保存則)に違反した結果を出力している。論文1では、この理由として、古典PINNsに対しては、❶全体損失関数における各損失関数(連続の式、NS方程式、初期条件･境界条件)のスケール間の不一致、❷基礎となるニューラル ネットワークの表現力の制限、❸ポイント･サンプリング戦略、または❹オプティマイザーの選択、を上げている。量子化PINNsに対しては(加えて)、データ符号化戦略(!?)の影響かもしれない、としている。物理法則に違反しないという制約条件･拘束条件を課すことによって、より良い結果を生み出すはずのPINNsにおいて、こういう結果を対象として良いのかは、疑問が残る。
なお、下記Ⅲ-2では、量子化PINNsのオプティマイザーはAdamよりSGDが良い、と言う結果が示されている。
3⃣ 高速性
論文1の量子化PINNsは、古典PINNsに量子回路を追加して構成している。このため(同じエポック数なら)学習時間は、古典PINNs以下にはならない。100エポック･量子化PINNs対750エポック･古典PINNsの差は、如何ほどだろうか?･･･実際は、あまり意味はないにしても、興味が湧く。
4⃣ 量子誤りの深刻性
残念ながら、量子シミュレータ≠NISQマシンは深刻らしい。量子誤り緩和は、どの程度使われているのだろうか。
5⃣ 流体解析
ちなみに、PINNsあるいは量子機械学習という枠をはめずに、量子計算×流体解析を行っているスタートアップを上げると、以下の通り。
❶ Quanscient(フィンランド)・・・有限要素法を使用。流体解析において量子優位性が生じる境界として、必要な量子ビット数62～82、必要なCNOTゲート数578～750と推定。
❷ Classiq technologies(イスラエル)・・・英ロールスロイスと協業して、航空エンジン設計用の流体解析に計算基盤を提供。具体的には、HHLアルゴリズム内の関数に最適化された量子回路を生成。
❸ BosonQ Psi(インド)・・・当面は、量子インスパイアード･アルゴリズムを使って、流体解析(等)を行うようである。
【1】概要
論文2では、量子光方式の連続変数量子コンピューティングを使用することを想定している。量子化PINNs は、(加ザナドゥの)Strawberry Fields で実装された。つまりは、量子シミュレータにて実行された(PINNsを実行するためのライブラリには、Keras/Tensorflowを使用)。評価するために、1 次元ポアソン方程式を使用する。境界条件はディリクレ条件で、ゼロに固定。 ソース項(ポアソン方程式の右辺)として、二次関数と正弦関数を採用する。量子回路パラメータは、ゼロを中心とする正規分布と標準偏差 0.05 で初期化する。
確率的勾配降下法(SGD)オプティマイザーと適応型オプティマイザーの2種類を使って、学習結果を比較する。適応型オプティマイザーには、以下を使用した：RMSprop、Adam、Nadam(Adam with Nesterov motion)、Adadelta、SPSA(同時摂動確率近似)、Limited memory Broyden-Fletcher-Goldfarb-Shanno(L-BFGS-B)。学習結果の比較とは、損失関数の値及び、最終誤差の比較である。損失関数は、｢近似解故に0にならない損失＋境界条件の場所ズレ故に0にならない損失｣で、L1ノルムにて計算される(コロケーションポイントにわたって平均化する)。最終誤差は、L2ノルムで計算される(コロケーションポイントにわたって平均化する)。
二次関数(二次ソース項)の場合、学習率は0.01。正弦関数(正弦波ソース項)の場合、学習率は、0.0001 。確率的勾配降下法(SGD)の場合、学習中は、固定学習率を使用する。オプティマイザーの反復は、500 回実行する。コロケーション･ポイントの数は、バッチサイズと同じで、バッチサイズは32を選択した。
ちなみに、量子PINNs の学習には、約 20分を要した。
【2】オプティマイザーの調査
古典PINNs では、Adam および L-BFGS-B オプティマイザーが広く採用され、成功している。論文2では最初に、量子化PINNsでも同じことが成立するのか、を調査した。結果は、SGDが適応型オプティマイザーよりも優れていた。適応オプティマイザーは、学習ランドスケープ内で終了せずに、極小値に収束する傾向があった。 Adadeltaにいたっては、解へ収束しなかった。
SGDオプティマイザーは、極小値を回避し、最適化環境をより適切に探索できた。SPSAも、SGDと同様の動作を実現できる可能性があるものの、SPSA のパフォーマンスは SGD よりも優れているわけではなかった。二次ソース項では、すべてのオプティマイザーが、比較的うまく機能した。一方、正弦波ソース項では、適応型オプティマイザーは、解の正弦波の性質を回復できず、最終的に重大なエラーを起こした。
古典PINNs では、PINNs の収束を加速するために、Adam後に L-BFGS-B が使用されることが多い。また古典PINNs では、L-BFGS-B は、学習開始時直ぐから使用されることはない(理由は、学習ランドスケープの極小値からエスケープすることなく、すぐに収束するため)。そこで、古典PINNsと同様に、量子化PINNsも、80 回SDGを反復した後に、L-BFGS-Bを使用した。L-BFGS-B は、量子・古典どちらの場合でも、損失関数の値を減らすことができた。対照的に、最終誤差は、SGDと、SGD＋L-BFGS-Bとで、ほぼ同じであった。
追加試験として、マルチステップの SGD＋L-BFGS-Bも実装したが、最終的なパフォーマンスの向上は達成されなかった。
【3】量子ニューラルネットワークの深さの調査
次に、代理ネットワーク内の量子ニューラルネットワークの、ユニット数を変えることにより、量子ニューラルネットワークの深さ (量子回路の深さ)が、計算に及ぼす影響を調査した。
全体として、4 個と 8 個の量子ユニットを持つ量子化PINNs は、同等の最終誤差を持つ一方で、8 つの量子ユニットを備えた”深い”PINNs では、L-BFGS-B が誤った解法で、学習ランドスケープの極小値に収束するため、重大なエラーを引き起こす。論文2では、｢L-BFGS-B への切り替えが早すぎた｣ことを、この原因としている(SGD反復85(ママ)回後に、L-BFGS-B へ切り替えている)。
【4】バッチサイズの調査及び、自動微分と有限差分の比較調査
最後にバッチサイズの影響を調査した。バッチ サイズの増加により、オプティマイザー反復ごとのコロケーション ポイントの数と関連する損失関数が増加する。このため、バッチ サイズを増やすと量子化PINNsが改善されると予想されるが、実際には、バッチ サイズの重大な影響は認められなかった。
導関数導出(微分計算)に対して、PINNsの特徴である自動微分と使った場合と、有限差分を使った場合も比較している。損失関数並びに最終誤差ともに、自動微分を使った場合の結果が優れていた。
【5】まとめ
オプティマイザーの選択が PINNsの精度に最も影響を与えていた。量子化PINNsでは、SGD が、適応型オプティマイザーよりも、正確で安定した解をもたらす。量子ニューラル ネットワークの深さは、PINNs のパフォーマンスに影響する。4 層が 2 層量子ネットワークよりも高い精度と、8 層ネットワークと同等のパフォーマンスを提供することがわかった。バッチ サイズを増やしても、量子化PINNsの精度は改善されなかった。
Ⅳ-0 DeepONetとフーリエ･ニューラル演算子の紹介
【0】本稿の概要
本稿は、代表的な3つのニューラル演算子(Neural Operator)の中で、DeepONet(Deep Operator Network)とフーリエ･ニューラル演算子(FNO)を扱う。グラフ･ニューラル演算子は割愛した。
参考資料は、PINNsやニューラルODEを幅広く含むレビュー論文[*34]、ニューラル演算子のレビュー論文[*46]及び[*48]、DeepONetのほぼオリジナル論文[*47]とFNOのオリジナル論文[*45]、(新しい)ニューラル演算子のレビュー論文[*59]である。
まずDeepONetとFNOについて、それぞれの拡張版を含めて、細かく見ていく。続いて、DeepONetとFNOの制限をまとめ、Physics-Informedニューラル演算子(PINO)を簡単に整理する。最後に、ニューラル演算子の拡張及びPINOの未解決問題で、締めくくる。
【1】はじめに
(1) ニューラル演算子の背景、動機、経緯
1⃣ アカデミアの研究者が書いた[*34]、[*45～*47]には、ニューラル演算子を欲する動機が、明確には書かれていない(から、思考の流れが共有できない)。ニューラル演算子は、現場ニーズとマッチしそうだ･･･ということは推論できる。具体的には、メッシュ非依存と解像度不変がニューラル演算子の特徴であるから、計算効率の向上(計算コストの低減)という現場ニーズとマッチしそうであることは、推論できる。
メッシュ非依存であれば、”使い回し性”が増す。つまり、(初期条件や境界)条件(を与える場所)を変えても、計算対象物が物理的に同じであれば、再シミュレーションする必要はないはずである。
また、(空間)解像度不変であれば、部分的にメッシュを細かくしても(スケールダウンしても)、再シミュレーションする必要はないから、計算効率は上がる(はずである†1)。
一方、営利企業であるNVIDIA(とカリフォルニア工科大学)が書いたレビュー論文[*48]には、しっかりと、Motivationという章立てがあり、次のように記述されている：標準的なニューラルネットワーク(NN)には、科学的モデリングには適さない、根本的な制限がある。学習と予測は、学習データの解像度に限定される。NNは、離散化された学習データしか利用できない一方、多くの科学現象は連続領域で発生する。
[*48]は、｢離散化されたデータ(離散化ベクトル)間の関係を学習することから、連続関数間の関係を学習することが、従来型科学技術計算用深層学習†2パラダイムの制限を突破する指針である｣ことが、容易に想像できるような書き方をしている。この指針に沿った新しいパラダイムが、ニューラル演算子である。重要なことは、NNに対して、次の｢演算子の普遍近似定理｣が成立するので、NNは、連続関数間の関係(演算子)を学習できる、ということである：
『単一の隠れ層を持つニューラル ネットワークは、任意の非線形連続汎関数(関数空間から実数への写像)または、非線形演算子(関数空間から別の関数空間への写像)を正確に近似できる』。
最初のニューラル演算子であるDeepONetは、この｢演算子の普遍近似定理｣に触発されて、開発された。
※ NNが演算子を精確に学習できるという事実は、NNが量子系を高精度に表現できることを数学に的裏付ける。物理屋さんによっては、I knew it!であろうか?
†1 計算コスト低減が期待されたものの、ニューラル演算子は、｢積分演算子｣の評価にコストがかかるため、CNNやRNNに匹敵する効率的な数値アルゴリズムを生み出していない[*45]。この問題意識から、フーリエ･ニューラル演算子が開発された。積分演算子については、【3】フーリエ･ニューラル演算子の項を参照。
†2 畳み込みニューラルネットワーク(CNN)や再帰型NN(RNN)といったプレーンバニラなNN([*34]では、科学技術計算に使った、これらのNNをPgNNと呼んでいる)及び、PINNsやPeNNを指す。
2⃣ ([*48]と同じく、NVIDIAとカリフォルニア工科大学の研究者が書いた)[*59]は、ニューラル演算子に対して、見通しの良い統一的な見方を提示している。具体的には、(すべての)ニューラル演算子を｢線形積分演算子と非線形活性化関数の合成として定式化｣する。そして、ニューラル演算子は(現時点で)、離散化不変性と普遍近似の両方を保証する唯一モデルクラスである、と表明している。また、ニューラル演算子が、(従来型の)ニューラルネットワーク(NN)のアーキテクチャを踏襲している理由として、｢NNが非局所性を組み込むことができる｣という性質をあげている。
[*59]では、ニューラル演算子として、グラフ･ニューラル演算子(GNO)、低ランク･ニューラル演算子(LNO)、多重極GNO†、FNOを取り上げている。そして(先述の通り)、積分演算子(積分カーネル)を如何に近似表現してるか？という統一的な見方を提示している(ちなみに、DeepONetに対しても提示されている)。GNOは、積分カーネルをニストロム(Nyström)近似する。LNOは、積分カーネルを低ランク近似する。多重極GNOは、GNOとLNOの組み合わせとして定式化する。FNOは(趣が異なり･･･)積分カーネルを、畳み込み積分に置き換える。
† Multipole Graph Neural Operator。正式な和訳は、未だないだろうが、古典的なmultipole methodにインスパイアされているそうなので、multipole methodの和訳｢多重極展解法｣から、多重極GNOとした。
(2) PgNN、PINNs、PENN、PeNNとの比較による特徴出し
ニューラル演算子は、連続関数間の関係(演算子)を学習しているという事実から、逆問題への対応が可能で、外挿も可能であることは、[本稿の記述では、数学的なセットアップが甘いものの]自然な帰結 として受け入れられるだろう。以下、ニューラル演算子の特徴出しをした。なお(曖昧な表現ながら)、ニューラル演算子の学習は簡単ではなく、熟練を要する、とされる。
1⃣ メッシュ依存性･･･ニューラル演算子、PENN、PeNNは非依存。PgNN、PINNsは依存。
2⃣ 解像度不変性･･･ニューラル演算子は可。PgNN†5、PINNs†1、PENN、PeNNは不可。
3⃣ 逆問題への対応可否･･･ニューラル演算子†2、PINNsは可。PgNN、PENN、PeNNは不可。
4⃣ 外挿の可否･･･ニューラル演算子†3、PINNsは可。PgNN、PENN、PeNNは不可。
5⃣ 物理法則の遵守･･･ニューラル演算子†4、PgNNは遵守しない。PINNs、PENN、PeNNは遵守する。
†1 不均質材料(heterogeneous material)の｢変位と応力｣を、低解像度で学習した後、高解像度で｢変位と応力｣を再構築したPINNs、physics-informed deep learning based super-resolution framework(PhySRNet)が存在する[*49]。[*59]によれば、PINNsには、解像度不変性を持たせることができる。
†2 教科書的には、PINNsでは可能とされているill-posedな逆問題に、ニューラル演算子が対応できるかは、分からない(難しいと推量)。PINNsとの融合により、可能となるかもしれない。PINNsと融合したニューラル演算子はPINOと表記される
†3 入力データ範囲外のデータに対して、外挿可能。学習中に見られない条件への外挿は不可。PINOでは、可能になると思われる。
†4 PINOでは、物理法則を遵守させることができる。
†5 [*59]によれば、CNNは特殊であり、内挿と組み合わせることで、解像度不変にできる。
(1) 概略的な説明
以下の表記は、[*47]に倣った。DeepONet は 、ブランチ(枝)･ネットとトランク(幹)･ネットという、2 つのサブネットワークから構成される。センサー†xi(i＝1,...,m)における入力関数uを符号化するサブネットワークが、ブランチ･ネット。出力関数G(u)の位置yを符号化するサブネットワークが、トランク・ネットである。Gはuを入力関数とする、非線形演算子G：u→G(u)である。
分かり易く表現すると、ブランチ･ネットは、u(xi)を受け取り、u(xi)の線形和を出力する。この出力を、シンプルに、bkと表記する(ここでk＝1～p)。トランク･ネットは、yを受け取り、yの線形和を出力する。この出力を、シンプルに、tkと表記する。DeepONet全体で言うと、{u(xi)及び、y}を入力として受け取り、非線形演算子Gを学習し、G(u)(y)を出力する。為念：G(u)(y)は、実数値であり、関数ではない。
† DeepONetでは、入力関数に対する有限個の固定入力点xiを、センサーと呼称している。
(2) 突っ込んだ説明
DeepONetの出力G(u)(y)は、この表記よりも、Gu(y)やG(y|u)という表記が、分かり易いかもしれない。改めて、G(u)(y)を数式で表現すると、bkと、tkとの積
G(u)(y)≃∑kbktk＋b0
で表現される(あくまで近似であるから≃としている)。バイアスb0は、汎化誤差を低減させるために、加算されている。
唐突だが、上の式は、変数分離形式となっている。つまりG(u)(y)＝G(y|u)＝(uだけの式)×(yだけの式)という形になっている。ここで、敢えてG(y|u)という表記を挟み込んだ。こうすると、変数分離形式が成立するという仮定は、かなり都合の良い仮定なのでは？との疑念が生じる。ただ、この仮定は成立する。なぜなら、普遍近似定理がこの形式で成立しているからである。成立はするものの、変数分離形式で、出力を構築するDeepONetの表現力は、多くの非線形偏微分方程式に対して、制限される可能性がある[*48]、と認識されている†。
なお、G(u)(y)は、連続量(実数値)yを"引数”とする関数であり、トランクネットの入力もyと表記される。ただし、ニューラルネットワークの連続量を入力することはできないから、適当な数量を入力する。もちろん位置は、xiと無関係であり、ランダムに与えられる。このように、
入力値の個数が異なるため、別々のサブネットワークが用意されている。
† 例えば、不連続性を持つ偏微分方程式で記述される物理系に対しては、DeepONetの精度が低い[*51]。
(3) 補足説明
0⃣ DeepONetは、米ブラウン大学、MITの研究者他によって考案された(2019年)。
1⃣ 離散化ベクトル間の写像を学習する従来のNNに対する数学的舞台装置は、有限次元ベクトル空間で済む。これに対し、関数間の写像を学習するDeepONetに対する舞台装置は、無限次元関数空間である必要がある。具体的に言うと、演算子の普遍近似定理は、バナッハ空間に対して成立している。後に、ルベーグ空間でも成立することが示されている。
DeepONetの入力関数uは、バナッハ空間のコンパクト部分空間で定義される連続関数のクラスに属する。出力関数G(u)は、有限次元ベクトル空間のコンパクト部分空間で定義される連続関数のクラスに属する。
2⃣ ブランチ･ネットワーク並びにトランクネットワークはアーキテクチャを問わない。つまり、畳み込みニューラルネットワーク(CNN)、でも再帰型NN(≃LSTM)でも、グラフNNでも構わない。
3⃣ DeepONetの利点は、次の2つとされる[*34]：
① 小さな汎化誤差→ ｢DeepONets が、なぜ小さな汎化誤差を引き起こす可能性があるのか、理論的にはまだわかっていない[*47]｣。PINNsも小さな汎化誤差を保証する可能性があるが、それは十分に大きなデータセットが必要と指摘されている[*34]。ニューラル演算子は、比較的小さなデータセットで、小さな汎化誤差を保証するとされている。
② 学習データの量に関する、学習誤差及びテスト誤差の急速な収束→ ｢深層学習で、誤差の指数関数的収束が観察されたのはこれが、初めてです[*47]｣。ただし、データセットが大規模になると、指数的収束ではなく、代数的収束になる[*34]。
4⃣ 出力関数G(u)の基底の数pは、少なくとも10個である[*47]。pは、G(u)(y)＝∑k･･･に現れる、添え字kの個数でもある。
5⃣ 学習の種類で言うと、DeepONetは、継続学習に該当する(PINNs等はインスタンス学習)。学習は、事前定義された入力空間内において、オフラインで実行される。このため、新しい条件が入力空間内にある限り、それ以上の学習は必要ない。入力空間外の任意の入力の場合は、再学習が必要となるが、入力空間が十分にサンプリングされていれば、これは比較的軽いと考えられる[*46]。
継続学習において指摘される、破滅的忘却も回避されていると思われる。
6⃣ オリジナルの普遍近似定理は、隠れ層1の浅いネットワークに対して成立している。DeepONetでは、表現力を増すため、より深いネットワークに拡張されている。
7⃣ オリジナルのDeepONetは、入力が固定位置で与えられているが、これが任意の位置(連続領域)に拡張された[*48]。
8⃣ 複数入力のDeepONet[*46]。演算子の普遍近似定理は、バナッハ空間上の単一の入力関数に対して定義されている。 DeepONetを現実的なセットアップ、つまり複数の入力関数と多様なアプリケーションに適応させるために、複数入力の演算子定理が理論的に定式化された。複数入力の演算子は、バナッハ空間のテンソル積として定義される。
なお、サンプルとして、患者固有の情報とともに、心臓の収縮期および拡張期における幾何学形状のグレースケール画像を使用して、胸部大動脈瘤患者における機械生物学的傷害の初期分布と範囲を予測した例が、提示されている。
9⃣ DeepONetは、対象の支配方程式(多くは、偏微分方程式)についての事前の理解(知識)を必要としない。ただし、支配方程式の知識を使う(PINNsと融合する)ことで、パフォーマンスは向上する(【1】(2)を参照)。
(4) デモンストレーション
1⃣ 周囲液体圧力の時間変化に応じた単一気泡の形成[*34]
DeepONetとLSTMを使用して、周囲液体圧力の時間変化に応じた単一気泡の形成を表すデータセットをモデル化した事例。まず、巨視的モデルとしてレイリー･プレセット方程式を、微視的モデルとして散逸粒子動力学を使用して、データセットを生成した。次に、ガウスランダム場を使用して、この動的システムの入力信号として機能する、様々な圧力場を生成した。 結果として、学習データがどれほど疎らであるかに関係なく、DeepONet は、液体の圧力軌跡を予測する際に LSTM を上回るパフォーマンスを発揮した。
さらに、入力が学習範囲内に含まれない場合、つまり圧力場の相関長が、学習範囲外にある場合についても検討された(外挿機能のチェック)。このケースでは、当初は正確な予測を行うことはできなかったが、事前学習された DeepONetトランク ネットワークに学習を転送し、数個のデータを追加することで問題は軽減できた。最後に、計算時間が、48時間(CPU時間)から、数分の 1 秒に短縮できることも実証された。なんと、5桁の改善である。
(5) 拡張
1⃣ 重みの自動適応[*46]
最適化プロセス中、制約(初期条件、境界条件)を満たすために、クエリポイント†の一部に他よりも、多くのペナルティを課すべき場合がある。このような場合、不均一なクエリポイントの重みを適切に設計すると、精度が向上する。これらのペナルティ･パラメータは、学習中に適応的に決定する必要がある。ペナルティ･パラメーターは、ネットワーク･パラメーターと並行して勾配降下法によって更新できる。
† DeepONetでは、評価点をクエリ･ポイントと呼ぶ(らしい)。PINNsでは、コロケーション･ポイントと呼ぶ。
👉 [*59]には、FNOが一番速い、と書かれている。指標はエポック当たりの学習時間(秒)で、比較対象はGNO、多重極GNO、LNO。Nvidia V100 GPUで実行した。ちなみに、FNOは4秒。
(0) 積分演算子、積分カーネル
FNOのアーキテクチャを説明する上で、 積分演算子・積分カーネルの説明は、重要である。[*45]では、標準的なNNのアーキテクチャを｢線形乗算と、非線形活性化関数とを組み合わせることによって、非線形関数を近似する｣と捉え、そのアナロジーとして、ニューラル演算子のアーキテクチャを『(大域的な)線形積分演算子と、(局所的な)非線形活性化関数とを組み合わせることによって、非線形の演算子を近似する』と捉える(なお[*48]でも、同じ立場がとられている)。
標準的なNNでは、離散化ベクトル間の写像を対象とするので、線形乗算で十分であるが、ニューラル演算子では、連続関数間の写像を対象とするので、線形乗算ではなく、積分演算子が必要となる。このアプローチは、演算子の普遍近似定理に準拠したDeepONetのアプローチとは一見、全く異なって見えるが、実際は、DeepONetも同じ枠組みで捉えることができる。
完全に天下り的であるが、この積分演算子を、積分カーネルκ(x,y)と入力関数a(y)による積分で定義する。積分カーネルは、出力領域の任意の点xと、入力領域の任意の点y との間の、学習可能なカーネルを示す。
(1) 概略的な説明
標語的に言うと、FNOは、｢積分カーネルを、フーリエ空間で直接パラメータ化することによって、演算子を定式化し、表現力豊かで効率的なアーキテクチャを実現する｣。具体的に言うと、積分演算子はフーリエ空間で定義された、畳み込み積分演算子で置き換える。そのために、積分カーネルをκ(x,y)＝κ(x-y)†として、畳み込み積分に置き換える理由は、変数分離方式を使いたいからである。畳み込み積分をフーリエ変換すると、うまく分離が行える。フーリエ変換を行う演算子をℱで表し、数式で示すと、
ℱ[κ(x-y)a(y)]＝ℱ[κ(x)]ℱ[a(y)]
のようになる。畳み込み積分が、フーリエ変換を使うと、(フーリエ空間では)積算に置き換わっていることがわかる(実際には、高速化するために、高速フーリエ変換(FFT)を用いる)。言わずもがな、ネットワーク パラメーターは、物理空間ではなくフーリエ空間で学習される。出力関数のフーリエ係数は、データから学習される。
† 積分カーネルに並進対称性を仮定すれば成立＝メッシュ間隔に等周期性を仮定しているから自然に成立？
(2) 図式的な説明
FNOをブロック図的に表現すると、次のようになる：まず㊀入力関数が、局所変換Pによって高次元表現にリフトされる。続いて複数の㊁フーリエ層において、反復的に更新が行われる。㊂最後に更新された入力関数が、別の局所変換Qによって射影され、出力となる。㊀～㊂について、以下にて詳しく述べる。
㊀局所変換Pは、線形変換あるいは、浅い全結合ニューラル･ネットワークによってパラメータ化される変換である。入力関数a(x)が、 P によって、高次元表現P(a(x))＝u0(x) にリフトされる。リフトすることで、ネットワークの表現力を高めている。
㊁フーリエ層は、入力ui(x)が、2系統に分かれて処理される。1系統目は、ui(x)がフーリエ変換、線形変換R、逆フーリエ変換を受ける。2系統目は、ui(x)が局所線形変換Wを受ける。1系統目と2系統目は合成されたあと、活性化関数で活性化されて、ui+1(x)となる。最後のui(x)は、uT(x)である。
㊁フーリエ変換及び、逆フーリエ変換は、FFTによって高速に実行される。FFTはフーリエ級数を使うため、積分カーネルκに、周期的という制約を加える。さらに(フーリエ級数を使うため)、離散化は均一でなければならない。これは、入力関数が、等間隔の格子グリッド上(の点)で定義されていなければならない、という制約をもたらす(補足説明2⃣を参照)。
㊂出力vは、最後のuT(x)に対して、Q(局所変換)を施すことによって得られる。出力は、入力関数を与えた等間隔の格子グリッド上(の点)で、得られる。
(3) 補足説明
0⃣ FNOは、カルテック(カリフォリニア工科大学)(とパデュー大学)の研究者によって、開発された(2021年)。その後も継続的に、カルテックが引っ張っている様子。cf.[*48]、[*59]、[*60]など。
1⃣ 他の手法とは異なり、FNO モデルの誤差は、入力解像度と出力解像度に関係なく一貫している。たとえば、標準的な CNN 手法([*34]ではPgNNということになる)では、解像度が上がるにつれて誤差が大きくなる。
2⃣ 入力関数が等間隔の格子グリッド上で定義されていない場合、FFTは使えない。つまり、FNOの高速性は期待できない(というより、大きな問題には、対応できない)。これを克服するために、㊀不規則な物理グリッドを、FFT が適用される規則的な潜在グリッドに変換する手法[*60](あるいは、(5)4⃣も参照)、㊁グラフ･カーネル、が提案されている。
別のアプローチとしては、時間領域で規則的なグリッドを維持しながら、グリッドが不規則であることが多い、空間領域でカーネル積分を分解するというアプローチがある。例えば、空間領域にグラフ･ニューラル演算子を適用し、時間領域に FNO を適用する。さらに、その他のアプローチとして、フーリエ領域を効率的にパラメータ化するためのテンソル化フーリエ･ニューラル演算子がある[*48]。なお【6】(1)も参照。
3⃣ 入出力が、格子グリッド上で定義されている必要があるFNOは、メッシュフリーか？という疑問が生じる。同じ計算対象物の入力条件を不規則グリッドで再入力するならば、メッシュフリーとは言えないだろう(現実味があるか否かは別として)。FNOの出力は、入力を与えたのと同じ格子グリッド上に出されるが、FNOは解像度不変性を有しているから、メッシュフリーと言って良いだろう。
4⃣ FNOは、❶高周波数モードで、エネルギー減衰が遅い、非線形性の高い偏微分方程式で生じる、複雑な演算子を近似できる、❷グラフベースのニューラル演算子が収束しない乱流領域における、ナヴィエ･ストークス方程式(ファミリー)の、解像度不変解演算子を学習する最初の手法(深層学習フレームワーク)である、❸ゼロショット超解像で乱流をモデル化することに成功した最初の手法である、❹従来のPDE(偏微分方程式)ソルバーと比較して、最大で3桁高速である、❺固定解像度の下で、従来の学習ベースのソルバーと比較して優れた精度を達成する[*45]。
5⃣ FNO は空間的および時間的両方で解像度のダウンスケーリングを実行できる唯一の技術である[*34]。これは、気候モデルのシミュレーションでは、極めて有用な性質と言える。
6⃣ 通常の(プレーンバニラな)FNO では、普遍的な近似特性を保証するために、フーリエ層ごとに異なる学習可能パラメーターが使用される。したがって、ネットワークが深くなるにつれて学習可能なパラメーターの数が増加する。そのため、学習プロセスが困難になり、潜在的に過学習が発生しやすくなる。また、隠れ層 L の数を増やすと、勾配消失により FNO の学習が困難になる[*46]。これらの問題を回避するために、陰的FNOが提案されている(→(5)2⃣を参照)。
7⃣ FNOは、すべての層がパフォーマンスを犠牲にすることなく、同じパラメーターを共有する再帰型(リカレント)ネットワークとして自然に定式化できる反復構造を持っている[*45]((→(5)2⃣を参照)。
8⃣ 粘性係数1.0×10−4 のナヴィエ・ストークス方程式を学習するには、10,000個の学習データを入力(aj)と出力vjの、それぞれで生成する必要がある(つまり2万個)。[粘性係数(原文でviscosityだが,記号はν)は、動粘度(kinematic viscosity)を指しているかもしれない。]
(4) デモンストレーション
1⃣ 解像度を 64×64 に固定した場合でも、既存のすべての深層学習手法を常に上回る。バーガース方程式では 30% 、ダルシー流れでは 60% 、ナヴィエ･ストークス方程式(粘性係数 ν = 1.0×10−4の乱流領域)では 30%、 低いエラー率を達成した。
256 × 256 グリッドでは、ナヴィエ･ストークス方程式を解くために使用される擬似スペクトル法の 2.2 秒と比較して、FNOの推論時間はわずか 0.005 秒(440倍高速)。FNOは、速度に大きな利点があるにもかかわらず、ベイズ逆問題で使用しても精度が低下しない[*45]。
2⃣ 二酸化炭素の回収･貯留(CCS)
最大圧力上昇と二酸化炭素プルームのフットプリントについての厳密な確率的評価を、FNO モデルは、わずか2.8 秒で実行できた。数値シミュレーターでは、これにほぼ 2 年かかった[*48]。つまり、400万倍(6桁!)、高速化されている。
以下は、等間隔の格子グリッド上で、入力関数を定義できないケース[*34]。つまり、FFTが使えず、高速化は期待できない。高速性を回復するために、並列化というアプローチを採っている。具体的には、ドメイン分割に基づいた FNO の並列化バージョンを使用して、地下の不均質貯留層における CO2プルームの一時的な時間発展をシミュレートした事例である：空間解像度60×60 ×64、30 秒の時間解像度で、並列化された FNO モデルが、従来の多孔質媒体ソルバーよりも(GPU を利用せずに)271 倍高速でありながら、同等の精度を達成できることがわかった[*34]。やはり、控えめである。
(5) Physics-Informed以外の拡張
1⃣ U-FNO[*34]
多孔質媒体内の混相流、具体的には広範囲の貯留条件、注入構成、流量、および混相流特性を備えた不均一媒体を通る CO2－水の混相流をシミュレートするために、U-FNOが提案された。
FNOにおけるフーリエ層が{カーネル積分変換＋線形変換}で構成されるのに対して、U-FNOは、{カーネル積分変換＋線形変換＋U-Net}で構成される。U-Net は元々、バイオメディカル画像セグメンテーション用畳み込みニューラルネットワークとして提案されたアーキテクチャである。ちなみに[*45]には、｢2 次元の畳み込みと逆畳み込みを含む、 4 つのブロックで構成される画像間回帰タスクの一般的な選択肢 ｣とある。U-FNOはU-Netを加えることで、高周波成分の学習性能が向上すると主張している([*34]で引用されている[*50])。
U-FNO を FNO および CNNと比較し、U-FNO アーキテクチャが、非常に不均一な地層におけるガス飽和と圧力上昇の両方の予測において、より良いパフォーマンスを提供することが示された。また、U-FNO アーキテクチャは元の FNO の学習精度を向上させるが、複数の離散化での学習とテストの柔軟性を自然に実現できるわけではないことも示された。
2⃣ IFNO[*46]
FNOの安定性を向上させるために陰的FNO(IFNO)が、開発された。IFNO では、各反復層に同じパラメータ セットが使用され、更新が行われる。ここで、学習可能なパラメータは、層に依存しないとみなされているため、層の数に応じて学習可能パラメータの数は増加しない。従って、オリジナル FNO(プレーンバニラなFNO)が遭遇する、過学習問題の主要なボトルネックが軽減される。
IFNO では、非常に深いネットワークを通過する順方向経路は、陰的問題として PDE 解を取得することに似ており、収束する不動点方程式が存在する限り、普遍的な近似能力が保証される。提案されたアーキテクチャは FNO の修正として構築されているため、積分カーネルをフーリエ空間で直接パラメータ化し、高速フーリエ変換 (FFT) を利用して積分演算子を効率的に評価する。したがって、IFNO は、解像度の独立性と効率性に関する FNO の利点を継承しながら、安定性が強化されただけでなく、深いネットワークの限界における精度も向上した。
INFOは、超弾性材料、異方性材料、脆性材料などのケースで性能が実証されている。ガラス･セラミックスに生じた亀裂の一時的な伝播に関して、IFNOとFNOを比較した事例では、変位場の予測において、精度及び計算コストの点で、INFOが優れていた[*34]。
3⃣ dFNO+: ドメインが異なる問題設定の解決策[*46]
この機能は、入力関数のドメインが出力関数のドメインと異なる問題設定のために実装されている。 たとえば、初期条件を解の空間的および時間的時間発展にマッピングしたい場合を考える。このような場合、入力空間は空間ドメイン上で定義されるため、出力を空間ドメインと時間ドメインに同時にマッピングすることは困難である。この目的のために、2 つのアプローチが提案されている。
最初は、追加の時間要素を持つ新しい入力関数を定義するというアプローチ。2 番目は、再帰型NN(RNN)を使用して、出力空間を定義するというアプローチである。2 番目のアプローチでは、解演算子(ニューラル演算子が学習する、関数空間間の非線形演算子G)が、一連の演算子に分解され、G:u(x, t)→ u(x, t + Δt)であるような、時間発展スキームを使って反復的に、各時間ステップの解が得られる。あるいは、領域全体にわたって定義された解に、境界条件をマッピングすることを試みながら、入力空間を出力空間の部分空間として定義することもできる。
4⃣ gFNO+: [複雑形状の解決策][*46]
FNO は高速フーリエ変換(FFT)を採用しており、格子グリッド上での入出力関数の定義が必要である。現実の複雑な形状で定義された問題の場合は、通常、非構造化メッシュが使用される。このため、❶非デカルト領域と、❷非格子メッシュの 2 つの問題に対処する必要がある。
❶ 非デカルト領域で定義された入出力関数に関連する問題を処理するために、境界ボックスを定義し、境界での連続性を維持するために「最近傍」によって入力空間と出力空間を射影する。
❷ 非格子メッシュに関連する問題については、非構造メッシュと格子グリッド･メッシュの間で、補間を実行する。
5⃣ ウェーブレット･ニューラル演算子[*46]
ウェーブレット･ニューラル演算子(WNO)は、周波数と空間の両方で局所化された、ウェーブレット空間内のネットワーク演算子を学習する。画像や信号のパターンをより効果的に学習できる、とされる。具体的には、FNOのフーリエ積分がウェーブレット積分に置き換えられる。 WNO は滑らかな幾何学形状と複雑な幾何学形状の両方を持つ領域を処理できることが示され、解領域と境界に不連続性や急激な変化がある高度に非線形な PDEファミリー(族)の解演算子の学習に適用されている。
【4】DeepONetとFNOとの対比を通じた、制限についての説明[*34]
(1) FNOはDeepONetとは異なり、等間隔メッシュ上の点毎の評価によって、入力関数と出力関数の両方を離散化する。したがって学習後、FNOは入力関数と同じメッシュ内でのみ、解を予測できる。DeepONet は任意の場所で予測を行うことができる。 DeepONetとFNOは、学習に全域のデータを必要とするが、DeepONetはより柔軟である。
(2) DeepONet と FNO は、入力データの次元が大きく、ネットワークの重みの数が大幅に増加するため、通常は、2次元若しくは、小さな3次元問題に限定される。また、学習可能なパラメータの数が増加するにつれて、過学習する傾向があり、学習プロセスがより困難になる。
IFNO はこの課題にある程度対処している。 IFNO では、解演算子(ニューラル演算子が学習する非線形演算子G)は、最初に陰的に定義された写像として定式化され、次に固定点としてモデル化される。後者は、深い層の場合のネットワーク学習の課題を克服することを目的とし、前者は学習可能なパラメータの数とメモリコストを最小限に抑える。
【参考】[*59]で示されているFNOに関する解析結果から、以下のことが言えるだろう。
1⃣ モデルの学習可能パラメータ数を増やすと、FNOは、急速に精度が上がる傾向を示す(対照的にGNOなどは線形で精度が上がる)。なお注意すべき点として、(FNOは一般的な状況で、精度が高いが)パラメータの数が少ない状況では、GNOなどに対して精度が劣後することもある。
2⃣ バーガース方程式、移流方程式、ダルシー流れ、及びナヴィエ･ストークス(NS)方程式に対する解析結果に対して：
㊀ ダルシー流れとNS方程式に対しては、学習データにノイズがあってもなくても、学習エラーに大きな差はない。ただし、バーガース方程式と移流方程式に関しては、差が出る。
㊁ ノイズありデータで学習した場合、ダルシー流れとNS方程式に対しては、テスト･データにノイズがあってもなくても、テスト･エラーに差は、ほぼない。バーガース方程式と移流方程式に関しても、差は大きくない。
㊂ ノイズなしデータで学習した場合の学習エラー≒ノイズなしテスト･データのテスト･エラー(ただし、ダルシー流れ除く)。ノイズありデータで学習した場合の学習エラー≒ノイズありテスト･データのテスト･エラー(ただし、ダルシー流れ除く)。
㊃ FNOは、NS方程式は得意とする反面、バーガース方程式と移流方程式は、苦手とする。
👉 [*59]は、｢FNO は、移流方程式やバーガース方程式では堅牢性が低い｣、｢移流問題は、不連続性があるため、FNO にとっては困難な問題｣とある。
【5】Physics-Informedニューラル演算子(PINO)[*34]
FNO やその他のニューラル演算子は、ゼロショット超解像を実現し、高い周波数まで外挿できることを経験的に示す。ただし、サンプルの解像度が限られているため、モデルは学習中に、より細かいスケールの情報にアクセスできない。
科学モデリングでは、通常、システムを完全に制約する偏微分方程式や、システムを部分的に指定する保存則や対称性などの物理制約にアクセスできる。 Physics-Informedニューラル演算子(PINO)は、データと物理損失の両方を組み込んでおり、より高い周波数へのさらに優れた外挿につながることが期待できる。
(1) PI-DeepONet
疎らなデータ領域の PINNs からインスピレーションを得ることにより、既知の微分方程式を損失関数に組み込みながら、DeepONets を非常に疎なラベル付きデータセットで学習することもできる。 このアプローチにより、Physics-informed DeepONets (PI-DeepONets)が生まれた。
拡散反応、バーガース方程式、移流方程式、アイコナール方程式などのベンチマーク問題に PI-DeepONets を採用し、バニラDeepONet と比較した結果、｢予測精度、汎化性能、データ効率｣が大幅に向上した。さらに、PI-DeepONets は、入出力･学習データのペアなしで、解演算子を学習できるため、従来のソルバーよりも最大 3 桁速く、計算力学における非線形および非平衡プロセスをシミュレートできるようになる。
(2) PI-V-DeepONet
DeepONetをPhysics-informedにして、かつ変分定式化したPI-V-DeepONetが、脆性破壊力学に使用された。PI-V-DeepONet の学習は、変分形式の支配方程式と、いくつかのラベル付きデータを使用して実行された。PI-V-DeepONetは、準脆性材料の脆性破壊における破壊経路、破壊ゾーン、および破壊に沿った損傷を特定した。また、初期の亀裂構成と荷重ステップに対する解決策を、迅速に予測できることも示された。
(3) PI-FNO
PI-FNOは、PINNs の課題 (特にマルチスケール動的システムの汎化と最適化)並びに、FNOの課題(高価で非現実的な大規模な学習データセットの必要性)に対処するために提案された。このモデルは、いくつかのベンチマーク問題 (コルモゴロフ流、蓋駆動キャビティ流など)に展開され、他ソルバーに対して、 FNOの高速性を維持しながら、PINNs及び FNO モデルを上回るパフォーマンスを発揮できることが示された。
【6】ニューラル演算子の拡張及び、PINOの未解決問題[*48]
(1) 線形積分カーネル演算子の非線形への拡張
積分カーネルκ は、非線形、つまり κ(x, y, a) に拡張できることが示唆されている。一例は、トランスフォーマー･アーキテクチャで、不規則グリッドでの予測を可能にするさまざまな試みが行われている。たとえば、演算子トランスフォーマー(OFormer)、メッシュ独立ニューラルオペレーター(MINO)、および 一般ニューラル演算子トランスフォーマー(GNOT)などである。
ただし、注意機構は、入力解像度に応じて2次関数的にスケールされるため、計算的に扱いにくくなる。これを軽減するために、ビジョン･トランスフォーマーでは、パッチを使用して次元を削減している。しかし、固定サイズのパッチを選択すると、モデルが固定解像度に制限され、演算子が生成されない。
演算子学習用のスケーラブルで表現力豊かなトランスフォーマー･アーキテクチャを開発することは、未解決の問題のままである。
(2) 確率的マッピングへの拡張
関数空間での確率的マッピングを近似するニューラル演算子は、敵対的生成ネットワーク(GAN)または、関数空間に拡張された拡散モデルに基づいている。特に、拡散ニューラル演算子は、ガウス･ランダム場を入力として受け取り、指定された任意の解像度でサンプルを生成できるスコア演算子を学習する。これは、火山活動や地震活動、ナヴィエ・ストークス方程式の統計、解像度のない視覚センサーなどの問題に適用されている[56-58]。
(3) 未解決問題
データと物理的制約を組み合わせた(平たく言えば、Physics-Informedと融合した)ハイブリッド型の演算子の学習には、未解決の問題がいくつかある。
❶ PDE 制約からの勾配は、autograd、フーリエ連続法†、有限差分法等の様々な方法で、コストを増加させながら、さまざまな忠実度レベルまで計算できる。 最適なトレードオフを達成するための、マルチ忠実度アプローチの使用は、未解決の問題である。
❷ データが限られた領域や高次元での最適化の課題を克服するには、適応サンプリングなど、PINNs ですでに検討されている方法と統合されたモンテカルロ･サンプリングなどの新しい統合が必要である。さらに、関数空間の不確実性を効率的に定量化する方法の開発は、特に高次元において依然として未解決の問題のままである。
† Fourier continuationの和訳(だが、正式な和訳は、未だないと思われる)。非周期関数のフーリエ展開を精確に実行する方法論らしい[*52]。
解像度不変性を持つデータ駆動型ソルバーとして、ニューラル演算子が知られている。DeepONetとフーリエ･ニューラル演算子(FNO)は、代表的なニューラル演算子である。解像度不変性が有効に機能するユースケースとして、真っ先に思い浮かぶのは、気候モデルであろう。解像度不変であれば、例えば100kmのメッシュで分割した領域に対して解析した結果を、25kmのメッシュで細分化した領域に、解析結果をサイズに適応させた上で適用することができる。気候モデリングの文脈では、これを、ダウンスケーリングと呼ぶ(画像解析の文脈とは逆である)。
(カナダ･ケベック州モントリオールにある)ミラ･ケベックAI研究所[*55]の研究者他は、(力学的)ダウンスケーリングにおいて、より有効なフーリエ･ニューラル演算子(FNO)を提案した論文[*56](以下、本論文)を発表した(2023年5月30日[第2版]@arXiv)。ダウンスケーリングに FNO を使用したのは、世界初と主張している。
【1】本論文の主張
本論文では、独自のFNOモデルである、ダウンスケーリングFNO(DFNO)モデルが提案されている。このDFNOモデルは、低解像度学習データを使っても、高解像度学習データを使ったFNOより、(ゼロショット予測でさえも)精度が高いと主張する。
【2】事前整理
(0) 気候モデルにおけるダウンスケーリング
気候モデルにおけるダウンスケーリングには、力学的ダウンスケーリングと統計的ダウンスケーリングがある。本論文は、力学的ダウンスケーリングに該当する。一般的には、部分的により細かくメッシュを切って、数値解析を行うことが力学的ダウンスケーリングである。詳細な計算が要求されるため、計算コストが高い。計算コスト対策として、ニューラルネットワークの利用が期待され、実際に、行われつつある。
本論文によれば、畳み込みニューラル ネットワーク(CNN)と敵対的生成ネットワーク(GAN)が、ダウンスケーリングのためのアーキテクチャとして人気がある。これまで、降水量、含水量、温度といった気象学的に重要な物理量が、ダウンスケーリングされている。
日本は、統計的ダウンスケーリングの取り組みが、主であるように思える。
(1) 保存則の適用による精度の向上
1⃣ DFNOモデルには、[*57]で提案されている”制約層”が組み込まれてる。その理由は、｢制約層を使用してネットワークを学習すると、誘導バイアスが導入され、ネットワークがより正確なダウンスケーリング結果を得るのに役立つ｣ためである。
ここで言う”制約”とは、保存則の縛りをつけるという意味である。保存則については、2⃣で述べる。ネットワークアーキテクチャ的には、最終出力の後に、制約層が追加される。[*57]では、3つの制約層が提案されているが、本論文では、その中で、ソフトマックス制約層が採用されている。
制約層を組み込まないDFNOモデルの出力を、説明の便宜上、非制約出力と呼ぶ。同様に、制約層を組み込んだDFNOモデルの出力は、制約出力と呼ぶ。ソフトマックス制約層は、低解像度の非制約出力に、高解像度の非制約出力で生成される｢ソフトマックス関数｣をかけることで、制約出力を算出する(非制約出力を制約出力に更新するというイメージ)。ここでいうソフトマックス関数は、もちろん活性化関数で用いられる、ソフトマックス関数である。
2⃣ 保存則について
”制約層”における保存則とは、小難しい物理的保存則ではなく、単純な収支関係である。つまりPhysics-Informedほど大げさではない。具体的に言うと、㊀高解像度の各値が非負であること及び、㊁高解像度の各値の合計が、低解像度の値に等しいこと、を指している。各項補足すると、次のようになる：
㊀→ 気候モデル･気象学に現れる、水分量のような物理的な量は、通常、マイナスにはならないので、非負という制約をかけている。
㊁→ 高解像度の各メッシュにおける値を合計すると、低解像度メッシュの値になる、という至極当然の収支関係を要求しているに過ぎない(が、計算機は”当然”が理解できないので、当たり前のことを制約として課す)。
【3】本論文のアイデア
(1) 哲学的な説明
FNOは、低解像度と高解像度の区別をしない。共通の関数を学習するから、解像度の区別が必要ないという哲学である。理想的には、そうであるが、所詮は近似の世界である。理想的に十分なデータが入手できれば、ニューラルネットワーク(NN)は任意の精度で近似できるが、それは現実的ではない。そこでDFNOは、低解像度と高解像度は別の世界であると割り切っている(と思われる)。ある意味、ダウンスケーリングに特化している。低解像度と高解像度は別の世界と割り切った上で、DFNOは、低解像度で成立している関数を、高解像度で成立している関数に、写像する演算子を学習する。そうすることでダウンスケーリングにおいては、精度が高くなるか？というアイデアである(と理解している)。
(2) アーキテクチャ的な説明
従って、アーキテクチャ的には複雑になる。まず前処理として、低解像度のデータから、低解像度の関数を構築する。ここで使用するのはNNであって、ニューラル演算子ではない。つまり、有限次元離散化ベクトル空間間の写像を学習している。具体的には、残差畳み込みネットワーク(いわゆるResNet)を使用している。
次に、ニューラル演算子は、入力が(連続)関数であるから、離散化ベクトルを(連続)関数に変換する必要がある。本論文では、バイキュービック補間関数を使って、(連続)関数に変換している。単純な関数であるが、2つの理由から、それで十分であると説明されている：①ResNetは、高次元の埋め込みを学習することができるため、単純な補間であっても、低次元のターゲットに対して高い表現力が維持される、②FNOは、入力が単純であっても、高度に非線形な演算子を学習することができる。
バイキュービック補間関数を使って変換された連続関数が、DFNO(この処理に関しては、DFNO＝FNOである。つまりアーキテクチャ的には、FNOと同じ)の入力で、出力が高解像度の関数である。ここで、分かり易さを優先して、一旦、高解像度の関数と書いた。正確には、｢低解像度関数の出力を補間して、高解像度関数の出力とする補間関数｣である。
(3) まとめ
全体で見ると(改めて言うまでもないが)、低解像度入力→高解像度出力になっている。DFNOの学習データセットは、低解像度データセットと高解像度データセットがダブルで必要となる(はず)。ただし一度、高解像度の補間関数を学習すれば、さらに高解像度である｢超解像度｣の予測は、学習なしで可能(いわゆるゼロショット予測が可能)。
【4】DFNOの比較検証
DFNOは、1⃣CNN及びGAN、2⃣FNO、と比較されている。
(0) 制約層の有無
本論文では、制約層のない場合とある場合に分けて、ベースライン(後述)との比較が行われているが、煩雑さを避けるため、本稿では、制約層がある場合に限定する。結果として、制約層ありが、高精度であることが示されているからである。
(1) 学習データセット
2種類用意されている。一つは、⓵ナヴィエ･ストークス(NS)方程式[*58]の数値解である。解像度64(km?)×64で、1万個が準備された。7千個が学習用、2千個が検証用、1千個がテスト用に使われる。平均プーリング層を介して、解像度 32×32 および 16×16の、アップスケーリングが行われる。気候モデルの文脈では、高解像度化＝ダウンスケーリング、低解像度化＝アップスケーリングである(画像解析の文脈では、逆のようである)。
もう一つは、⓶ERA5気候データセットである。これは、欧州中期天気予報センター(ECMWF)が提供している様々な気象データで、本論文では大気の総含水量が選択されている。総含水量＝水蒸気、雲水、雲氷を含む大気中の水分含有量合計の鉛直積分、である。高解像度サンプルの平均プーリングによって32×32 および 64×64 の低解像度サンプルが生成された。また、4万個が学習用に、検証用とテスト用に1万個が準備された。
(2) ベースラインと評価指標
1⃣ CNN、GANとの比較
ベースラインは、単純なバイキュービック補間である。ベンチマークには、CNNとGANが用いられる。つまり、ベースラインからの改善度合を、DFNOとCNN、GANとの間で比較する。ベンチマークとして、この2つが選ばれた理由は、これまでの気候モデルのダウンスケーリングで、使用されてきたからだろう(故に、選択としては妥当であろう)。正確に言うとベンチマークは、32×32の解像度データで学習されたCNN-2、GAN-2。さらに、64×64の解像度データで学習されたCNN-4とGAN-4である。
さらに、例えば、CNN-2の出力は32×32である。CNN-2における16×16や64×64の出力は、2つの方法で作成される。一つは、バイキュービック補間を使って、16×16や64×64の出力が作成される(CNN-4やGAN-2、GAN-4も同様)。この手法を、便宜上、❶｢補間アプローチ｣と呼ぼう。もう一つの方法は、平均プーリングを介して、NNの学習プロセスで、他解像度の出力を生成する方法である。この手法を、便宜上、❷｢プーリング･アプローチ｣と呼ぼう。
評価指標として、平均自乗誤差(MSE)と平均絶対誤差(MAE)、及びピーク信号対雑音比(PSNR)と構造類似性指数測定(SSIM)、が使用された。
2⃣ FNOとの比較
ベースライン(この場合、グランド･トルゥース†)は、NS方程式の数値解。解像度32×32と、64×64の2種類で比較する。DFNOは、16×16と32×32で学習したDFNO-2(出力は32×32)と、16×16と64×64で学習したDFNO-4(出力は64×64)がエントリー。なお、DFNO-2にとって、64×64はゼロショット予測になる。同じく、DFNO-4にとって、32×32はゼロショット予測になる。
FNOは、 解像度32×32で学習したFFNO-32と、解像度64×64で学習したFFNO-64がエントリー。
評価指標は、平均自乗誤差(MSE)と平均絶対誤差(MAE)である。
† 元々、気象学の分野で使われていた用語らしい。
(3) 比較結果
1⃣ CNN、GANとの比較
データセットとして⓵NS方程式の数値解を使った、❶ 補間アプローチ(本論文では、表3)、❷ プーリング･アプローチ(本論文では、表5)を見てみよう。以下、⓵と❶のケースは、本論文と同じく表3と呼ぶ。すなわち⓵と❷なら表5と呼ぶ。ちなみに、表2や表4は、制約層なしの場合である。
細かいことを抜きにして言えば、DFNOがぶっちぎっている。
細かく言うと、以下の通り。表3で、DFNOは一つを除いて、ベストパフォーマンスである。つまりMSE･MAE･PSNR･SSIMの全項目、16×16～64×64の全ケースで、PSNRの16×16を除いて、最善である。PSNRの16×16で最善はCNN-2で、次善がGAN-2。全体的に言うと、GANよりもCNNが優秀である。表5でも、PSNRの16×16を除いて、DFNOが最善。PSNRの16×16で最善は、ベースラインで、DFNOは次善であった。CNN＞GANの傾向は変わらず、特にCNN-4の性能が高い。
次に、データセットとして⓶ERA5気候データセットを使った、❶(本論文では、表7)と❷(本論文では、表9)を見てみよう。表7 でも、DFNOが優秀であるという傾向は同じ。PSNRの16×16を除いて、最善(PSNRの16×16の最善は、CNN-2で次善はGAN-2)。ただ、CNN≃GANで、GAN-2がCNN-4よりも、少し優秀という結果がでている。表9では、PSNRの16×16含めて、DFNOがベスト。蛇足ながら、やはりGAN-2が優秀。
2⃣ FNOとの比較
本論文では表1に比較結果が示されている。簡単に言うと、DFNOの圧勝。MSEだと桁が1つ小さい。具体的に言うと、例えば、DFNO-2の32×32で0.0004。FFNO-32は0.0101、FFNO-64は0.0113という具合である。MAEだと桁は同じであるが、例えばDFNO-4の64×64が0.0168に対して、FFNO-32が0.0788、FFNO-64が0.0739という具合である。
ただし、より本質的には、DFNOのゼロショット予測が、FFNOの予測を上回っていることである(と本論文は主張する)。例えば、DFNO-2にとって64×64はゼロショット予測であるが、MSEは0.0018であり、FFNO-32の0.0136及び、FFNO-64の0.0118を大きく上回っている。MAEでも傾向は同じ。同様に、DFNO-4にとってのゼロショット予測である32×32のMSEは0.0012に対して、(記述の通り)FFNO-32は0.0101、FFNO-64は0.0113である。MAEでも傾向は、同じ。
なお、当然ながら、DFNOにおいて、学習済予測はゼロショット予測よりも、精度が高い(つまり、DFNO-2の32×32は、DFNO-4の32×32よりも精度が高い)。また、同じゼロショットでも、｢外挿に相当する｣64×64に対するゼロショット予測よりも、｢内挿に相当する｣32×32に対するゼロショット予測の精度が高いことは、納得性が高い。つまり、DFNOは、正常にワークしているという納得性が高いと思われる。
【5】考察など
(1) 確かに精度は高いかもしれないが、DFNOは、計算時間が長い(電力消費量が多い)のではないか？と懸念される。推論時間という意味では、DFNOは事前学習しておけば良いのだから、問題とはならないだろう。
(2) PI-FNOと、DFNOを比較したら、どうなるのか興味がある。
新型コロナが記憶に新しいパンデミックの発生や、激甚災害の増加など、社会や自然における極端な出来事は、壊滅的な結果をもたらす可能性がある。そのため、極端現象を(効率的に)予測できれば、メリットは大きいと考えられる。しかしながら、極端現象は滅多に発生せず、一見良性の状態から生じ、複雑でしばしば未知の無限次元システムに属しているため、特徴付けるのは困難である。
米MITと米ブラウン大学の研究者は、極端現象の定量化(発見と予測)に、実験計画法(＝能動学習)を適用するという枠組みを提示した論文[*61](以下、本論文)を発表した(22年9月20日@arXiv)。本論文では、極端現象の例として、次の3つを扱っている：パンデミック、巨大波、船舶設計のための構造疲労。
能動学習を適用するという枠組みであれば、標準的には、ガウス過程が使われるであろうが、本論文では、ニューラル演算子が使われている(ところが面白い)。
【1】本論論文の主張
本論文は、極端現象の定量化において、実験計画法(BED)＋ガウス過程(GP)ではなく、ベイズ的(Bayesian)BED＋ニューラル演算子を適用すると、 いくつもの優れた結果が得られたと主張する。
定量化とは、具体的には、注目する確率的物理量の確率分布関数(PDF)を定量化することを指している。｢優れた｣には、効率的という意味が含まれる。ここで、効率的であるとは、サンプル複雑性(能動学習の文脈では、ラベル複雑性とも呼ばれる)が低いという意味である。もっと平たく言うと、極端現象の発見タスクの学習に必要なサンプル数が、より少ないという意味である。
【2】事前整理[*62]
(1) 言葉の整理
1⃣ 実験計画法と能動学習
予め実験プロトコルを設計しておくことにより、最小限のコストで必要な情報を得る方法論を、｢実験計画法｣と呼ぶ。(逐次的な)実験計画法は、機械学習の文脈では、能動学習(active learning)と呼ばれる。ベイズ的実験計画法(つまり、ベイズ的能動学習と呼んでも同じ意味)は、ベイズ推定の枠組みを使った実験計画法である。具体的には、必要な情報を予測する予測器の確率分布(事前確率分布)を設定し、事後確率分布を学習する｢予測器が確率的に選択される枠組み｣を、 ベイズ的実験計画法(BED)と呼ぶ。
2⃣ ベイズ最適化と能動学習
ややこしいことに、能動学習は、ベイズ最適化を含む広い意味で用いられることがある。また、ベイズ最適化と同一視されることもある。ベイズ最適化は(狭義には)、｢少数の試行から最適な実験設定を探索する技術｣である。能動学習は(狭義には)、｢少数の学習データからより良いモデルを得るための技術｣である。
なお、ベイズ最適化は、GPの採用を前提として語られることが多いが、GPを採用しないベイズ最適化も存在する。
(2) 能動学習における獲得関数
1⃣ 不確実性サンプリング
どのサンプルを選択するかを判断するために使用される獲得関数は、ベイズ最適化においても、重要な要素である。能動学習の多くは、何らかの意味で｢情報量の多い｣サンプルか、｢代表的な｣サンプルを選択するような獲得関数を利用する。
現在の予測器にとって、最も不確かなデータを選ぶという考え方に基づく獲得関数(の選択基準)を、不確実性サンプリング(uncertainty sampling；US)と呼ぶ。｢最も不確かなデータを選ぶ｣そのココロは、｢現状の予測モデルによる予測が、最も不確かなものにラベルを与えればモデルの不確かさが最も減るであろう｣という期待に基づいている。
本論文によれば、US の人気は、以下の3点によるものである：｢実装の容易さ、低コストの評価(小規模データセットの場合)、勾配ベースのオプティマイザーの使用が可能｣。
2⃣ 他のサンプリングに基づく獲得関数
US以外の獲得関数には、予測分布のエントロピーが高いサンプルを選択する｢エントロピーサンプリング｣、一番確率が高いラベルと二番目に確率が高い予測ラベルの確率の差が最大になるようなサンプルを選ぶ｢マージンサンプリング｣、予測確率の最大値が最小のサンプルを選ぶ｢最小確信度サンプリング｣などがある。
(3) 獲得関数の設計方法
1⃣ ハイブリッド･アプローチ
(2)で上げた1⃣及び2⃣は、どちらも｢情報量の多い｣サンプルを選択する獲得関数であった。これに対して｢代表的な｣サンプルを選択する獲得関数として、k近傍距離による重み付けを用いることで、代表的なサンプルを選択するという方法論がある。ただし、代表的なサンプルを選択するアプローチは、情報量の多いサンプルを選択するアプローチと比較して、多数のサンプルが必要となると言われている。このため、両方のアプローチを合わせること(ハイブリッド･アプローチ)で、効率的かつ効果的な獲得関数を構成する方法も研究されている。例として、ミニマックス最適化問題を解くことで、代表的かつ情報量の多いサンプルを選択する方法が、あげられている。
2⃣ メタ能動学習
(ハイブリッドか否かを問わず)獲得関数を人手で選択するのではなく、データから直接学習する、メタ能動学習が提案されている。通常は、強化学習モデルを使って、獲得関数が学習される。
(4) ニューラルネットワークの不確実性を定量化する3つの方法
本論文では、ニューラルネットワークの不確実性を定量化するための技法を3つに分類している：❶単一の決定論的ネットワーク、❷ベイジアン･ニューラルネットワーク推論アプローチ、および❸アンサンブル手法。❸は、❶と❷の中間に位置する。
❶は、決定論的ネットワークの1回の前方経路を使用して、ラベル付けされた出力の平均と分散の両方を学習する。1つのモデルのみを使用することで、学習と評価を安価に行うことができる。ただし、物理的な不安定性/極端な事象を定量化する回帰問題には、適さない。
❷は、ベイズ理論を深層ニューラルネットワーク(DNN)が享受する表現力、スケーラビリティ、予測性能と組み合わせている。ベイズ理論は、このモデルの成功可能性に信頼を与えているが、複雑さという重大な欠点がある。
【3】本論文のオリジナリティ
(0) 前説
本論文の発想は、以下のようなものであろうと推測する：極端現象＝レアイベント＝少ないデータ。少ないデータから予測したい→｢少数の学習データからより良いモデルを得るための技術｣＝能動学習。→ 極端現象の定量化(予測)に能動学習(＝実験的計画法)を適用。極端現象の定量化に能動学習を適用することは珍しいと思うし、故に、面白いと思った。
本論文で提案されているモデルを【2】で整理した知識を使って、細部を表現すると、｢ハイブリッド型の獲得関数を、メタ能動学習(強化学習)で学習するモデル｣となるだろう。強化学習の部分を、GPからニューラル演算子に置き換えた箇所も、もちろんオリジナリティである。こうすることで、極端現象の定量化がより効率的になった等、と主張している。そして、その直感的説明として、｢パラメータ空間において、学習に大きく寄与する領域のみが抽出されることによって、効率的なサンプル選択が行われる(ので、定量化が効率的になる)｣が与えられている。
(1) ベイズ的実験計画法＋ニューラル演算子という座組
本論文では、｢ベイズ的実験計画法(BED)とガウス過程(GP)回帰＝GPを適用した能動学習、による不確実性予測は、うまく機能する｣一方で、｢無限次元システムに適用したり、大規模な学習セットにスケールしたりすることはできない｣と評価する。解決策として、BED＋ニューラル演算子という座組(ニューラル演算子を適用した能動学習)を提案している。
この座組の優位点として、強化学習における探索(exploration)と活用(exploitation)を実行する空間が、適切に分離されていることが、あげられている。GPは、探索と活用とで、同一の｢確率過程のパラメータ空間｣を使う。これに対して、関数空間間の非線形連続演算子を学習するニューラル演算子を使う本論文のフレームワークでは、探索はパラメータ空間で、活用は関数空間で行われる。本論文は、こうすることで、パラメータ空間全体にわたって汎化性能が向上すると主張する。さらに、データ･サイズの3乗にスケールするGPと異なり、ニューラル演算子は、本質的にビッグデータに適している、と述べている。
なお、ニューラル演算子は、特定のタイプに限定されないが、とりあえずDeepONetが想定されている。ちなみに、本論文はMITとブラウン大学の研究者が作成しているが、DeepONetもMITとブラウン大学の研究者が提案した(フーリエ･ニューラル演算子は、カリフォリニア工科大学)。
(2) 実際に採用された獲得関数
本論文では、ハイブリッド型の獲得関数を採用している。極端現象を取り扱っているので、”代表的なサンプル”として、極端現象を選択する。具体的には、分散(→不確実性サンプリング)に、尤度比を掛け合わせた、｢尤度加重不確実性サンプリング｣を獲得関数とする。ここで尤度比は、起こり得る現象と、極端現象の比(重み)である。
(3) 確率分布関数(PDF)を定量化するスキーム
1⃣ 概要
改めて、整理すると、モンテカルロ･サンプリングを使っても、観測データ･セット{xi, yi}に対して、PDFを得ることができる。しかし、これでは効率が悪いので(効率良くPDFを得るために)、写像G：x→G(x)＝yを、代理モデル＝GPやニューラル演算子、で近似することを考えるのであった。完全に天下り的であるが、｢平均モデル｣μ(x)とxのPDFがあれば、yのPDFを得ることができる。
平均モデルμ(x)は、GPだと自然に導出される。ニューラル演算子(ニューラルネットワーク)の場合は、【2】(4)で示したニューラルネットワークの不確実性定量化を使って導出できる。ここでは、結論だけを言うと、本論文では、ランダムな初期重みを使用したアンサンブル･アプローチを採用している(詳しくは、下記2⃣を参照)。
2⃣ アンサンブル･アプローチ
｢少ないデータから(効率的に)良いモデルを得る｣という能動学習の思想･作法を、ニューラルネットワーク･アーキテクチャに埋め込む一つの方法は、ニューラルネットワークの不確実性を利用することである。換言すれば、ニューラルネットワークの不確実性を定量化するための適切な方法を選択することは、能動学習(＝実験的計画法)にとって重要な要素である。
【2】(4)で示した通り、ニューラルネットワークの不確実性を定量化する手法はいくつかある。本論文では、｢不確実性定量化におけるアンサンブル･アプローチの有用性は、依然として議論の的である｣としながらも、結論として、❸アンサンブル･アプローチを採用している。さらに、アンサンブルを作成するアプローチも、ランダム重み初期化、データ･シャッフル、データ拡張、バギング、ブート･ストラップ、スナップショット･アンサンブルなど複数存在する。本論文では、評価精度と分布外検出に関して、❷ベイジアン･ニューラルネットワーク･アプローチと同等以上のパフォーマンスを発揮することが判明した、｢ランダム重み初期化｣を採用している。
(3) 注意事項
本論文の結果は、経験的なものであり、堅牢な保証はなく、実装は慎重かつ計画的に実行する必要がある、と警告されている。また、次のような注意喚起が続く：ニューラル演算子の学習は、簡単ではなく、熟練を必要とする。獲得関数は、予期せぬ病理学的ケースの餌食になる可能性がある。
本論文の付録Hには、一般的な実装上の懸念事項及び解決策が提示されており、参考のため、下記に転記した。
参考：ニューラル演算子と尤度加重サンプリング実施のヒント
㊀ ニューラル演算子(DeepONet)への入力関数は、－1から1の範囲で変化するようにスケーリングされるべきである。
㊁ 注目している物理量(すなわち出力)も、－1から1の間の値にスケーリングされ、正規化されるべきである。
㊂ ユーザーは、ニューラル演算子が実際に出力に適合していることを確認する必要がある。出力を調整する方法として、学習エポック、層、ニューロンの幅を増やすという方法がある。適合していなければ、上記のスケーリングが適切に行われていない可能性がある。
㊃ 病的なケースでは、尤度加重獲得関数の重み＝尤度比は、出力である確率分布関数(PDF)によって、悪影響を受ける可能性がある。出力PDFと、その重みへの影響を監視することが推奨される。
㊄ 同じカーネルでTensorflowの学習と評価を連続して実行すると、速度が低下することがわかった。これに対処するために、学習、評価、および次の獲得サンプルの決定ラウンドが完了した後、Pythonは終了し、次のステップのために新しいインスタンスをロードする。
【4】極端現象を予測する能力を評価した結果
(0) 前説
0⃣ 概要
下記(1)では、ニューラル演算子を使うと、少ないサンプル数でモデルが構築できると述べている。(2)では、GPと直接比較をしており、ニューラル演算子を使うと、色々良いことがあると、述べている。
1⃣ 共通セットアップ
観測データ･セット{xi、yi}のxiは、ラテン超方格サンプリングを使った105個のテスト･サンプルである。ニューラル演算子を使う本論文の座組では、写像：u→yが学習される。ここで、uは、カルーネン･レーベ展開†u=∑ixiφiである。GPを使うケースであれば、写像：x→yが学習される。ここでxiは、あくまで、代理モデルである｢学習済のニューラル演算子｣に対する、サンプル入力データである。
本論文では、誤差指標 としてlog-PDF誤差と呼ぶ指標を使っている。具体的には、近似出力分布の対数(底は10)と真の出力分布の対数(底は10)との差の絶対値をとり、その絶対値を積分した値をlog-PDF誤差としている。真の出力分布とは、ガウシアン･カーネル密度推定(例えば、Pythonのモジュール、scipy.stats.gaussian_kde)を使って得られた出力分布を指している。
† カルーネン･レーベ展開≒主成分分析である。
2⃣ サンプル選択アルゴリズム
新たな入力データが、獲得関数を最大化した場合に、サンプルとして追加する。最大化は、モンテカルロ手法を採用した。勾配降下法を使用するオプティマイザーと比べて、より優れた結果が得られたためである(後述)。その理由として、獲得関数の非凸性を上げている。
(1) 例1：パンデミック
観測データ･セットのyi(注目する確率的物理量)は、感染者割合である。パンデミックを表現する感染症数理モデルとして、SIRモデルを採用している。SIRモデルは、人口集団を感染のステージにより、感受性(susceptible:免疫を持っていないため、今後感染する可能性がある状態)、感染性(infectious)、隔離/回復(removed/recovered：回復とは、感染後に免疫を獲得して、以後は感染しないと目される状態)の3つに分け、感染に係る状態の時間的な変化をボトムアップに記述するモデルである[*63]。SIRにおける感染率(標準的にはβで表される)は、一般にスカラーであるが、本論文では確率変数を持つ関数としている。βは、β(u)である。本論文では、ニューラル演算子としてDeepONetがイメージされているが、その描像を踏襲すると、ブランチネットワークでは、G：β→G(β)が学習されている(※そう明示されているわけではない)。
ニューラル演算子(DeepONet)の出力は、感染者割合であり、感染者割合が従う確率分布関数が、log-PDF誤差に使われる。計算条件は、以下の通り：初期感染者50人、総人口1億人、計算ステップ0.1日、計算期間45日。
結論は、少ないサンプル数で、良いモデルが得られるである。少ないサンプル数、の部分は、ガウシアン･カーネル密度推定(KDE)に要したサンプル数107が、本論文の座組では102に減った、ということのようである(105個から、102個が効率的に選択された)。良い＝精度が良い、の部分は｢わずか 2 つのニューラルネットワーク(今の場合は、ニューラル演算子)のアンサンブルを使用したケースで、反復回数50(正確には47) までに収束し、近似出力PDFとの誤差は 10−3未満になる｣である。
(2) 例2：巨大波(rogue wave；暴れ波とも言われる)
巨大波とは、｢周囲の平均的な波の高さの、2倍以上の波｣が巨大波と定義されている[*64]。観測データ･セットのyi(注目する確率的物理量)は、波の高さの最大値である。本論文では、巨大波のモデルとして、1次元波動乱流について提案されたMMTモデル(分散非線形モデル)を採用している。uは、波の高さ自体である。
1⃣ GPとの比較
初期条件の入力点のサイズ(2、4、6)と、3つの獲得関数(ラテン超方格サンプリング、不確実性サンプリング、本論文の尤度加重不確実性サンプリング)に対して、GPとニューラル演算子を比較している。概ね、｢尤度加重不確実性サンプリングは、他の獲得関数に勝っている｣、｢ニューラル演算子はGPに勝っている｣で良いであろう。優劣は、logｰPDF誤差で測っている。
2⃣ バッチ処理とモンテカルロ法による最適化
ニューラル演算子を使った場合、高次元での計算コスト削減に有利な結果をもたらしている、と述べている。その理由として、サンプルのバッチ処理と、アンサンブル･サイズ＝2の使用の両方が、損失なく実行される、ことを上げている。バッチサイズ25を選択した場合でも、300 サンプルを超えても、フレームワークのパフォーマンスが低下しないことが示されている。
バッチ化されたサンプルは、局所最適化のいくつかの領域から取得されている。最適化手法として、本論文では、モンテカルロ法と勾配降下法によるオプティマイザーL-BFGS-Bとを比較して、モンテカルロ法がより効率的である、と結論している。
3⃣ 二重降下の回避
本論文で採用した、獲得関数＝尤度加重不確実性サンプリングは、二重降下†を回避していることが示された、としている。この理由として、まず｢二重降下は、観察されたシステムのダイナミクスに決定的に寄与するデータのみを選択することによって回避される｣ことを上げ、尤度加重不確実性サンプリングが、それに資するため、二重降下が回避されたと結論している。
† 多くの深層学習モデルでは、一度誤差が減少した後に、再び増加する現象が発生する。その場合、モデルサイズ、データサイズ、学習時間を増やすことによって、精度が再度向上するといったことが生じる。この現象を二重降下と呼ぶ。
4⃣ アンサンブル･サイズN＝2が最も、精度が高いことについて
アンサンブル･サイズN＝2は、N > 2 よりも一貫して優れたパフォーマンスを示している。この理由としては、｢小さな N を使用すると、Thompson サンプリングと同様の方法で、貪欲な探索が課される｣と考えることができる、述べている。
(3) 例3：船舶設計
割愛
【5】まとめ等
(1) ”効率良く”精度の高いモデルが得られるという特徴の他に、次のような優位点が上げられているので、改めて整理する：①メンバー 2つだけの浅いアンサンブルが、最も優れたパフォーマンスを発揮する。②初期データの状態に関係なく(つまり、初期データに極端現象が含まれていなくても)、極端現象が発見される。③二重降下を排除する。④ステップ･バイ･ステップのグローバル最適化と比較して、次善の獲得サンプルのバッチを使用しても、BEDのパフォーマンスは妨げられない。
(2) GP→ニューラル演算子よりも、獲得関数に尤度加重不確実性サンプリングを採用した事が効いているのではないか、と思う。つまり、極端現象の定量化に対して、ガウス過程を用いたベイズ的実験計画法という枠組みでも、かなり有効ではないかと思う。
(3) 巨大波の例●初期条件6点入力のケースで、ラテン超方格サンプリング ＋ニューラル演算子が、相当良い結果を出しているのは、なぜだろうか。
【0】はじめに
大阪大学と神戸大学の研究者が、｢物理現象における保存則を発見するニューラルネットワーク･アルゴリズムを開発した｣とする研究成果を国際会議で発表した(23年5月3日)。成果をまとめた論文[*A-1](以下、本論文)は、23年3月28日に、arXivへ投稿されている。以下、本論文のエッセンスを理解可能な範囲でまとめた。
【1】本論文の訴求ポイント
(1) 保存則を”保存”しながら力学系を学習できるので、データから保存則を発見できる。
(2) ハミルトン系以外の幅広い物理系で、保存則を発見した。
(3) 保存則を”保存”するようにシミュレートすることで、誤差の蓄積を回避できた。
【2】事前整理
(0) 第一積分を保存するニューラル常微分方程式
本論文において、紹介されているニューラルネットワークは、｢第一積分を保存するニューラル常微分方程式(First Integral-preserving Neural Differential Equation:FINDE)｣と名付けられた。ニューラル常微分方程式については、後述する。なおFINDEには、連続時間バージョンと離散時間バージョンがある。
物理で保存量と呼ばれる量を、数学屋さんは、エレガントに｢第一積分(first integral)｣と呼ぶ。保存量を保存するだと却って分かりにくいので、第一積分という文言を採用したのかもしれない。もっとも、可積分系からのアプローチだと、第一積分という文言は、物理屋さんにも馴染みがあるのかもしれない。
(1) 第一積分と、第一積分の保存
方程式の解xに対して、定数であるような特定の関数Fを、第一積分(あるいは保存量)と呼ぶ。つまりdF/dx＝0である。
本論文では、第一積分を保存しながら物理系を学習可能な方法として、射影法と離散勾配法を導入している。射影法は連続時間バージョンに使用され、離散勾配法は離散時間バージョンに使用される。
射影法は、｢誤差等によって保存則が破られる時に、強制的に保存則を満たす状態に戻す方法｣である。離散勾配法は、｢保存則を破るような誤差が生じないように、物理方程式を離散化する方法｣である。
(2) シンボリック回帰
第一積分を発見するアルゴリズムとして、 シンボリック回帰を使用する。シンボリック回帰よりも、関数同定問題という文言の方が、意味するところが伝わりやすい一方で、関数同定問題だと誤解を招きかねない。
一般に、データから関数を推測することをシンボリック回帰と呼ぶ。シンボリック回帰では、一方で加減乗除を用意し、他方で解析関数を用意する。それらをランダムに組み合わせることで、データの入出力関係に適合する関数を見つける(もちろん"関数"は、入力変数のn次多項式を含む)。実際のアルゴリズムとしては、遺伝的アルゴリズムやニューラルネットワークが使用される。本論文では、(sckit-learnライブラリの)gplearnという遺伝アルゴリズムを使って、第一積分を｢発見する｣シンボリック回帰が実施されている。
(3)ニューラル常微分方程式
FINDEは、ニューラル常微分方程式(Neural Ordinary Differential Equation:NODE)の拡張バージョン(NODE3.0)という捉え方ができる。NODEは元々、計算量を減らす等の工夫として、｢微分方程式を解く｣という行程を導入したニューラルネットワークである。本稿では、これを便宜上、NODE1.0と呼ぶ。NODE1.0は転じて、｢微分方程式を学習するネットワーク｣となった。便宜上、NODE2.0と呼ぶ。NODE2.0は、さらに、第一積分を保存するように拡張された。便宜上、NODE2.5と呼ぶ。なお、”保存”の意味をニュアンス的に説明すると、シンプレクティック数値積分法がハミルトニアンを保存する、という場合に近い(英語では、preserveが当てられている)。
(4) 改めて、FINDEについて
FINDEがNODE3.0とみなせる根拠は、その野心にある。(本稿で便宜上、呼ぶところの)NODE2.5では、事前知識として第一積分を与えていた。FINDEでは、データから第一積分を発見するという野心があるため、3.0と考えられる。ここで、FINDEは、なぜ第一積分を保存するように設計されているのか、をPhysics-Informedニューラルネットワーク(PINNs)との比較で表現してみよう。ちなみに、本論文の緒言において、物理系をニューラルネットワークを使ってモデル化する手法として、PINNsとNODEがあげられており、｢今回、我々は NODE に焦点を当てる｣と書かれている。
PINNsは、(右辺＝0という形式にした)支配方程式を損失関数として用いる最適化問題として、定式化された。つまり、支配方程式に従うという物理的な拘束条件を課すことで、ニューラルネットワークの出力が、物理的に意味のない解となることを回避しようとする。
一方、NODE(NODE2.5)は、第一積分に関する事前知識を組み込むことで、物理系(ターゲットシステム)を正確に学習するように拡張された。PINNsとの対比で言い表すと、NODE2.5は、保存則に従うという物理的な拘束条件を課すことで、ニューラルネットワークの出力が、物理的に意味のない解となることを回避しようとする[*A-2]。FINDEはNODE2.5の精神を引き継ぐが、さらに野心的であった。
ニューラルネットワークがターゲットシステムを学習する一般的な状況では、ターゲットシステムの第一積分が未知ということもあるはず。そのような場合、ニューラルネットワークが第一積分を学習できるとは限らない。そのような事態を避けるため、FINDEでは、データから未知の第一積分を見つけて、保存することを目指した。
【3】FINDEの詳細説明
(1) 定式化
FINDEは、微分方程式を学習するニューラルネットワーク(本稿の表現ではNODE2.0、本論文では、｢基本モデル｣と呼ばれている)に加えて、未知の第一積分を出力するニューラルネットワークとして定式化される。基本モデルには、NODE(NODE2.0)やHNN(ハミルトニアン･ニューラルネットワーク、本稿の表現ではNODE2.5)等々が用いられる。基本モデルが第一積分を保存するように設計されている場合(例えば、エネルギー(ハミルトニアン)を保存するHNNの場合)、FINDEはエネルギー以外の第一積分を発見する(そして、ニューラルネットワークに組み込む)。
損失関数で表現すると、NODE2.0あるいはNODE2.5の損失関数に｢保存量である第一積分の微分は0｣という損失関数が加えられる。
(2) 連続時間バージョンのFINDE(continuous FINDE)
先述の通りcFINDE は、FINDEの連続時間バージョンであり、第一積分を保存するために射影法が適用されるのであった。射影法では、第一積分が変化しない方向の空間にダイナミクスを射影する。このようなメカニズムによって第一積分を保存しながら力学系を学習することを可能にし、その結果、データから未知の第一積分を見つけることができる。
基本モデルと比較してcFINDEでは、未知の第一積分を出力するニューラル ネットワークの追加計算、複数の行列乗算、および逆演算が上乗せされる。しかし本論文は、｢第一積分の数をKとした場合の逆演算の計算コストはO(K3)である。K が小さければ、計算量コストは、それほど高くない｣、｢先行モデル(ラグランジュ･ニューラル ネットワーク、ニューラル シンプレクティック形式、及びCHNNのような多くの先行モデルも同様に、逆演算を必要とする｣ため、大きな問題ではないと結論している。
(3) 離散時間バージョンのFINDE(discrete FINDE)
cFINDE は、第一積分で時間離散化誤差が発生する。dFINDEは、この誤差をなくすことを目的としている。dFINDEは、離散時間で定義される離散勾配を使用する。本論文は、｢微分系の離散アナログを、離散時間で積分することで、第一積分を、丸め誤差を含めて、正確に”保存”する｣と主張する。ただし、｢dFINDE は陰解法であるため、比較的計算コストが高くなる｣。
離散勾配法を使用して、第一積分を保存するモデルとして、(本論文の研究者によって開発された)DGNetも存在する。ただしDGNetが保存する第一積分は、ハミルトン系でのハミルトニアン(つまりエネルギー)に限定される。FINDEは、非ハミルトン系の幅広い第一積分を保存する。なお、離散勾配は、自動離散微分アルゴリズム(ADDA)によって得られる。ADDAは、PINNsで重要な役割を果たす自動微分の離散バージョンと考えられる。
【4】セットアップ
本論文では、FINDEと基本モデルの比較評価を行っている。
(1) 評価対象となった物理系、及び基本モデル
評価対象の次の4つである：❶重力2体問題、❷KdV方程式、❸二重振り子、❹フィッツヒュー･南雲方程式。
❶は、典型的な(正準形式の)ハミルトン系である(ハミルトンベクトル場で定まる力学系をハミルトン系と呼ぶ)。ハミルトニアンに加えて、(ハミルトン系であるためネーターの定理から)空間の並進対称性と回転対称性に対応して、運動量と角運動量という第一積分が存在する。❷は、代表的なソリトン方程式である。非正準形式のハミルトン系で、多くの第一積分を持つ。❸は、ポアソン系(縮退したハミルトン系)である。ハミルトニアン以外に、4つの第一積分を持つ。❹のフィッツヒュー･南雲方程式は、生体ニューロン(神経細胞)の定性的な振る舞いをモデル化した方程式である。2つの第一積分を持つ。
比較対象となった基本モデルは、HNNと(やや曖昧であるが)NODEである。
(2) 環境
1⃣ パッケージライブラリ
Scipy(Pythonの科学技術計算用ライブラリ)、Pytorch(Pythonの機械学習用ライブラリ)、torchdiffeq(Pytorch用のNODE実装ライブラリ)、functorch(自動微分を含む関数変換をPytorchに追加するライブラリ)、および gplearn(シンボリック回帰を実行する遺伝的アルゴリズムを実装した機械学習ライブラリ)を使用した。
2⃣ その他
㊀ハードウェア・・・ NVIDIA A100
㊁数値積分法・・・Dormand-Prince法
㊂活性化関数・・・双曲線正接関数
㊃バッチ サイズ・・・ 200
㊄オプティマイザー・・・Adam
㊅学習率とスケジューリング・・・10−3から、コサイン･アニーリングでゼロに減衰。
(3)評価指標
1ステップ誤差と有効な予測(VPT)を採用している。ただし、1ステップ誤差はかませ犬であり、誤解を招く悪い指標として採用されている。1ステップ誤差は、予測値と実測値の平均自乗誤差(MSE)として定義されている。 この指標は、値が低いほど良い。
VPTは、(初期値問題において)予測状態のMSE が、最初に所定のしきい値を超える時系列の長さ で除算された時点を示す。 この指標は、高いほど良い。
【5】本論文の成果
(1) 基本モデルとの比較
❶重力2体問題の基本モデルは、HNN。それ以外(❷～❹)の基本モデルはNODE。FINDEは、cFINDEとdFINDEの両方が比較対象。VPT(試行5回の中央値)を評価指標として、❶～❹の全てで、dFINDEが最も性能が高かった。基本モデルとは、大差を付けている。❸の二重振り子では5倍以上の差を付けた。
(2) 第一積分の発見
第一積分の(候補)数Kに対応するVPTを並べていき、最も高いVPTを発生させるKが、第一積分の(真の)数であるという解釈をしている。
❶重力2体問題の場合、Ｋ＝２のVPTが最大であり、運動量と角運動量の２つの第一積分を発見できた、と結論している。❷KdV方程式の場合は、第一積分は多数あるが、cFINDEで３つ、dFINDEでは４つ見つかったことになる(本論文では、１つ目若しくは２つ目の第一積分が全質量で、３つ目の第一積分がエネルギーであろう、と結論している)。
❸２重振り子ではＫ＝5のVPTが最大かつK=6でVPTが急落しているので、5つの第一積分を、全て発見できたとしている(急落の度合はcFINDEでおよそ1/117、dFINDEでおよそ1/15)。❹フィッツヒュー南雲方程式は、Ｋ＝2のVPTが最大かつK=3のVPTが急落しているので、2つの第一積分を、全て発見できた、としている(急落の度合はcFINDEでおよそ1/62、dFINDEでおよそ1/9)。
(3) ダイナミクスの再現
❶、❸、❹で予測状態が実際のダイナミクスを再現することを確認している。❶では、重心位置と質量位置が、HNNに比べれば、実測値とFINDEで同じような挙動を示すことを確認(HNNはchaoticな動きになる)。❸では、2つの振り子の位置が、実測値とFINDEで同じような軌道を辿ることを確認(NODEはややchaoticな動きになる)。❹では、周期軌道に収束する実測値とFINDEで同じ挙動を示すことを確認している(NODEは収束しない)。
【6】感想
(1) 可積分系
FINDEは、関数的に独立な第一積分を考えている。｢関数的に独立｣とは、関数の勾配が、稠密な開部分集合上で一次独立であることを言う。実は、関数的に独立かつ、包合的な第一積分が存在すれば、可積分である。もう少し正確に言うと、ℝ2nの領域Dで定義されたハミルトニアンを持つハミルトンベクトル場が、n個の｢関数的に独立｣かつ、｢包合的｣な第一積分を持つとき、(リウヴィル)可積分である(もっと正確に言うと、ユークリッド空間ℝ2nに限定されず、シンプレクティック多様体に対して成立する)。｢包合的｣とは、任意のポアソン括弧式が恒等的にゼロになること(ポアソン可換)を言う。
ハミルトン系であれば、シンプルな条件で可積分であることが言える。ハミルトン系を含む幅広い物理系を対象とするFINDEが、もし可積分系の研究を押し広げることができれば、興味深い。
(2) その他分野への展開
FINDEは元々、未知システムへの応用をモチベーションとしているようであるから、金融分野への展開なども興味深い(PINNsは、既に、金融において使われている)。
(3) 分かりにくい?
❶重力2体問題は、K=2でVPT最大となっているが、K=3でVPT最大でなければおかしいのでは?・・・と一瞬思うが、❶の基本モデルはHNNで、既にエネルギー(ハミルトニアン)を保存するように設計されているので、K=2で良い。しかし、分かりにくい?
【0】はじめに
科学計算総合研究所(RICOS)及び筑波大学の研究者は、物理埋め込み型ニューラルネットワーク(PENN)という計算手法を提案･発表した(23年3月23日、arXivにて論文(以下、本論文)[*A-3]公開)[*A-7]。
PENN誕生までの過程を記すと、(1)グラフベースの畳み込みニューラルネットワーク(GCN)→(2)同変性を満たすGCNであるIsoGCN→(3)境界条件を厳密に満たす＋大域的相互作用を捉えたIsoGCNであるPENN、となる。
【1】本論文の訴求ポイント
(1) PENNは、ディリクレ並びにノイマン境界条件を厳密に満たすので、信頼性が高い。
(2) PENNは、長時間経過後の状態を正確に予測することができる。
【2】事前整理
(1) 先行研究の整理
1⃣ Pythics-Informedニューラルネットワーク(PIN)との比較
PINNsの欠点を、次のように記述している。
❶ PDE(偏微分方程式)の解は、領域の形状や境界条件が変わると、全く異なるものになる可能性がある。これは、モデルの汎用性を大きく制限する。
❷ PINNsモデルによる予測は、物理(則)から外れることがある。
一方PENNは、同変性(後述)を満たすので、物理量は座標変換に対して不変(座標変換と同じ変化をするので不変＝同変)という｢物理の前提｣を満たすため、高い汎用性を実現している、と主張している。
2⃣ GNNを用いたPDEソルバーとの比較
GNN(グラフ･ニューラルネットワーク)は、境界条件を入力としてモデルに与えることで、境界条件を考慮することができる。しかし、本論文は、｢ディリクレ条件のような難しい制約を満たす保証はない｣と主張する。対照的に、PENNはディリクレ条件を(含む境界条件を)満たすことを保証する。※ディリクレ条件は、境界上の値そのものを与える、という境界条件である。
また、｢多くのGNNは、固定数のメッセージ･パッシングによる局所的な接続を利用しており、大域的な相互作用の考慮が欠けている｣が、PENNは、GNNに大域的な接続を組み込むことができる、と主張している。
(2) PENNに至るまで
1⃣ GCN以前
PENNの大元を成す計算アルゴリズムは、有限要素法(FEM)である。工学的には、FEMは構造解析(より正確には航空機の設計分野)で、最初に提案された。最初の論文が発表されたのは、1956年と言われている。その当時、FEMは、業界の長老たちからcompletely wrongと指弾された、と噂されている。流体解析には、60年代末～70年代始にかけて、最初の研究発表が行われている。
FEMは解析領域を、大きさが様々な単純図形(≃三角形)を数多く使って、分割近似する。この単純図形の集合体は、網目に見えることから(通常)、メッシュと呼ばれる。昔は、メッシュを切る作業が、大変だった。
メッシュを(ジュリウス・ピーターセンが命名したところの)グラフと考えると、メッシュを数学的な対象にも、機械学習の対象にもできる(ちなみに、無向グラフ)。機械学習の文脈で言うと、グラフニューラルネットワーク(GNN)あるいは、グラフ畳み込みニューラルネットワーク(GCN)の枠組みで、メッシュを扱うことができる。
2⃣ GCN：グラフ理論的整理
グラフは、頂点(英語ではvertexヴァーテックスあるいは、nodeノード)の集合と、(頂点のペアを繋ぐ)辺(枝とも言う。英語ではedgeエッジあるいは、arcアーク)の集合から構成される。グラフを表現するデータ構造は、接続行列と隣接行列である。
無向グラフの接続行列は、行が頂点集合で、列が辺集合に対応している。頂点が辺の端点であれば(通常)、対応する行列要素(あるいは成分)は1。端点でなければ0とする。例えば、頂点1が、頂点2と辺1を、頂点3と辺2のみを作っている場合、1行1列の行列要素は1。1行2列の行列要素も1。それ以外の1行の列要素は0となる(言わずもがな、文字で書くと分かり辛いが、図で書くと分かり易い)。尚、有向グラフの場合、向きがあるので、端点でも始点と終点を区別する。そのため、行列要素には、1,0,-1の３つが割り振られる。
無向グラフの隣接行列は、行と列の、それぞれに、頂点集合が対応している。行列要素は(通常)、頂点間に辺が存在していれば1。存在していなければ0とする(有向グラフの場合、辺が存在するのみならず、｢向き｣まで考慮する)。例えば、1行1列の行列要素に該当する頂点は、(番号1を付けた)頂点となる。(同じ頂点1の間に辺は、ないから)1行1列の行列要素は当然0になる(つまり対角要素は常に0)。1行2列の行列要素に該当する頂点は、頂点1と頂点2である。頂点1と頂点2が繋がっていれば、行列要素の値は1となり、繋がっていなければ0となる。"頂点1と頂点2が繋がっていれば、頂点2と頂点1はつながっている"のだから、無向グラフの隣接行列は、対称行列となる。
一方、有向グラフは、向きまで考えるから、行列要素(i,j)の値は(j,i)の値と同じにはならない。というより、スタビライザー符号のように、(i,j)の値が1なら、(j,i)の値は0にmなる。(i,j)の値が0なら、(j,i)の値は1になる。0の場合は、もちろん、変わらず0である。言わずもがなであるが、行列要素(i,j)はi行j列の行列要素という意味である。
接続行列を使うにせよ、隣接行列を使うにせよ、行列形式でグラフ(今の場合、正確にはメッシュ)を扱えるようにすれば、機械学習(今の場合、GCN)の土俵に乗せること自体は、難しくはない。通常は、隣接行列を使う。
余談ながら、接続行列の領域計算量はO(mn)、隣接行列はO(n2)である。ここで、mは辺の数、nは頂点の数である。
3⃣ GCN以後：同変性
GCNからIcoGCNのステップアップは、同変性の考慮である。これは、｢物理現象は、座標変換に伴って変化しない｣から、それをシミュレーションにも反映させる、ことを意味している。手続き論的に素朴に表現すると、シミュレーションを実施した後に、座標変換(例：並進変換、回転変換)を行った結果と、座標変換を行った後にシミュレーションして得られる結果は、同じであることを意味している。
本論文によると、GCNは、同変性を満たさない。一方、畳み込みニューラルネットワーク(CNN)は並進変換については、同変性を満たす[*A-4](本論文にも、同変性を満たすようにCNNを改良して、流体解析に適用した例が示されている)。CNNは畳み込み操作を工夫することで、同変性を導入する。GCNではメッセージ･パッシングを使って、畳み込み操作を行う。このため、IsoGCNではメッセージパッシングにおいて工夫を施すことで、同変性を導入していると考えられる。IcoGCNでは、同変性を満たすために、頂点の絶対位置で隣接行列を構成するのではなく、頂点の相対位置で隣接行列を構成している[*A-5]。
また、この隣接行列(重み付き隣接行列)は、最小二乗移動粒子半陰解法(LSMPS法)の微分作用素(勾配∇)に対応している。このため、様々な空間微分作用素に対応するIsoGCNを構成することができる。結論として、IsoGCNはPDEの学習に適している、と主張する[*A-6]。
なお、PENNの名称に現れる｢埋め込み｣とは、同変性が考慮されていることを意味している(と思われる)。
(3) PENNの特徴－IcoGCN以後
本論文によると、IcoGCNの課題は以下の通り：
㊀境界条件の扱いが不十分。ディリクレ境界条件の厳格な履行は不可欠である。
㊁局所的に接続されていることは、汎化性能を高めるためにはプラスに作用する。その反面、局所的であるが故に、大域的な相互作用を見逃す可能性がある。長時間後の状態を予測するために、大域的相互作用を捉える必要がある。
1⃣ 境界条件の厳格な履行
本論文によれば、｢ディリクレ境界条件を満たすソルバーを得るためには、変数と境界条件の間で、同じ符号化空間を確保する必要がある。同じ符号化空間を確保するために、変数とそれに対応するディリクレ境界条件に対して、同じエンコーダ(本論文では、境界エンコーダと呼ばれている)を使用する｣。また、｢符号化された境界条件が、元の物理空間における境界条件と一致することを保証するデコーダー(本論文では、疑似逆デコーダーと呼ばれている)は、広い範囲のニューラルネットワークで構成可能｣と主張している。
ノイマン境界条件に関しては、LPMPS法によるノイマン境界条件モデル(壁面境界モデル)を用いて、IcoGCNを修正している。IcoGCNが、LSMPS法の微分作用素に対応して(構成されて)いるから、LPMPS法のモデルを使用している。※ノイマン条件は、境界上の値の微分値(≃勾配)を与えるという境界条件である。
2⃣ 大域的相互作用を捉える
今更のように、PDEをニューラルネットワークで解くとは、PDEの解法を最適化問題の解法に置き換えることである、というお題目を思い出そう。PENNもニューラルネットワークである。最適化問題は、お決まりの勾配降下法を使って解く(正確には、少ない反復回数で収束することが知られている｢Barzilai–Borwein法｣と呼ばれる勾配降下法を使う)。
この勾配降下法(Barzilai–Borwein法)のステップサイズは、メッシュ上の内積を使って、計算することができる。内積はメッシュ全体で計算されるため、本論文では｢勾配降下のステップサイズの計算は、大域的プーリングに相当する」と主張する。結果として、大域的な相互作用を捉えることを可能にした、と主張する。
【3】実験及び成果
実験は2種類行われている。PENNが、㊀与えられたスカラー場から勾配場を予測できるか、という実験及び、㊁複雑な形状を流れる非圧縮性流体のダイナミクスを学習できるか、という実験である。
以下では、㊀は割愛して、㊁のみ整理した。
(1) 評価対象モデル、テストデータ、並びに評価指標
比較対象モデルは、メッセージ･パッシング･ニューラルPDEソルバー(MP-PDE)である。また、テストデータは、プレーンなソルバーを使って、以下のようにして作成している。
❶ソルバー＝OpenFOAM
❷タイムステップ＝10-3
❸初期条件＝ポテンシャル流の解
❹レイノルズ数＝1000
❺アウトプット＝4秒後の流速、圧力
❻テストデータ＝アウトプットを、ランダムに並進変換あるいは、回転変換させて作成。
❼評価指標：流速と圧力の平均自乗誤差(MSE)
(2) 結果および、アブレーション分析
PENNは、並進変換あるいは回転変換前後のテストデータで、流速、圧力共に、評価指標(MSE)が変化していない。MP-PDEでは、変換前後で、MSEが大きく変動している。
アブレーション分析の結果、大域的相互作用の取り込み、が一番効いていた。
【4】考察
(1) 使い回し性
PENNはサロゲートモデル(代理モデル)を前提として(かつ、流体解析にフォーカスして)いるから、汎用性、俗な言い方をすると｢使い回し性｣に拘っている。使い回し性を高めるためには、確かに、同変性が必要と考えられる。メッシュ作成は、時間コストが高い。解析対象を回転させたくらいで、改めてメッシュを切り直していたら、サロゲートモデルの意味がない。ただ、時間がかかっても(初期投資が重くても)、使い回しさえできれば、回収はできる。
一方PINNsは、領域の形状や境界条件が変わるたびに、学習仕直す必要がある。使い回し性が劣ることは否めない。その反面、(例えば、金融分野が該当すると思われる)使い回し性が要求されないケースでは、PENNの優位性は減じられるであろう。また、本論文でも指摘されているように、PINNsが得意とされる逆問題は、PENNでは対応できない。その意味でも、金融においては、PENNよりPINNsが有用性が高いだろう[*A-8]。
(2) PINNsとの融合
またPINNsでは、学習データを与える時空間座標{ti,xi}におけるu(ti,xi)と、学習データの値uiが異なる影響を、損失関数を最小化することで、低減させている。それが不十分で、PINNsが物理則から外れるというケースは、あるだろう。なお、ここでuは、PDEの近似解である。もちろん、同変性の考慮だけで、PENNが物理則から外れることを回避(して、シミュレーションで高い性能を発揮)できるわけではない。
(非圧縮性流体を解析する場合に必要と言及されている)大域的相互作用の取り込みに関しては、PINNsと同じ取り扱いをしている。つまり、(右辺＝0の形式での)支配方程式を損失関数とする最適化問題に置き換えて、ニューラルネットワークに取り込んでいる。その意味では、PINNsとPENNの精神的融合が、行われていると言えるかもしれない。
いずれにしても、数値流体解析が必要な現場では、PENNは有力なツールと考えられるだろう。
【尾注】
*1 Salvatore Cuomo et al.、Scientific Machine Learning through Physics-Informed Neural Networks: Where we are and What’s next、https://arxiv.org/pdf/2201.05624.pdf
*2 George Em Karniadakis et al.、Physics- informed machine learning、https://www.brown.edu/research/projects/crunch/sites/brown.edu.research.projects.crunch/files/uploads/Nature-REviews_GK.pdf
*3 Raissi M, Perdikaris P, Karniadakis GE、Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations、https://www.sciencedirect.com/science/article/pii/S0021999118307125
この論文は3442回引用されており、PINNs論文の中で、ダントツで引用数が最も多い。
*4 Tim De Ryck, Siddhartha Mishra、Error analysis for physics informed neural networks (PINNss) approximating Kolmogorov PDEs、https://arxiv.org/abs/2106.14473
*5 人工知能技術は一般に、内挿は得意だが、外挿は不得意であることが多い。例えば、物性値予測での外挿とは、学習用のデータと分子構造が大きく異なる化合物の物性予測などである。外挿予測は新規の材料や薬剤の開発に極めて重要である。
出典：https://www.aist.go.jp/aist_j/press_release/pr2020/pr20201111/pr20201111.html
*6 Lu Lu et al.、DeepXDE: A Deep Learning Library for Solving Differential Equations、https://epubs.siam.org/doi/epdf/10.1137/19M1274067
*7 https://www.riken.jp/press/2022/20221201_2/index.html
*8 同じ地球物理分野で、PINNsに対して｢物理情報ニューラルネットワーク｣という訳語をあてている記事もある：磯真一郎、地球科学分野における機械学習の近年の研究動向について、深田地質研究所年報、No.22、p.155-164(2021)
*9 https://www.jsap.or.jp/docs/pressrelease/JSAP-2022autumn-chumoku-09.pdf
*10 欠番
*11 Andreas Louskos、Physics-Informed Neural Networks and Option Pricing、https://math.dartmouth.edu/theses/undergrad/2021/Louskos-thesis.pdf
*12 岡野豊明、アメリカン・オプションの理論価格の計算法、UNISYS TECHNOLOGY REVIEW第61号, MAY1999、https://pr.biprogy.com/tec_info/tr61/6106.pdf
*13 Johannes Ruf、Weiguan Wang、Neural networks for option pricing and hedging: a literature review、https://arxiv.org/pdf/1911.05620.pdf
*14 Soohan Kim et al.、PHYSICS-INFORMED CONVOLUTIONAL TRANSFORMER FOR PREDICTING VOLATILITY SURFACE、https://arxiv.org/pdf/2209.10771.pdf
*15 Ramy Tanios、Physics Informed Neural Networksin Computational Finance: High Dimensional Forward & Inverse Option Pricing、https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/491555/1/Thesis.pdf
Webで閲覧できるバージョンは未完成版のようである(ETHZ公式であろうhttps://www.research-collection.ethz.ch/aboutから辿っても、同じ結果)。
*16 CEVモデルとヘストンモデルで、インプライド･ボラティリティを算出している。
*17 この修論のポイントは、PINNsが、高次元(100次元)ブラックショールズモデル及びヘストン･モデルのヨーロピアン･コールオプションを評価できたことである。つまり次元の呪いを避けられた、との主張である。
*18 Zaharaddeen Karami Lawal et al.、Physics-Informed Neural Network (PINNs) Evolution and Beyond: A Systematic Literature Review and Bibliometric Analysis、https://www.mdpi.com/2504-2289/6/4/140
*19 Alexandr Sedykh et al.、Quantum physics-informed neural networks for simulating computational fluid dynamics in complex shapes、https://arxiv.org/pdf/2304.11247.pdf
*20 幾何形状のバリエーションを、シミュレーションの再実行ではなく、転移学習によってカバーしている。具体的には、Y字型流路のY字部分の角度を変更する場合、1000エポック・Adamで事前学習したPINNsを、100エポック・L-BFGSで再学習する。
*21 Lucas Friedrich and Jonas Maziero、Quantum neural network cost function concentration dependency on the parametrization expressivity、https://www.nature.com/articles/s41598-023-37003-5
*22 Stefano Markidis、On physics-informed neural networks for quantum computers、https://www.frontiersin.org/articles/10.3389/fams.2022.1036711/full
*23 https://sumitomoelectric.com/jp/press/2023/08/prs103
*24 Ryu Shirakami et al.、QTNet: Theory-based Queue Length Prediction for Urban Traffic、https://dl.acm.org/doi/pdf/10.1145/3580305.3599890
*25 主に、以下を参考にした：挺屋 友幹他、時間・空間的相関を考慮した人流の時系列予測手法、https://ipsj.ixsq.nii.ac.jp/ej/?action=pages_view_main&active_action=repository_action_common_download&item_id=209591&item_no=1&attribute_id=1&file_no=1&block_id=8&page_id=13
*26 Guangyin Jin et al.、Spatio-Temporal Graph Neural Networks for Predictive Learning in Urban Computing: A Survey、https://arxiv.org/pdf/2303.14483.pdf
*27 Haakon Robinson et al.、Physics guided neural networks for modelling of non-linear dynamics、https://www.sciencedirect.com/science/article/pii/S0893608022002854
*28 元データにおける情報に対する、特徴マップ上の対応範囲が広い、ことを受容野が広いという。元々は、脳の視覚野で、視覚刺激に対して応答する範囲を受容野と呼ぶところから来ている。CNNは画像認識から始まり、その文脈で使われたことから、時系列データをCNNで処理する場合は、全て受容野という用語が使われているようである。
*29 岡田 治他、AIによる洪水予測と水位データの異常検知、http://www.river.or.jp/report_sympo_h2905.pdf
*30 緒方陸他、交通速度の予測におけるGraph Neural Network学習の効率化に関する研究、https://www.jstage.jst.go.jp/article/pjsai/JSAI2023/0/JSAI2023_1M3GS1002/_pdf
*31 DCRNN(Diffusion Convolutional Recurrenct Neural Network)は、交通量予測のために考案されたGNN。交通量を拡散過程と関連付けることで、空間依存性をモデル化する。参照先は、https://db-event.jpn.org/deim2019/post/papers/315.pdf
*32 AGCRN(Adaptive Graph Convolutional Recurrent Network)は、｢ノード固有のパラメータ空間を生成し、データからノード埋め込みを推論することで、ノードの隣接行列を学習によって得るモデル｣である。参照先は、[*25]。
*33 MegaCRNは、メタ･グラフ学習とエンコーダー-デコーダー型グラフ畳み込み再帰型ネットワークとを組み合わせた、ニューラルネットワーク。トヨタと東京大学が共同で開発した。参照先は、https://www.toyota-tokyo.tech/news/pdf/20230201pressrelease_jp.pdf及びRenhe Jiang et al.、Spatio-Temporal Meta-Graph Learning for Traffic Forecasting、https://arxiv.org/pdf/2211.14701.pdf
*34 Salah A. Faroughi et al.、Physics-Guided, Physics-Informed, and Physics-Encoded Neural Networks in Scientific Computing、https://browse.arxiv.org/pdf/2211.07377.pdf
*35 https://www.jamstec.go.jp/j/about/press_release/20231011_2/
*36 以下を参考にした：松田･宮武、数値解析と確率･統計による不確実性定量化、https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2167-07.pdf
*37 以下を参考にした：Liu Yang et al.、B-PINNss: Bayesian Physics-Informed Neural Networks for Forward and Inverse PDE Problems with Noisy Data、https://arxiv.org/pdf/2003.06097.pdf
*38 例えば、Qiang Liu et al.、Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm、https://arxiv.org/pdf/1608.04471.pdf
*39 例えば、宮本、変分推論の理論、https://www.ccn.yamanashi.ac.jp/~tmiyamoto/img/variational_bayes1.pdf
*40 Sifan Wang et al.、AN EXPERT’S GUIDE TO TRAINING PHYSICS-INFORMED NEURAL NETWORKS、https://arxiv.org/pdf/2308.08468.pdf
*41 少なくとも日本では、同大学は、ウォートン･ビジネス･スクールの存在によって有名である。なお、2023年のノーベル生理学・医学賞を受賞したカタリン・カリコ氏とドリュー・ワイスマン氏は、同大学の非常勤教授と教授である。
*42 https://developer.nvidia.com/ja-jp/blog/enhancing-digital-twin-models-and-simulations-with-nvidia-modulus-v22-09/
*43 Stanのサポートに加えて、ニューラル演算子の実装が拡張されている。ニューラル演算子は、DeepONet、フーリエ･ニューラル演算子、Physics-Informedニューラル演算子、である(グラフ･ニューラル演算子はカバーされていないようである)。英Graphcoreや米SambaNovaが、NVIDIAに対抗してGPUを開発しているが、ソフトウェア･ライブラリの充実度やカバー範囲を見る限り、勝負にならないだろう。
*44 Raghav Gnanasambandam et al.、Self-scalable Tanh (Stan): Faster Convergence and Better Generalization in Physics-informed Neural Networks、https://arxiv.org/pdf/2204.12589.pdf
*45 Zongyi Li et al.、FOURIER NEURAL OPERATOR FOR PARAMETRIC PARTIAL DIFFERENTIAL EQUATIONS、https://arxiv.org/pdf/2010.08895.pdf
*46 Somdatta Goswami et al.、Physics-Informed Deep Neural Operator Networks、https://browse.arxiv.org/pdf/2207.05748.pdf
*47 Lu Lu et al.、DeepONet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators、https://arxiv.org/pdf/1910.03193.pdf
*48 Kamyar Azzizadenesheli et al.、Neural Operators for Accelerating Scientific Simulations and Design、https://arxiv.org/pdf/2309.15325v1.pdf
*49 Rajat Arora、PhySRNet: Physics informed super-resolution network for application in computational solid mechanics、https://arxiv.org/pdf/2206.15457.pdf
*50 Gege Wen et al.、U-FNO - an enhanced Fourier neural operator-based deep-learning model for multiphase flow、https://arxiv.org/pdf/2109.03697.pdf
*51 Samuel Lanthaleret al.、NONLINEAR RECONSTRUCTION FOR OPERATOR LEARNING OF PDES WITH DISCONTINUITIES、https://openreview.net/pdf?id=CrfhZAsJDsZ
*52 Faisal Amlani & Oscar P. Bruno、An FC-based spectral solver for elastodynamic problems in general three-dimensional domains、https://ui.adsabs.harvard.edu/abs/2016JCoPh.307..333A/abstract
*53 小山田耕二、ニューラルネットワークを用いた偏微分方程式導出・求解を支援する視覚的分析、https://axies.jp/_files/conf/conf2021/paper/WD5-2.pdf
*54 https://www.hpcwire.com/off-the-wire/multiverse-computing-wins-uk-funding-to-improve-flood-risk-assessment-with-quantum-algorithms/及びhttps://thenextweb.com/news/quantum-startup-multiverse-computing-predict-floods
*55 ヨシュア･ベンジオが設立した研究所。ヨシュア･ベンジオは、｢深層学習のゴッドファーザー｣の一人。ゴッドファーザーは3人いて、他の二人は、ジェフリー･ヒントンとヤン･ルカン。3人は、2018年に(計算機科学のノーベル賞とも言われる)チューリング賞を受賞。アンドリュー･エンをいれて4大グルとも呼ばれる。アンドリュー･エン(Ng)は、アンドリュー･ングとか、アンドリュー･ンとも表記される。
ジェフリー･ヒントンは、自身が設立したスタートアップがグーグルに買収(2013年)された経緯から、2013～2023年までグーグルに席を置いた。ヒントンの弟子筋にあたるヤン･ルカンは、NECの北米研究所に所属していた時期もあるが、2013年以降、メタに席を置いている。
ヨシュア・ベンジオは、英国で11月1～2日にかけて開催されたAI安全サミットが、作成する｢AIの能力と将来のリスクについて共通の理解構築に向けた報告書｣の作成に協力する、と報道されている。
*56 Qidong Yang et al.、Fourier Neural Operators for Arbitrary Resolution Climate Data Downscaling、https://arxiv.org/pdf/2305.14452.pdf
*57 Paula Harder et al.、Physics-Constrained Deep Learning for Climate Downscaling、https://arxiv.org/pdf/2208.05424.pdf
*58 正確には、2次元非圧縮性(粘性)NS方程式で、かつ渦度形式のNS方程式を用いている。渦度形式のNS方程式とは、速度ではなく渦度(勾配∇と速度の外積)の時間発展を記述するNS方程式である。粘度は、10-4で設定されている。
NS方程式で非圧縮性とするのは、特殊なことではない(本論文のセットアップが、リアリティを欠いた設定ではない、という意味)。またNS方程式は粘性流体を前提としているので、通常、わざわざ粘性は付けない。
*59 Nikola Kovachki et al.、Neural Operator: Learning Maps Between Function Spaces With Applications to PDEs、https://www.jmlr.org/papers/volume24/21-1524/21-1524.pdf
*60 Zongyi Li et al.、Fourier Neural Operator with Learned Deformations for PDEs on General Geometries、https://arxiv.org/pdf/2207.05209.pdf
*61 Ethan Pickering et al.、Discovering and forecasting extreme events via active learning in neural operators、https://arxiv.org/pdf/2204.02488.pdf
*62 日野英逸、特集｜能動学習：問題設定と最近の話題、日本統計学会誌 第50巻, 第2号, 2021 年3 月,317頁 ∼ 342頁、https://www.jstage.jst.go.jp/article/jjssj/50/2/50_317/_pdf
*63 鈴木絢子･西浦博、感染症の数理モデルと対策、https://www.naika.or.jp/wp-content/uploads/2020/11/nichinaishi-109-11-article_4.pdf
*64 https://www.nikkei.com/article/DGXZQOUD122RY0S2A610C2000000/
葛飾北斎の富嶽三十六景･神奈川沖浪裏に描かれている大きな波濤は、巨大波ではないか？とも言われているらしい。
*65 他は、中国の曁南大学と米パデュー大学。
*66 Ehsan Kharazmi et al.、Identifiability and predictability of integer- and fractional-order epidemiological models using physics-informed neural networks、https://www.nature.com/articles/s43588-021-00158-0.pdf
*67 鈴木清樹、疫学モデルを用いた病害の動態解析と防除対策、植物防疫 第63巻 第10号(2009年)、https://jppa.or.jp/archive/pdf/63_10_10.pdf
*68 杉本信正、整数階でない微積分法について 、https://www.mathsoc.jp/publication/tushin/2104/2016sugimoto.pdf
*69 岩山隆寛、非整数階の微積分、https://www.se.fukuoka-u.ac.jp/iwayama/teach/fractional_calculus/tsuuron_fractional_deriv.pdf
*70 Supplementary Information、https://static-content.springer.com/esm/art%3A10.1038%2Fs43588-021-00158-0/MediaObjects/43588_2021_158_MOESM1_ESM.pdf
*71 Ben Moseley et al.、Finite basis physics-informed neural networks(FBPINNs): a scalable domain decomposition approach for solving differential equations、https://link.springer.com/article/10.1007/s10444-023-10065-9
FB-PINNsのカジュアルな資料として、以下がある：https://maths4dl.ac.uk/wp-content/uploads/2023/07/Moseley.pdf
*72 Zhongkai Hao et al.、PINNACLE: A COMPREHENSIVE BENCHMARK OF PHYSICS-INFORMED NEURAL NETWORKS FOR SOLVING PDES、https://arxiv.org/pdf/2306.08827.pdf
*73 Zongren Zou et al.、Correcting model misspecification in physics-informed neural networks(PINNs)、https://arxiv.org/pdf/2310.10776.pdf
*74 Zhongkai Hao et al.、PINNACLE: A COMPREHENSIVE BENCHMARK OF PHYSICS-INFORMED NEURAL NETWORKS FOR SOLVING PDES、https://arxiv.org/pdf/2306.08827.pdf
*75 Jiachen Yao et al.、MultiAdam: Parameter-wise Scale-invariant Optimizer for Multiscale Training of Physics-informed Neural Networks、https://arxiv.org/pdf/2306.02816.pdf
*76 Zhiwei Fang et al.、ENSEMBLE LEARNING FOR PHYSICS INFORMED NEURAL NETWORKS: A GRADIENT BOOSTING APPROACH、https://arxiv.org/pdf/2302.13143.pdf
*77 Johannes Muller ＆ Marius Zeinhofer、Achieving High Accuracy with PINNs via Energy Natural Gradient Descent、https://proceedings.mlr.press/v202/muller23b/muller23b.pdf
*78 福永竜世他、Physics-Informed Neural Networks による不飽和浸透特性パラメータの逆解析、AI･データサイエンス論文集4巻3号,2023、pp.100-108、https://www.jstage.jst.go.jp/article/jsceiii/4/3/4_100/_pdf/-char/ja
*79 甘利俊一、解説:特集 適応･学習制御システムの新展開| 自然勾配学習法-学習空間の幾何学、計測と制御、第40巻第10号 2001年10月号、https://www.jstage.jst.go.jp/article/sicejl1962/40/10/40_10_735/_pdf
*80 https://seiya-kumada.blogspot.com/2018/03/fisher.html
*81 Alex Kaltenbach & Marius Zeinhofer、The Deep Ritz Method for Parametric p-Dirichlet Problems、https://arxiv.org/pdf/2207.01894v1.pdf
*82 長瀬准平･長沼大樹、自然勾配法を用いたニューラルネットワークの学習のスキップ接続による影響、https://www.ieice.org/publications/conference-FIT-DVDs/FIT2021/data/pdf/F-018.pdf
*83 横田理央、近似行列分解と分散深層学習、https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2167-03.pdf
*84 Shengze Cai ea al.、Physics-informed neural networks (PINNs) for fluid mechanics: A review、https://arxiv.org/pdf/2105.09506.pdf
・・・[*85]以降は、PINN2.htmlへ・・・
(*87 Katsiaryna Haitsiukevich & Alexander Ilin、IMPROVED TRAINING OF PHYSICS-INFORMED NEURAL NETWORKS WITH MODEL ENSEMBLES、https://arxiv.org/pdf/2204.05108.pdf)
(*88 Jose Florido et al.、Investigating Guiding Information for Adaptive Collocation Point Sampling in PINNs、https://arxiv.org/pdf/2404.12282)
*A-1 Takashi Matsubara and Takaharu Yaguchi、FINDE: NEURAL DIFFERENTIAL EQUATIONS FOR FINDING AND PRESERVING INVARIANT QUANTITIES、https://arxiv.org/pdf/2210.00272.pdf
*A-2 流体力学においては(幅広く)、保存則と基礎方程式(支配方程式)に対応関係が成立している。最もプリミティブには、質量保存：連続の式、運動量保存：運動方程式(ナヴィエ･ストークス方程式)である。
*A-3 Masanobu Horie and Naoto Mitsume、Physics-Embedded Neural Networks: Graph Neural PDE Solvers with Mixed Boundary Conditions、https://arxiv.org/pdf/2205.11912.pdf
*A-4 GCNも並進変換に限れば、同変ではないかと思うが、そう単純ではないのだろうか?
*A-5 絶対位置とは異なり、相対位置なら、並進変換や回転変換で変化しない、という素朴な理解で良いのだろうか。
*A-6 堀江正信･三目直登、物理現象の対称性をを持つグラフニューラルネットワークによる流動現象の学習、https://www.nagare.or.jp/download/noauth.html?d=41-6_403_tokushu7.pdf&dir=168
*A-7 PENNは、Physics-Embeded NNである。他方で、Physics-Encoding NN(PeNN)も存在する[前出*34]。PeNNは、物理法則をハード制約として課すことが特徴である。つまりPINNsのように、罰則項として物理学を課す(支配方程式をなるべく満たすようにソフト制約を課す)のではなく、物理法則を絶対に満たすようにガチガチに制約をかける(ハード制約)。このような性質から、PeNNのアーキテクチャ及び実装は複雑で、学習は困難となる。これが、PeNNのデメリットである。一方、PeNNのメリットとして｢データ不足に対する堅牢性、高い汎化性能｣が上げられる。｢収束率、安定性、スケーラビリティ、サンプル･サイズ｣に対する課題は、PINNsと同じである。
[*34]では、PeNNの代表例として、再帰型畳み込みニューラルネットワーク(RCNN)と、ニューラル常微分方程式(Neural ODE)が上げられている。なお、不変性(保存則)をハード制約するPENNと、物理法則をハード制約するPeNNとは、思想が似ている。
なお、学習のカテゴリーで言うと、PINNsはインスタンス学習(つまり、都度学習)に分類され、PeNNは継続学習に分類される(できる)。継続学習により、使い回し性を担保している、という見方ができるだろう。
*A-8 PeNNの範囲で言うと、時系列データの作成に、ニューラルODE(Ordinary Differential Equation)が使われているらしい(出所：https://www.jstage.jst.go.jp/article/jsaisigtwo/2022/FIN-028/2022_78/_pdf)
