神戸大学統合研究拠点 「惑星科学研究センター」
平成24年6月7日(木)アンサンブルベース逐次データ同化とHPC
樋口知之

1/52

（統計数理研究所）

朝日新聞平成14年2月17日(日)2/52

14版第２社会面

ＪＲ東 安全研との共同研究

4/52

帰納の論理●論理学者 vs. 数学者 『数学者の厳密性の甘さかげんに我慢がならない』●物理学者：『６０という数はすべての数で割りきれる』
1,2,3,4,5,6 OK. “勝手にとったところの”10,12,15,20,30,...
●技術者：『すべての奇数は素数』

例を使って考えるなといっているのではない．
1,3,5,7 OK. 9? “ひとまず横において”，11,13 OK.
むしろそうしろと言っているのだ．
よって，“これは測定の誤差”

個々の特殊な事例を調べると，そこに共通する性質や関係から，一般的な命題や法則を推論できることがある．これを帰納という．帰納の結果が誤りに導くこともいかに多いかをわれわれは知っている．しかし，注意深く，感性を研ぎ澄ましてかかれば，帰納がすばらしい結果に導くことがあるということもわれわれは数々の例からしっているのである．
Ｇ．ポリア著，「帰納と類比」から（実際は，金出武雄著「素人のように考え，玄人として実行する」から）

5/52データ中心主義の時代へ

100Gb/day/sequencer
データ格納
データ処理

データ産出

6/52

内挿と外挿問題
データ

7/52－得意と不得意－データ無しの領域

CPSとは• Center for Planetary Science• Cyber Physical System

8/52

予測能力の向上に集中！
二つの機能の性能の合成
フォワード（前向き）計算モデルの記述力＋対象の現状態（現況）を捉える認識力＝ 予測能力計測手法のイノベーションと直結

地球や生命体のような複雑なシステムの理解と制御においては、対象に関する知識は常に不完全であることを前提に、現象の予測能力でもって研究の進め方を評価し、修正する方策が有効

9/52

つなぐ：データ同化

データ同化

The Fourth Paradigm:
Data-IntensiveScientific Discovery
By Tony Hey, Stewart Tansley,and Kristin Tolle, Eds.
Microsoft Research, 2009.

10/52

データ同化の目的：気象・海洋学の観点から[1] 予報を行うための最適な初期条件を求める。これは既に、現業の天気予報で実用化されていることである。
[2] シミュレーションモデルを構成する際の最適な境界条件を求める。連成現象を取り扱う際の適応的な境界条件設定もこの作業に含まれる。
[3] スケールが異なるシミュレーションモデル間の橋渡しを行うスキーム内に含まれる諸パラメータの最適な値を求める。経験的に与えられるモデル内のパラメータ値の検証も一つの具体例である。
[4] シミュレーション（物理）モデルにもとづいた、観測されていない時間・空間点における観測値の補間を行う。この作業は再解析データセットの生成とも呼ばれる。このデータセットから新しい科学的発見をもくろむ。
[5] 時間・経費を節約できる効率的な観測システムを構築するための仮想観測ネットワークシミュレーション実験や感度解析を行う。
（参考文献： 蒲地 他、「統計数理」、54(2), 223-245, 2006.)
1011/52

I SM

シミュレーションモデルの構成 (1)（Tsunami simulation model in Japan Sea）

x cx 2  t

PDE : Partial differential equation

xSupposet  F (a xt 1 )where we conduct a two-dimensionalcase(time varying)simulation experiment for understanding
the flow ofshallow water such as the tsunami.

vt

xt  F ( xt 1 , vt )12/52

physical variable vector  m,t isassigned at each grid point.

ξ m,tmm 1

ξ m 1,tシミュレーションモデルの構成 (2)Two-dimensional
water flow vector
Sea surface height which ismeasured from the averagesea surface
Depth is taken from thebottom topography data set.

In fact, this data set is known to be erroneous.

13/52

Vectorシステムモデルとしてのシミュレーションモデル
（simplified meteorological model around Japan）

State Vectorx cx 2  t

PDE : Partial differential equationxt  f t ( xt 1 )代入計算

xt  f t ( xt 1 , v t )14/52(time varying)vt ξ 1,t    ξ m ,t x t   ξ m  1,t    ξ M ,t  θデータ同化と一般状態空間モデル
State Vector（Simulation variables）

システムモデル

Stochastic simulation modelx t  f t ( x t  1 , v t ),v t ～ p ( v | θ sys )y t  h t ( x t , w t ),w t ～ p ( w | θ obs )気象・海洋のデータ
同化の枠組みtyt  Ht xt  wt , wt～N (0, Robs ) t : sampling time of observatio ns
t : simulation
time step
Missingvalue (Outlier) t  1  t

Observation model
Measurement model

観測モデル

t
time integration

15/52

条件付き分布

逐次ベイズ計算 --- 日次株価データを考えると --predictive density: p( xt | y1:t 1 ) きのうまでのデータに基づく今日の状態
予測分布

filterdensity:p(x|y)
t1:t
フィルタ分布p(x|y)
smootherdensity:
t1:
T
平滑化分布
jp( x j | y1i: )i

y1:T  { y1 ,  , y T }p( xt | y1:t 1 )filteringp( xt | y1t: )p( xt 1 | y1:t )p( xt 1 | y1:t 1 )３つの条件付分布と３つの操作

smoothing

16/52

数年後，データをすべて得たもとで振り返った今日の状態

時刻Tまでのデータをまとめたベクトル

predictionp( xt 1 | y1:t 1 )

今日までのデータに基づく今日の状態

p( xt | y1T: )p( xt 1 | y1:T )p( xT | y1T: )15

フィルタリングとベイズの定理
xt : 時刻 tの興味のある対象
yt : データ

x1:T  x1 ,  xT y1:T  y1 ,  , yT 観測モデル

フィルタリング

予測分布p( yt | xt ) p( xt | y1:t1)p( xt | y1:T ) x p( yt | xt ) p(xt | y1:t1)
フィルタ分布tT

データの尤度p( y1:T | )  p( yt | y1:t1,) 一期先尤度
t 1

パラメータの事後分布

事前分布p( y1:T |  ) p( )p( | y1:T )  p( y1:T |  ) p( )事後分布1617/52

ベイズの定理がなぜ今役立つのか？４つの理由
イギリスの牧師・数学者（1702 – 1761年）
1763年に発見

x : 興味のある対象

２．対象の特徴をとらえるセンサー性能の向上
高精度センサーのコモディティ（日用品）化
４．高速（無線）イン
ターネット網の整備

y : データ
ベイズの反転公式

p( x | y) p( y | x ) p( x ) p( y | x ) p( x )１．膨大な数の積分（和）操作には高速な計算機が必要
コンピュータの性能向上

３．対象の細かい情報を不確実性を含めて数値化。個人の情報を網羅的に収集ストレージの廉価化
1718/52

ベイズの定理と情報循環
p( y | x )  p( x )p( x | y)  p( y | x )  p( x )p( y)Posterior

Likelihood

Prior

Improved knowledge

Feasibility of realization of yfor given x

Beliefabout values of xxの空間
確率分布

19/52

尤度関数

about values of x

確率分布

逐次データ同化のアルゴリズム
逐次データ同化では観測を得るたびに確率変数の分布または値の推定を行う

xn（ p ( xi | y1:k )  p ( xi | y1 , y2 , , yk ) ）y n 1p ( xn 1 | y1:n 1 )p ( xn | y1:n 1 )yn

シミュレーションp ( xn | y1:n )シミュレーション

20/52

p ( xn 1 | y1:n )y n 1

逐次（アンサンブル）vs. 非逐次（最適パス）
データ同化のイメージ
データ空間にシミュレーション解を射影した時

ht ( x t )逐次（オンライン）型： 集団の時間発展を追う。つまり、Swarm Filter

非逐次（オフライン）型：ベスト初期値をもつパスを求める

代表例： EnKF (Ensemble KalmanFilter)

代表例： ４次元変分法（Adjoint法）

21/52

予測のステップ (EnKFとPFで共通)
State xf t ( xt(i )1|t 1 , v t( i ) )

x



N(i )
t 1|t 1 i 1xt(|it )1

x (i ) Nt |t 1 i 1xt(|t2)1xt(11) |t 1xt(21)|t 1
simulation: ensemble member of predictive PDFxt(|1t )1xt(N1|)t 1t 122/52xt(|tN1)

t: ensemble member of filtered PDF

Prediction stepTimePFでのフィルタリングのステップ
Statex: ensemble member of predictive PDF

xx ((ii)) N N
tt|t|t 1i i11 likelihoodx( 2)t |t 1: ensemble member of filtered PDF

t(|t2 )Observation:ytCalculate(i ) p ( yt | xt |t 1 ) likelihood for t(|ti )  ( j) each particle  p ( yt | xt |t 1 ) x

(1)t |t 1(N )t |t 1xt
23/52

j

Filtered by a resample proportional to likelihood

t(|t1)t(|tN )Filtering stepTimeFilterEnKFでのフィルタリングのステップ
statext

xx ((ii)) N N
tt|t|t 1i i11xt(|t2)1Linear Gaussian
observation modelyt  H t x t  w t: particle for predictive pdf: particle for filtered pdf

Sample-based Variance
Covariance Matrix :Vˆt|t 1

observation:xt(|1t )1xt(|tN)1yt

Kˆ t  Vˆt |t  1 H tT ( H tVˆt |t  1 H tT  R t )  1xt(|it)  xt(|it )1  Kˆ t ( yt  wt(i )  H t xt(|it )1 )Kalman Gain

Filtering

t
24/52

time
2323

選択と集中策時刻

t-1

p ( xt 1 | y1:t 1 )  X t 1|t 1

一つのシナリオvt(i ) システムノイズボール

Ft ()p ( xt | y1:t 1 )  X t |t 1yt

時刻t

選択と集中策による計算機資源の有効利用

クローン

× ×
時刻

25/52

t+1

××p ( xt | y1:t )  X t |t

４次元変分法とベイズ統計
システムモデル予測xˆt  f (xˆt1)

システムノイズを入れたくない（保存則が破れることを避けたい）

フィルタリングp( yt | y1:t 1)   p( yt | xt ) p( xt | y1:t 1)dxt p( yt | xˆt )Tp( y1:T | )  p( yt | y(1:kt)1,)

初期値のみ最適化p( xˆ | y1:T )  p(ty1 1:T | xˆ 0 ) p( xˆ | y:0 )
(k )
0(k )
0xˆ 0(k 1)  xˆ 0(k )パラメータの事後分布p( | y1:T )  p( y1:T |  ) p( )(k )( k 1)



(k )

(k )

Importance Sampling
MCMC
2530/52

4DVar vs. アンサンブルベースデータ同化法
逐次型

非逐次型
4DVar
力学的バランスを重視（システムノイズ無し）

Pros

・状態ベクトルの次元が非常に大きい、最大規模のシミュレーションモデ
ルが取り扱える。
・感度解析が可能（○）

・ベクトル計算向き

Cons

・時間を遡るシミュレーションコードを書きおろす必要があるため、人的労力の負荷が高い。

（×）

31/52

・統計モデルでないので、統一的な視点や基準でもってモデル解析ができない。

EnKF
簡便＆安心

・実装が容易PF

原型

超簡便＆原理的には万能・実装が著しく容易・退化現象（分布表現能力の減 ・観測モデルが非線形の場合少。）が原理的におきない。にも自然に対応可能・並列計算向き・共分散行列の更新ステップの計算コストが高い。

・時不変パラメータ推定問題の場合は、退化現象がおきる。

・分布がガウスから大きく逸脱し ・厳密な平滑化アルゴリズムのた時には誤った結果を導く。
実現が実質的に無理・超高次元状態ベクトルのシミュレーションモデルが取り扱えない。

Merging particle filter (MPF)

M個で重みつき平均M=3の場合

32/52

27

非線形フィルタリング間の関係

2833/52

並列計算機への粒子フィルタの実装• リサンプリングの際に，node間通信が頻発．
• 任意の2 nodes間でほぼat randomに通信が発生し，並列化しにくい．

34/52①階層構造をはじめから意識1ノード（CPU） ＝ ８コア
16GFlops x 8 = 128 Gflops
ノード数は８万以上
SPARC64 VIIIfx[第２階層]ノード内
OpenMP：スレッド並列[第１階層]
コア内：ＳＩＭＤ化（コンパイラが対応）[第３階層] ノード間 MPI：プロセス並列
ノード

ノード

ノード

ノード

ノード

30

35/52

①階層構造をはじめから意識プロセス並列：MPI
第３階層：
CPU間

ユーザが指定（頑張る）

スレッド並列：OpenMP
第２階層：
CPU内

SIMD：コンパイラ
第１階層：
コア内

3136/52

MetaPF: Meta-particle filter
Nodeに割り当てられたensemble の subset (=超粒子)をresamplingする．

現状では，nodeに割り当てられた重みがどれか1つでも0.3を超えたらresamplingすることにしている．
37/52

②ネットワーク構造をはじめから意識
６次元メッシュトーラス : 3次元の直方体に配置したノードをそれぞれ6方向で結合し、各次元がそれぞれリング状に結合されるネットワーク構成

Q. ネットワークについて、三次元トーラスで隣り合うノードどうしの通信が速いとのことですが、隣り合わないノード間の通信はその通信回数分遅くなるのでしょうか。
出典：理研のホームページから

A. 隣り合わないノード間の通信は、基本的にはバケツリレー方式で行われるので、経由するノードの数が多ければ多いほど、通信時間（レイテンシ）は長くなります。
IBM Sequoia (BlueGene/Q)
Cray XT5(Jaguar) も３次元トーラス
33

38/52

Alternate lattice-pattern switchingノード

ここでは，上図のような2種類の格子パターンでノードをグループ化することを考え，2種類のパターンを毎回切り替えるものとする．

39/52

データ同化研究開発センターでの研究プロジェクト
・Coupled Ocean-Atmosphere model・Tsunami model・Ocean tide・ 3D structure of ring current
・Genome informatics
・Marketing (with multi-agent simulations)

40/52

シミュレーション (同化結果でない)• Simulation of Okushiri Tsunami– Simulation based on topographies made by different
organizations.
– It looks similar，but time series of sea surface
displacement at a point ( ) is ...

41/52(From Japan Coast Guard)海水面変位DBDBVSKKUにおける海水面変位（cm）
20
100-10-20102030

津波到着直前からの時間（分）

42/52

40

大きい相対的誤差：深さに関する４つのデータベース間の比較

Depth (m)43/52

潮位計設置点

観測データ
潮位計設置点• 潮位計データ– 設置地点付近の海面変位を反映– 得られるデータは１次元時系列• 設置地点は右図– 30点程度例）
北海道南西沖 40
地震津波20（韓国・束草）
潮位（cm）

0海上保安庁-20-40

44/52

050

100

150200250

津波到着直前からの時間（分）同化結果
浅深•大和海嶺周辺は４種類の海底地形データの平均よりもやや浅いと推定される•南斜面に平均よりも深いと推定される部分がある

45/52

“個”にマッチしたシミュレーション：
境界条件の設定機能をパーソナライズする“個”によって異なる形状，形態情報をシミュレーションモデルに取り込む 『メタシミュレーションモデル』
運動方程式：

連続式：

我々研究チームによる，潮汐シミュレーションの例

海底摩擦項

vv
v
 ( v   ) v  f  v   g    b
 AH  2 v
t
H
海底摩擦係数水深（地域依存性）   vH   0
tv： 水平（２次元）流速ベクトル
η： 海面水位
H： 水深，f： コリオリパラメータ

多品種少量生産を基本とする製品開発現場でのステップの簡略化や、患者一人一人に合った治療サービスの提供と期間の短縮化
アラスカ州
南東域

46/52

流動の原動力のパラメータ推定
国立遺伝学研究所・木村研究室と共同研究

細胞質流動（線虫）
・境界条件の推定： ミオシンがつくる剪断応力値を推定

細胞質流動の定量データ

ＰＩＶ（Particle Image Velocitmetry）計測

細胞質流動のシミュレーション

数値解法はＭＰＳ（Moving Particle semi-implicit）法を採用。連続体を多数の粒子で近似．計算格子を生成しなくていい．自由表面などの移動境界がある問題の扱いが簡単

47/52

東北地方太平洋沖地震（東日本大震災）に伴う地震波・地震音波伝搬シミュレーション
物理モデル• 地球中心から地表までの固体地球および地表から高度1000kmまでの大気を考慮した1次元地球構造モデル• 気象庁が決定した震源解（モーメントテンソル解）を長さ各計算グリッドにおける地震波（表面波）および音波の応答波形を150kmの断層に沿って配置した断層破壊モデル
ノーマルモード法（Kobayashi[2007]）によって計算

グリッド幅は緯度・経度方向0.1度および高度方向10km
断層を含む緯度・経度方向10度および高度方向120kmの範囲の可視化した結果を示す．
統計数理研究所が所有するスーパーコンピュータ（富士通Dシステム）を約4000ノード時間（512並列で約60時間）使用

津波よりも伝搬速度が速い音波を利用した将来型津波警報システムの構築
48/52

オーダーメイド医療・創薬に向けたPersonalized Simulation
東京大学・医科学研究所・宮野研究室と共同研究

データ同化
••

生体システムを統合的に理解・予測するための有望な方法論
実験科学と統計・シミュレーション科
学の融合

データ同化
Validation
DATAExperiments

10万コアの超並列計算に成功！

Parameter Estimation
Model Selection
Super Computer

Hypothesis
Design Creation

49/52 本スライドは東京大学（医科研）・山口類氏から借用，改変

Order-Made
Medicine,
Drug discovery
with
PersonalizedSimulation

Ａさん

Drug A

Ｂさん

Drug B“ものづくり”分野：

地盤シミュレーションのデータ同化
京都大学（農学研究科）・村上教授グループとの共同研究■実学的観点から： 地盤沈下の最終沈下量の予測・関西国際空港の埋め立てにおける沈下量予測など、途中までの計測を用いて最終沈下量の予測をしたい。
・データ同化の観点からは，シミュレーション予報のための初期値生成作業。
■計算手法の観点から： 弾性体であれば，線形性・復元性によりデータ同化は簡単。（カルマンフィルタ(KF)が有効。 実際の土は塑性をもつので、 変形の経路依存性を考慮しなければならなくなるので、ＫＦは適用不可。

載荷構造シミュレーション：

観測点

有限要素法、境界要素法

８点の観測（沈下量測定）を用いて，弾塑性パラメータを推定する(Wikipediaから）載荷過程

Embankment loading history

5・塑性変形も加わると，非線形性に伴うシミュレーション計算の破綻（保存則を満たさない等）の問題が出てきてしまい，KF では難しかったがSIS（逐次データ同化の一種） で解決できた。

Height （m）

432100100200300

Elapsed time (day)50/52400500(理研のホームページから）富を産む仕組み
前世紀： 物質（「もの」）を均質に大量に生産するシステム
21世紀： 個人化された情報サービスを提供するシステム
個人をターゲットにした商品・サービスの提供を効率的に行えるシステム“コ”－個人，個性，個別，固有ーが大切！

51/52

ご静聴ありがとうございました。
バーチャルにリアルを埋め込むデータ同化：
ベイズ、アンサンブル予測
ヒューマンモデリング、個人化情報サービス(2011年4月刊行）(2011年9月刊行）(2011年11月刊行）

データ同化研究開発センターの紹介ビデオ（１０分）を作成しました。YouTubeにアップしています。
52/52