内容 |
表示
不確実性評価が可能な新しい次元変分法New Fourdimensional Variational Method that Enables Uncertainty Quantiˆcation長尾大道(ながお東京大学地震研究所ひろみち)准教授伊藤伸一(いとう東京大学地震研究所しんいち)特任研究員. は じ め に理学分野及び工学分野において,物理や化学等の理論に基づいて構築されるモデルと,観測や実験によって得られるデータを比較検討することが研究推進の駆動力となることは,古今東西を問わず論をまたない。しかしな( a )従来の 4 次元変分法,( b )不確実性評価が可能な 4 次元変分法,( c )初期値依存性を排した 4がら,近年の数値シミュレーションモデル(以下,数値次元変分法のイメージ図―モデル)の大規模化,及び観測・実験データ(以下,観測データ)の大容量化に伴い,両者を単純に比較するこ系列データに対して事後分布を構成し,それを最大化すとすら困難となってしまっている。データ同化は,このる解を勾配法によって探索する手法である。事後分布のような大規模数値モデルと大容量観測データをベイズ統勾配を計算するためのアジョイントモデルを構成する必計学に基づいて融合する計算技術であり1),2),特に現代要があるため,数式展開やプログラム実装の手間が煩雑の数値天気予報はデータ同化抜きには成立しないと言っとなるが,計算量やメモリ量を数値モデルと同程度に抑て良い3)。データ同化は,系の状態やパラメータの推定えることが可能であるため,大自由度の数値モデルに対に加えて,将来予測のための定量的手段を与えるため,するデータ同化法として極めて有用であり,実際,数値地震学,生物学,あるいは材料科学等の様々な分野へと,天気予報における主力のデータ同化法として用いられてその応用範囲を広げている4)~6)。このように述べると,いる。しかしながら,従来の四次元変分法では,事後分データ同化は,従来の逆解析あるいは状態空間モデリン布を最大化する状態空間内の 1 点のみの情報しか得らグと何が違うのかという疑問が湧くかもしれない。デーれないため,推定値の不確実性を評価することは原理的タ同化の専門家の間でも見解が一致しているわけではなに不可能であった(図―(a))。例えば,台風の進路予いと思われるが,筆者は,データ同化の場合は大規模モ測の際に示される予報円は,台風の中心位置に関する不デルと大容量データの存在が前提となっていることが決確実性を表しているが,これは非逐次型データ同化と逐定的な違いだと考えている。有限の計算機資源によって,次型データ同化を組み合わせることによって算出してい現実的な計算時間内でデータ同化計算を実施するために,る。不確実性評価は,例えば興味あるパラメータを目的創意工夫を凝らした様々なデータ同化アルゴリズムが提の精度で推定するためには実験を何回繰り返せば良いか,案されており7),これは従来の逆解析や状態空間モデリあるいは場の量をより良く推定するためには,どこに計ングでは,十分には意識されていなかったと思われる。測センサーを設置すれば良いか等,観測・実験計画に対データ同化は,逐次型と非逐次型に大別される。アンしても重要なフィードバックを与える。大自由度の数値サンブルカルマンフィルタ8) や粒子フィルタ9) に代表さモデルに基づくデータ同化の需要が高まっていることかれる逐次ベイズフィルタを用いる逐次型データ同化は,ら,逐次型データ同化に頼ることなく,4 次元変分法の数値モデルの時間発展を計算しながら観測データを取り枠組みで不確実性評価が可能なアルゴリズム(図―込むことにより,状態変数を修正していく流儀であり,(b))が求められることは想像に難くない。多数の実現値(アンサンブルメンバー)の集合によって,本稿では,不確実性評価が可能な新しい 4 次元変分確率密度関数(分布)を近似表現する。そのため,アン法10) について紹介する。章で一般的な自励系モデルサンブルメンバー数を十分大きくすれば,分布の平均だに対する従来の 4 次元変分法について概説した後,けでなく,不確実性を表す標準偏差や,さらに高次の情章で 2nd order adjoint 法11) に基づく不確実性評価法を報を自然に得ることができる。しかしながら,全計算量定式化する。章では,構造材料内での組織形成を表現がアンサンブルメンバー数と数値モデルの計算量の積にする数値モデルとして用いられるフェーズフィールドモ比例するため,数値モデルの自由度が大きくなるにつれ,デル12),13) に基づく数値実験により,本手法の妥当性を現実的な計算時間内でデータ同化を実施するために多大示す。な工夫を要する。一方,非逐次型データ同化の代表的手法である 4 次元変分法2)は,あらかじめ得られている時2地盤工学会誌,―()論.説以上をまとめると,4 次元変分法による最適化は,次元変分法 式( 7 )~( 9 )により,x0 に適切な初期値を設定する,状態空間モデルが,次のように与えられたとする。 適当な勾配法を用評価関数の勾配 &J / &x0 を求める,&x(t)=f(x(t)) ………………………………………(1)&tり返す,という手順となる。ここで注意したいのは,y(ts)=h(x(ts))+w(ts),で得られる評価関数の勾配は「自動微分」とも呼ばれ14),w(ts)~r(・) ……………(2) 停止条件に至るまで と を繰いて x0 を更新する,式(1)はシステムモデルと呼ばれ,N 次元の状態 x(t)∈x0 の近傍における評価関数の差分で勾配を近似する数f R N → RN値微分に比して,計算精度及び計算量の両面において圧RNの時間発展を記述する。すなわち,関数は,数値モデルそのものである。また,式(2)は観測モデルと呼ばれ,状態 x(t) と時刻 t=t1, …, tn で得られたK 次元時系列データ y ( ts )∈ RK を比較することを意味倒的に優れていることである。.2nd order adjoint 法に基づく不確実性評価している。ただし,関数 h RN → RK は x ( ts ) から観測データと比較可能な量を抽出する観測演算子であり,また w(t) はある分布 r(・) に従う観測ノイズである。さて,4 次元変分法の目的は,事後分布p(x)r( y(ts)-h(x(ts)))………………(3)p( y) s_=1np(x|y)=を最大にする x(t) を探索することにある。まず,p( y)章で述べた 4 次元変分法により,事後分布を最大(評価関数を最小)にする初期状態 x0= âx0 を得たとする。この推定値の不確実性を評価するため,事後分布の x=âx0 の近傍を正規分布によってラプラス近似すると,p(x0|y)は所与の観測データに関する分布であり,単なる定数と1(2p)N/2|H-1|1/2×exp見なすことができる。また, p (x ) は事前分布であり,[-1(x0- âx0)TH(x0- âx0)2]……(10)系に関する事前知識を与える。ここでは簡単のため,一ここで,H∈RN×N は J に関するヘッセ行列 H=&2J/&x02様分布であると仮定するが,実問題においてこれは特別である。すなわち,事後分布を近似する正規分布は,に強い仮定というわけではないことが多い。したがって,N( âx0, H-1) となるので, âx0 の第 l 成分 âx0,l に関する不確事後分布最大化において本質的に意味を持つのは,式実性は, H-1 の第 l 対角成分の平方根 Hll-1 で与えら(3)の総積部分の尤度である。計算を易化するため,尤れる。C を 1 回のシミュレーションに必要な計算量とし度の対数を-1 倍した評価関数た と き , H-1 を 陽 に 求 め る た め に 必 要 な 計 算 量 はnJ=-∑ log r( y(ts)-h(x(ts))) ……………………(4)s=1O ( CN 2 + N 3 ) であり, N が大きい場合には現実的ではない。の最小化について考えることにする。式(1)からも分かそこで我々は, 2nd order adjoint 法11) を導入するこる通り, x ( ts ) は x0 = x ( 0 ) が与えられれば全て決定論とにより, H-1 を陽に求めることなく,不確実性評価的に得られるため, J の最小化にあたっては x0 の最適に必要な成分 Hll-1 のみを算出する手法を開発した10)。化を図れば良い。4 次元変分法は勾配法に基づいてこの2ndorder adjoint 法は,接線形モデル最適化を行うが,式(4)には x0 が陽に含まれていないため,勾配 &J / &x0 を直接求めることができない。そこで,式( 1 )の拘束の下で式( 4 )の最小化を図るため,ラグランジュ未定乗数 l を導入し,Tf{J¼ +lTJ=0(f(x)-&x&t)}dj &f= j ………………………………………………(11)dt &xj(0)=j0 ………………………………………………(12)及び 2ndorder adjoint モデルdt ……………………(5)( ) ( )dh&2f+jdt&x2Tl+&f&xTh+&2J¼j=0 ……………(13)&x2について考える。ただし, T は ts 以後の時刻, T は転h(T )=0 ………………………………………………(14)置を表し,J¼ は Dirac の d 関数を用いて,h(0)=Hj0 ……………………………………………(15)nの組み合わせによって構成される。章で述べた従来のs=14 次元変分法は,評価関数の勾配を計算するための技法J¼ =-∑ d(t-ts)log r( y(ts)-h(x(ts))) ……………(6)で与えられる。式(5)の変分をとり,停留条件を求めると,l に関する以下のアジョイント方程式が得られる。( )dl&f+dt&xTl+&J¼=0 ………………………………(7)&xであったが, 2nd order adjoint 法は,評価関数のヘッセ行列 H と任意のベクトル j0 の積を計算するための技法であることが分かる。目的の Hll-1 を得るために,次の線形方程式を反復法によって解くことを考える。Hq=dl ………………………………………………(16)&J………………………………………………(8)&x0ただし, dl は第 l 成分のみが 1,他の成分は 0 のベクトl(T )=0 ………………………………………………(9)ルを表す。これを 2nd order adjoint 法を用いて解くたl(0)=これらは,式( 9 )を初期条件として,式( 7 )を時間後退めの手順は,例えば勾配法に基づいて解く場合には,以方向に解くことにより, t = 0 のときの解として,目的 q の初期値を与える, 評価関数 J ′下の通りである。であった評価関数の勾配が得られることを示している。2/ &q = 2H=∥Hq - d∥l (∥・∥は L2 ノ ル ム ) の 勾 配 &J ′October, 20173論説( Hq - dl ) を, 2nd order adjoint 法を 2 回適用して計算 得られた勾配に基づき q を更新する, 停止する, と を繰り返す, 最終的に求めら条件に達するまでれた q の第 l 成分が目的の Hll-1 である。実際には,効率よく解を求めるため,我々は勾配法として共役残差法を採用している10) 。ここで特筆すべきは,この不確実性評価法の計算量がわずか O(C) であることであり,陽に H-1 を求める場合とは比較にならないほど,高速な計算が可能である。.図―小林フェーズフィールドモデルに基づく時間発展解フェーズフィールドモデルに基づく検証. 小林フェーズフィールドモデル章と章で述べた新しい 4 次元変分法を検証するための数値モデルとして,本稿では固体・液体間等の界面移動を記述するフェーズフィールドモデル( PFM )を用いる。PFM は系全体のエネルギーを定義した上で,界面を空間の連続場で表現し,そのダイナミクスを数値的に取り扱うモデルの総称である。 PFM は様々な物理現象を再現するために用いられ,金属材料の内部の構造相転移15) ,密度差の大きい流体のシミュレーション16),亀裂進展問題17) など多岐にわたって応用されている。PFM は自由度が高く,様々な現象を再現できるが,現象論的なモデルであるため,一般に多くのモデルパラメータを含む。また PFM は連続場のモデルであり,界面を精度よく記述するために空間離散化の際に格子点を図―小林フェーズフィールドモデルに基づく双子実験細かく取らなくてはならず,それに伴って状態の次元がによる( a )パラメータと( b )初期状態の同時推定,大きくなるため,アンサンブルベースの逐次型データ同及び( c )( d )パラメータ推定の不確実性のサンプリング間隔依存性化では計算量が劇的に増大する。したがって,提案手法の応用先として適したモデルの一つだと考えられる。まずは,単純な固液相転移問題を考え,固体相と液体ことにより,無制約最適化問題に転換することができ相が共存している状況で,固体相が成長する問題を考える10)。提案手法を適用した結果,図―に示す通り,mる。時刻 t,場所 x での固体相・液体相の存在確率をそと q(x, 0) の両者が反復を重ねるにしたがって真値に収れぞれ q(x, t), 1-q(x, t) とおく。2 つの相の競合的な束し,特に m については,観測データ量に応じた不確時間発展を表現する最も単純かつ基本的な小林 PFM12)実性 dmâ が適切に評価されていることが確認できる。観は,以下で与えられる。測データ量が十分な(サンプリング間隔 DT が十分短t()&q 2 211=e : q+q(1-q) q- +m , |m|< …(17)&t22ここで, :2 はラプラス演算子, t, e, m は定数であり,い)場合については dmâ ∝ DT の関係があり,大数の法則と整合している。. マルチフェーズフィールドモデル特に m は界面の移動速度を特徴づける重要なパラメー構造材料内における粒成長等,より現実的な系をタである。図―に, m = 0.1 の場合の二次元空間内にPFM によって取り扱うために,複数の相の存在を考慮おける小林 PFM に基づく時間発展解を示す。したマルチフェーズフィールドモデル(MPFM)提案手法の検証のための双子実験では,図―のような時間発展の写真が観測データとして得られた場合に,& qi=g&t(ne2:2qi+2qi2-2qi ∑ qj2j= 1)………………(19)提案手法によって推定されたパラメータ m と初期状態を導入する12)。図―は MPFM に基づき,提案手法をq(x, 0) が真値に収束することを示す。ただし,m が時双子実験によって検証した結果を示す6)。モビリティを間変化しないことを記述するモデル& m′=0,&t= m+m′1………………………………(18)2を導入し,式( 17 )とセットでシステムモデルを構成す表すパラメータ g 及び初期状態 q(x, 0) に対して真値を設定し,これらに基づくシミュレーションによって得られた擬似実験データ q (x, t) に提案手法を適用して g 及び q(x, 0) の再現性を調べたところ,それぞれ真値に収る。このように設定することにより,m と q(x, t) はと束することを確認した。特に g については,小林 PFMもに 0 と 1 の間で制約される量となる。このような制の時と同様,観測データ量に応じた不確実性が適切に評約付き変数の最適化については,簡単な対数変換を施す価されたことは特筆すべきであり,この情報を利用した4地盤工学会誌,―()論2)3)4)5)6)図―マルチフェーズフィールドモデルに基づく双子実験における( a )真の初期状態,( b )初期状態推定開始時の初期解,( c)推定された初期状態,及び7)( d )( e )パラメータ推定の不確実性のサンプリング間隔依存性8)実験デザイン最適化が期待できる。. お わ り に9)10)本稿では,様々な科学分野への応用が広がっているデータ同化手法の中でも,特に大規模モデルに基づくデータ同化の場合に威力を発揮する 4 次元変分法の原理,ならびにフェーズフィールドモデルへの応用例につ11)いて紹介した。本稿で紹介した不確実性評価が可能な新しい 4 次元12)変分法10) は,非常に汎用性の高いデータ同化手法であり,現在,我々は様々な実問題への応用を進めている。13)ただし,四次元変分法における最適化は単純な勾配法に基づいているため,現状では初期値依存性が極めて大きい手法である。初期値依存性を可能な限り排した四次元14)変分法へと高度化を図るために,我々はマルコフ連鎖モンテカルロ法( MCMC )の並列版とも言えるレプリカ15)交換モンテカルロ法18) の実装に取り組んでいる(図―(c))。このような四次元変分法が完成すれば,大規模数値モデルと大容量観測データの場合であっても,系の状態推定のみならず,不確実性評価から観測・実験計画16)の最適化に至るまで可能な,究極のデータ同化の方法論17)が実現するであろう。参1)考文献樋口知之・上野玄太・中野慎也・中村和幸・吉田 亮データ同化入門―次世代のシミュレーション技術―,朝倉書店,2011.October, 201718)説淡路敏之・蒲地政文・池田元美・石川洋一データ同化―観測・実験とモデルを融合するイノベーション―,京都大学学術出版会,2009.露木 義・川畑拓矢気象学におけるデータ同化,気象研究ノート,Vol. 217,日本気象学会,2008.Hoshiba, M. and Aoki, S.: Numerical shake predictionfor Earthquake Early Warning: Data assimilation, realtime shake mapping, and simulation of wave propagation, Bull. Seismol. Soc. Am., Vol. 105, No. 3, pp. 13241338, 2015.Niwayama, R., Nagao, H., Kitajima, T. S., Hufnagel, L.,Shinohara, K., Higuchi, T., Ishikawa, T. and Kimura, A.:Bayesian inference of forces causing cytoplasmic streaming in Caenorhabditis elegans embryos and mouse oocytes, PLoS ONE, Vol. 11, No. 7, doi:10.1371/journal.pone.0159917, 2016.Ito, S., Nagao, H., Kasuya, T., and Inoue, J.: Graingrowth prediction based on data assimilation by implementing 4DVar on multiphaseˆeld model, submitted to Sci. Technol. Adv. Mater.Miyoshi, T. and Kondo, K.: A multiscale localization approach to an ensemble Kalman ˆlter, SOLA, Vol. 9, pp170173, 2013.Evensen, G.: The ensemble Kalman ˆlter: Theoreticalformulation and practical implementation, Ocean dynamics, Vol. 53, No. 4, pp. 343367, 2003.Kitagawa, G.: Introduction to time series modeling, CRCPress, 2010.Ito, S., Nagao, H., Yamanaka, A., Tsukada, Y., Koyama,T., Kano, M. and Inoue, J.: Data assimilation for massiveautonomous systems based on a secondorder adjointmethod, Phys. Rev. E, Vol. 94, Issue 4, 043307,doi:10.1103/PhysRevE.94.043307, 2016.Le Dimet, F. X., Navon, I., Daescu, D. N.,: Secondorderinformation in data assimilation, Monthly Weather Review, Vol. 130, No. 3, pp. 629648, 2002.Kobayashi, R.: Modeling and numerical simulations ofdendritic crystal growth, Physica D, Vol. 63, Issues 34,pp. 410423, 1993.Steinbach, I., F. Pezzolla, B. Nestler, M. Seeßelberg, R.Prieler, G. J. Schmitz, and J. L. L. Rezende, A phase ˆeldconcept for multiphase systems, Physica D, Vol. 94, No.3, pp. 135147, 1996.久保田光一・伊理正夫アルゴリズムの自動微分と応用,コロナ社,1998.Shimokawabe, T., Aoki, T., Takaki, T., Endo, T.,Yamanaka, A., Maruyama, N., Nukada, A., andMatsuoka, S.: Proceedings of 2011 International Conference for High Performance Computing, Networking,Storage and Analysis, SC'11, pp. 3:13:11, 2011.Jacqmin, D.: Calculation of twophase NavierStokes‰ows using phaseˆeld modeling, J. Comput. Phys., Vol.155, Issue 1, pp. 96127, 1999.Aranson, I. S., Kalatsky, V. A. and Vinokur V. M.: Continuum ˆeld description of crack propagation, Phys. Res.Lett., Vol. 85, No. 1, pp. 118121, 2000.Hukushima, K. and Nemoto, K.: Exchange Monte Carlomethod and application to spin glass simulations, J. Phys.Soc. Jpn., Vol. 65, No. 6, pp. 16041608, 1996.(原稿受理2017.6.30)5
|