粒子フィルタの実装
中野 慎也 (統計数理研究所)1
導入
粒子フィルタの特徴: • 確率密度分布を多数のサンプルで近似 (モンテカルロ近似) する手法の一つ． • ブートストラップフィルタ，モンテカルロフィルタ，Sampling/Importance resampling (SIR)フィルタなどと呼ばれることもある． • Kitagawa (1993, 1996), Gordon(1993) によって独立に提案された． • EnKF と同様，確率密度分布を多数のサンプルで近似する． • EnKF が基本的に線形観測モデルを仮定しているのに対し，より一般的な時系列モ デルに適用可能． • 考え方が素直で，実装も容易． • 並列計算への書き換えも容易． • ただし，EnKF と比較してかなり多数の粒子を必要とする傾向がある． なお，以下では 1CPU あるいは 1core に複数の粒子を割り当てるられる程度の規模のあ まり大規模でない問題を想定して話を進める．2
粒子フィルタの復習
想定するのは，非線形状態空間モデル xt= ft(xt−1, vt) (1a) yt= ht(xt, wt) (1b) である．このようなモデルの下，PF では，各ステップの予測分布 p(xt|yt−1)を近似する アンサンブル{x(i)t|t−1}N i=1,フィルタ分布 p(xt|yt)を近似するアンサンブル{x (i) t|t} N i=1を以下 のようにして逐次的に求める:State x State x ᄫഐက࿑ Ӈҕӝұက࿑ ჰ๗ 図 1: 粒子フィルタの手続きの概念図． 1. 初期分布のアンサンブル{x(i)0|0}Ni=1 (x (i) 0|0 ∼ p(x0) )を生成する． 2. t = 1,· · · , T について a から d までのステップを実行する． (a) 各 i (i = 1,· · · , N) について I, II を実行する． I. システムノイズを表現する乱数 v(i)t ∼ p(vt) を生成する．
II. x(i)t|t−1 = ft(x(i)t−1|t−1, v(i)t )を計算する． (b) 各 i (i = 1,· · · , N) について λ(i)t = p(yt|x
(i)
t|t−1) を計算する．
(c) βt(i) = λ(i)t /(∑Ni=1λ(i)t ) を求める．
(d) アンサンブル{x(1)t|t−1,· · · , xt(N )|t−1} から各粒子 x(i)t|t−1が βt(i)の割合で抽出される
よう復元抽出し，{x(1)t|t,· · · , x(N )t|t } を生成する．
このうち，2a が一期先予測を行う手続きで，2b から 2d までがフィルタの操作である．フィ ルタの操作は，具体的に言えば尤度の低い粒子を破棄して代わりに尤度の高い粒子の複 製を増やすという手続きに対応する (図 1)．
3
非並列計算機における粒子フィルタ
3.1
まずやるべきこと
粒子フィルタを実装しようとする際にまずやるべきことは，シミュレーションコードで 使われているすべての変数を， • 各粒子がそれぞれ固有の値を持つ変数と • そうでない変数と に分類することである． 粒子フィルタでは，各粒子がそれぞれ別々の初期条件，境界条件，パラメータ値などを 与えた run に対応しており，それぞれの run を並行して走らせることになる．ここで，粒 子間で値が異なる変数を粒子ごとに別々のメモリ領域に値を格納しておかないと，複数 の粒子 (run) を走らせた場合に各 run が正しく動かなくなってしまう可能性がある． 各粒子が固有の値を持つ変数とは，要するに状態空間モデルにおける状態変数そのも のである．つまり，ここで言っているのは，シミュレーション変数の中でどれが状態変数 になるのかを見極める必要があるということである．世の中のシミュレーションコード は，どれが状態変数であるかなどということを意識せずに書かれているのが普通であり， 意外にその選別が手間のかかる作業になることがある． どれを状態変数とするか 1. 全粒子で共通の値を取るような変数は，各粒子ごとに別々のメモリ領域を用意する 必要はなく，全粒子で共有するようにすれば充分である．(MPI で並列化する場合 で言えば，粒子ごとではなくプロセスごとにメモリ領域を用意すればよいというこ とになる．) 2. 同じシミュレーションモデルを使った場合でも，どの変数に不確定性を仮定するか によって，状態変数にすべき変数が変わってくる．例えば，境界条件を推定する場合には境界条件は状態変数になるが，境界条件をある決まった値で仮定したときに は，全粒子で共通の値を取るから状態変数として扱う必要はない． 3. 粒子ごとに異なる値を取ると言っても，ある特定のサブルーチンの中で一時的に使 われるだけで使い捨てにされるローカル変数については，粒子固有の変数に含める 必要はない． なお，オブジェクト指向プログラミングの言葉を借りれば，状態変数は各粒子をシミュ レーションモデル (クラス) のインスタンスと見たときのインスタンス変数に対応すると 見ることができる．オブジェクト指向プログラミング言語を使ってデータ同化のプログ ラムを書くならば，実際にそのような設計にした方がその後のコーディングは楽になる のではないかと思われる． 筆者の場合，Fortran 90 というオブジェクト指向言語ではない言語でプログラムを書い ているが，シミュレーションモデルの状態変数をすべてまとめて，一つの構造体—Fortran 90の用語で言えば derived type—にまとめるということはしている．このようにすると， プログラム上で各粒子をオブジェクトとして認識しやすくなるだけでなく， • どの変数が状態変数扱いなのか見通しがよくなる • 後述するリサンプリングの手続きがコーディングしやすくなる などの利点がある．これはあくまでも筆者の趣味ではあるが．
3.2
初期分布の生成
初期分布のアンサンブルの生成は，上の記述に従えば，分布 p(x0)からサンプルを N 個取ればよいということになる．ただ，問題は p(x0)としてどのような形のものを仮定す るかである．x0の各成分ごとに独立な正規乱数で生成する (これは対角の分散共分散行列 を仮定したことになる) のは，楽な方法ではあるが物理的に不自然な状態が生成されてう まくいかないことが多い．そこで，例えば， 1. シミュレーション上の隣り合うグリッドでは強い相関を持つ分散共分散行列を考える，2. 状態変数のごく一部にだけノイズを与え，あとはシミュレーションモデル自体でノ イズの影響を周りに行き渡らせることにより，ばらつきを持つ初期アンサンブルを 生成する などの手段によって初期分布を生成する．
3.3
システムノイズの生成
初期アンサンブルを各粒子の値がばらけるようにきちんと生成していれば，初回は必 ずしもシステムノイズを加える必要はないが，次のステップ以降はシステムノイズをき ちんと与える必要がある．ここでも，p(vt)としてどのような形のものを仮定するかが問 題になるが，状態変数のごく一部にだけ与えるなどの形を取ることが多いようである．3.4
予測
予測は単にあるステップにおけるある粒子の値 x(i)t−1|t−1 (あるいはそれにシステムノイ ズを加えたもの) を初期値にして，シミュレーションを時刻 t まで走らせるという手続き によって行う．したがって，この部分はシミュレーションが用意されていれば，基本的に はそれをそのまま流用するだけで済む．3.5
フィルタリング
フィルタリングの手続きは， 1. 各粒子の尤度を計算する． 2. 尤度に基づいて，各粒子の重みを計算する． 3. 各粒子の重みにしたがって，どの粒子の複製をどの粒子と置き換えるかについての 対応テーブルを作る． 4. 対応テーブルに沿って，実際に粒子の分配操作を行う． の 4 つの段階に分けられる．3.5.1 尤度の計算 状態ベクトル xtの尤度 p(yt|xt)の与え方は観測モデルによって決まる．例えば， yt= Htxt+ wt のような線形の観測モデルを仮定し，観測ノイズ wtが共分散 Rtの正規分布に従うと仮 定するならば，粒子 x(i)t|t−1の尤度は λ(i)t = p ( yt|x(i)t|t−1 ) = √ 1 (2π)l|R t| exp [ −1 2 ( yt− Htx (i) t|t−1 )′ R−1t ( yt− Htx (i) t|t−1 )] となる．その他，観測モデルとして Poisson 分布などの分布を仮定した場合などでも，し かるべき形で求めることができる． 3.5.2 重みの計算 重み βt(i)は，各粒子の尤度から計算される．ガウス分布を仮定した場合などには，まず 対数尤度を求めたうえで尤度を計算することが多いが，データが高次元の問題では対数 尤度の絶対値が大きな負の数になるため，コンピュータで尤度を計算しようとしたとき にアンダーフローがしばしば起こる．場合によっては，ほとんどの粒子でアンダーフロー が起こってしまうこともあり，そうなると精度が著しく落ちてしまうので注意する必要が ある．ただし，実際に PF の重みの計算で必要なのは尤度そのものではなく尤度の比であ ることを考慮すれば，以下のようにしてアンダーフローの影響を回避できる． 1. 全粒子のうち，対数尤度 ℓ(i)が最も大きい粒子を選び，その粒子を x(K)，その対数 尤度を ℓ(K)とする．
2. 各粒子ごとに，ψ(i) = exp(ℓ(i)− ℓ(K))の値を計算する．
3. Ψ =∑iψ(i)を計算する．
4. β(i) = ψ(i)/Ψで重みが得られる．
実際の重みの計算は，例えば以下の subroutine のようにすればできる．この subroutine では，各粒子の対数尤度の値を持つ配列を xlamarr とし，その対数尤度の値から重み weightを求める．
プログラム 1: 重み計算ルーチン s u b r o u t i n e c a l c _ w e i g h t ( nptotal , xlamarr , w e i g h t ) i m p l i c i t n o n e re a l (8) , d i m e n s i o n ( n p t o t a l ) :: xlamarr , w e i g h t re a l (8) :: xlammax , w s u m integer , i n t e n t ( in ) :: n p t o t a l i n t e g e r :: i x l a m m a x = m a x v a l ( x l a m a r r ) ! 対数尤度 x l a m a r r ( i ) の最大値を求める ! 尤度比 (= 重み ) の計算 w e i g h t ( 1 ) = exp ( x l a m a r r (1) - x l a m m a x ) ws u m = w e i g h t (1) do i =2 , n p t o t a l w e i g h t ( i )= exp ( x l a m a r r ( i ) - x l a m m a x ) w s u m = w s u m + w e i g h t ( i ) end do ! 和が 1 になるように規格化 do i =1 , n p t o t a l w e i g h t ( i )= w e i g h t ( i )/ w s u m end do r e t u r n end s u b r o u t i n e c a l c _ w e i g h t 3.5.3 リサンプリングの対応テーブル 重みが決まったら，どの粒子の複製をどの粒子と置き換えるかについての対応テーブ ルを作るわけであるが，テーブルを作る前に，予測アンサンブルの粒子をそれぞれ何個 抽出 (あるいは複製) するかをまず決める．抽出する粒子の個数の決め方には以下のよう な方法がある． 最も素直な方法 予測アンサンブル{x(i)t|t−1} から各粒子が βt(i)の確率でランダム抽出され るものとして，復元抽出を N 回繰り返せば，フィルタ分布を近似するアンサンブルが得 られる．これが，恐らく最も素直なやり方． しかし，ランダム抽出を行った場合，たとえすべての粒子の重み βt(i)が同じであったと しても，各粒子が抽出される回数にばらつきが出てしまうことになり，分布の近似精度は 必ずしもよいとは言えない．
ランダム抽出を N 回行ったときに，ある粒子 x(i)t|t−1のが抽出される回数がどのくらい ばらつくかは，粒子 x(i)t|t−1の抽出が確率 βt(i)のベルヌーイ試行と見なせることから直ちに 見積もることができる．(二項分布の標準偏差の式からすぐに出る．) 見積もってみると √ N βt(i)(1− βt(i))程度となる．したがって，フィルタアンサンブル予測アンサンブルのメ ンバー x(i)t|t−1の複製が{x(i)t|t}N i=1} に占める割合は √ βt(i)(1− βt(i))/N 程度のばらつきを持 つことになる．このことから，特に N があまり大きくない場合には，ランダム抽出を用 いると，得られたアンサンブルによって表現される分布が本来の分布の形状とは大きく 異なったものになってしまう可能性がある． 比例配分 1 そこで，抽出される粒子数の比が尤度比になるべく近くなるようにするために，ラン ダム抽出を行わず，代わりに以下のような方法を取ることがある．
1. N βt(i)を整数部分 ¯m(i)t として，まず，各粒子 x(i)t|t−1をそれぞれ ¯m(i)t 個ずつ抽出する． これで，∑im¯(i)t 個が確定する． 2. 粒子数を N 個にするために，各粒子 x(i)t|t−1が (N βt(i)− ¯m (i) t )/(N− ∑ im¯ (i) t )の確率で 抽出されるようなランダム抽出によって，不足分 N−∑im¯(i)t 個の粒子を抽出する． 比例配分 2 また，次のようにまったく乱数を使わずに決める方法も一つの手段として考えられる (図 2)． 1. まず，{ζt(i)}N i=0, {η (i) t }Ni=1をそれぞれ ζ(0) = 0, ζ (i) = i ∑ j=1 βt(j), ηt(i) = (i − ε)/N (0 < ε≤ 1) のように定義する．
2. i = 1,· · · , N について ζt(i−1) ≤ ηt(j)< ζ(i)を満たす j を探し，x(j)t|t := x(i)t|t−1とする． 抽出する粒子の個数が決まったら，元々予測アンサンブルの粒子がいた席にリサンプ リングによって得られるフィルタアンサンブルの粒子をどのように配置するかを定めた テーブルを作る．
0 ࠖџဖѭѓ ѳѢോсщћ ࣤಂѿ =1 (1) t ζ (2) t ζ (3) t ζ (4) t ζ (5) t ζ (N) t ζ (2) t η (1) t η (N) t η (6) t η (5) t η (4) t η (3) t η 図 2: 乱数を使わなくても抽出される粒子の数の比が，尤度比に近くなるようにすること ができる．
Re
sa
mp
lin
g
図 3: 破棄された粒子のいた席に，増やした粒子を置く． 1. 抽出される個数が 0 の粒子，すなわち破棄される粒子のインデックスをリストアッ プする． 2. 抽出される個数が 2 以上の粒子，すなわちリサンプリングによって数が増える粒子 について，その複製を破棄する粒子の席に順次配置していく． なお，細かいことではあるが，このようにすると破棄 されない 粒子がいた席はそのまま 元々いた粒子 (の複製) が占めることになる (図 3)．破棄しない粒子の席に別の粒子が配さ れる可能性があるようなテーブルの作り方をしたとすると，元いた粒子の値を一旦に退 避させるための一時的な変数を用意する必要が生じるが，このようにしておけば余計な メモリ領域を必要としない． 以下は，リサンプリングのテーブルを作る subroutine の一例である．この subroutine では，抽出する粒子の個数を「比例配分 2」の方法で決めることにより，リサンプリングのテーブル isrc を生成する． プログラム 2: リサンプリングのテーブル生成ルーチン s u b r o u t i n e r e s a m p l e _ m a p ( np , weight , i s r c ) i m p l i c i t n o n e integer , i n t e n t ( in ) :: np i n t e g e r :: i , j , n z e r o r e a l (8) , p a r a m e t e r :: eps = 0 . 5 r e a l (8) , d i m e n s i o n ( np ) , i n t e n t ( in ) :: w e i g h t r e a l (8) , d i m e n s i o n ( np ) :: zeta , eta integer , d i m e n s i o n ( np ) :: isrc , n s a m p l e integer , d i m e n s i o n ( np ) :: i z e r o ! ! 重みから ze t a を計算 z e t a ( 1 ) = w e i g h t (1) do i =2 , np -1 ze t a ( i )= w e i g h t ( i )+ z e t a ( i -1) end do z e t a ( np ) = 1 . 0 ! ! eta を計算 do i =1 , np eta ( i )=( i - eps )/ np end do n s a m p l e ( : ) = 0 ! i 番目の粒子が複製される数 i z e r o ( : ) = 0 ! 破棄される粒子の番号を順次参照する n z e r o =0 j =1 do i =1 , np if ( j <= np ) t h e n ! ! z e t a ( i -1) - - z e t a ( i ) の間の粒子数を数える do w h i l e ( eta ( j ) < z e t a ( i )) n s a m p l e ( i )= n s a m p l e ( i )+1 j = j +1 if ( j > np ) e x i t end do end if ! ! 破棄される粒子へのラベル付け if ( n s a m p l e ( i ) = = 0 ) t h e n n z e r o = n z e r o +1 i z e r o ( n z e r o )= i end if end do do i =1 , np is r c ( i )= i ! どの粒子をどこに移すかのテーブル end do
! ! 破棄される粒子の代わりに持ってくる粒子への対応付け j =1 do i =1 , np do w h i l e ( n s a m p l e ( i ) > 1) ! 粒子の複製が余っていたら i s r c ( i z e r o ( j ))= i ! 破棄される粒子の場所を余剰で埋める j = j +1 n s a m p l e ( i )= n s a m p l e ( i ) -1 end do end do r e t u r n end s u b r o u t i n e r e s a m p l e _ m a p 3.5.4 粒子の再分配 上述のリサンプリングテーブルにしたがって，尤度の低い粒子を尤度の高い粒子の複 製に置き換えていく．ここで先に述べたように状態変数をすべて構造体にまとめておく と便利である．例えば，状態変数を m o d u l e s y s t e m _ m o d e l integer , p a r a m e t e r :: n d i m =50 type , p u b l i c :: s t a t e _ v e c t o r re a l (8) , d i m e n s i o n ( n d i m ) :: a re a l (8) , d i m e n s i o n ( n d i m ) :: p end t y p e s t a t e _ v e c t o r end m o d u l e s y s t e m _ m o d e l
というように，‘state vector’ という derived type で固めておくと，
t y p e ( s t a t e _ v e c t o r ) :: x ( np ) r e a l (8) , d i m e n s i o n ( np ) :: w e i g h t ! ! 中略 c a l l r e s a m p l e _ m a p ( np , weight , i s r c ) do i =1 , np if ( i /= i s r c ( i )) t h e n x ( i )= x ( i s r c ( i )) end if end do のような操作で，リサンプリングが済んでしまうので，プログラムを書くのが多少楽に なる．なおこの例では，リサンプリングテーブル isrc が破棄しない粒子のいた席には必 ず元々いた粒子が配置されるように作られているものと仮定している．
Rank 0 ↝ ⇽∓⇡⇟ Rank 1 ↝ ⇽∓⇡⇟ Rank 2 ↝ ⇽∓⇡⇟ Rank 3 ↝ ⇽∓⇡⇟ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ ᡫ̮ 図 4: MPI による並列計算の概念図．
4
MPI
による並列計算
並列計算を行う手段には幾つかあるが，最も普及しているのは MPI (Message Passing Interface)と呼ばれるものであろう．MPI を使った並列計算プログラムは，走らせると複 数のプロセスが並行して走るようになっている．但し，プロセスはそれぞれ独立で走る のではなく，互いに通信して情報を交換し合いながら走る (図 4)．MPI を使うときに新た に意識すべき部分は，このプロセス間の通信の部分であり，MPI ではこの通信のための subroutineが多数提供されている．MPI はライブラリの形で提供されており，ユーザー は MPI の subroutine をプログラム上で call することでプロセス間通信を指示することに なる．MPI のライブラリには，Fortran 用のものも C 用のものも提供されているが，ここ では特に Fortran の場合について説明をする．C の場合も使用する関数などは Fortran と ほとんど同じである． なお，MPI を使って書いたプログラムをどうコンパイルするか，そしてどう実行するか について，実行したい環境に依存してしまうのでここで述べることはしない．また，MPI には多種多様な subroutine が用意されているが，全部説明するのは大変なので，以下で は，粒子フィルタを用いたプログラムを書く際に必要な知識に絞って話を進める．MPI
についての詳細について知りたい場合には，理化学研究所の青山幸也先生の書かれた「並 列プログラミング入門 MPI 版」を読むとよいだろう．このテキストは，Web1) から無料 でダウンロードできる．印刷するのが面倒な場合は，理化学研究所で半年に一度くらい の頻度で開催されている講習会に参加すると製本されたものが手に入る．
4.1
MPI
の基礎
まず，MPI を使うにあたって最低限必要そうな知識を述べておく．以下は MPI のごく 基本的な subroutine だけを使って書いたプログラムである． プログラム 3: MPI のプロセス基本情報の取得 p r o g r a m M P I p r o g r a m i m p l i c i t n o n e i n c l u d e ’ m p i f . h ’ i n t e g e r :: nprocs , i r a n k i n t e g e r :: i e r r c a l l m p i _ i n i t ( i e r r ) c a l l m p i _ c o m m _ s i z e ( M P I _ C O M M _ W O R L D , nprocs , i e r r ) c a l l m p i _ c o m m _ r a n k ( M P I _ C O M M _ W O R L D , irank , i e r r ) p r i n t * , ’ Num ␣ of ␣ p r o c : ’ , nprocs , ’ ␣ R a n k : ’ , i r a n k c a l l m p i _ f i n a l i z e ( i e r r ) s t o p end p r o g r a m M P I p r o g r a m include ’mpif.h’MPIの subroutine を用いるプログラム，モジュールでは，必ず mpif.h というファ イルを include する必要がある．
call mpi init(ierr)
他の MPI の subroutine を使う前に一度だけ call する必要がある．MPI を使わない 処理ならばこの文の前に書いてあっても差し支えない．ierr にはこの subroutine が 正常に終了すると 0 が返る．
1)
call mpi comm size(MPI COMM WORLD,nprocs,ierr)
実行したときに同時に走るプロセスの数を取得するときに使う．nprocs にプロセス 数が返る．MPI COMM WORLD は mpif.h で定義されている定数だが，とりあえ ずこう書くものだと思っておいてよい．
call mpi comm rank(MPI COMM WORLD,irank,ierr)
Rankとは各プロセスに振られる通し番号のようなもので，0, 1, 2, . . . というように 順次振られていく．irank に rank の値が返る．
call mpi finalize(ierr)
MPIの subroutine を使い終わった後に必ず call する必要がある．この文の後に MPI に関する処理を書くことはできないが，MPI を使わない処理をこの文の後に書くの は問題ない． このプログラムをプロセス数 8 で実行すると， Num of p r o c : 8 R a n k : 0 Num of p r o c : 8 R a n k : 1 Num of p r o c : 8 R a n k : 2 Num of p r o c : 8 R a n k : 3 Num of p r o c : 8 R a n k : 5 Num of p r o c : 8 R a n k : 7 Num of p r o c : 8 R a n k : 4 Num of p r o c : 8 R a n k : 6 といったような結果が出る．
4.2
データの収集
次に，データを全プロセスからある特定の一つのプロセスに集める mpi gather という subroutineを紹介する．集められたデータは，各プロセスから受け取ったデータをすべて 並べた配列に格納される (図 5)． 以下のプログラムは mpi gather を使った例である． プログラム 4: mpi gather 使用例 p r o g r a m M P I g a t h T e s t i m p l i c i t n o n e i n c l u d e ’ m p i f . h ’Rank 0 ↝ ⇽∓⇡⇟ Rank 1 ↝ ⇽∓⇡⇟ Rank 2 ↝ ⇽∓⇡⇟ Rank 3 ↝ ⇽∓⇡⇟
MPI_GATHER
図 5: MPI gather を使うとデータは一つのプロセスに集められる．集められたデータは， 各プロセスから受け取ったデータをすべて並べた配列に格納される． i n t e g e r :: nprocs , irank , i r o o t i n t e g e r :: ierr , i s t a t i n t e g e r :: nn i n t e g e r :: i integer , p a r a m e t e r :: n =2 r e a l (8) :: x ( n ) r e a l (8) , a l l o c a t a b l e , d i m e n s i o n (:) :: z i r o o t =0 c a l l m p i _ i n i t ( i e r r ) c a l l m p i _ c o m m _ s i z e ( M P I _ C O M M _ W O R L D , nprocs , i e r r ) c a l l m p i _ c o m m _ r a n k ( M P I _ C O M M _ W O R L D , irank , i e r r ) do i =1 , n x ( i )= n * i r a n k + i end do nn = n * n p r o c s a l l o c a t e ( z (1: nn ) , s t a t = i s t a t ) c a l l m p i _ g a t h e r ( x , n , M P I _ R E A L 8 , z , n , M P I _ R E A L 8 , iroot ,& M P I _ C O M M _ W O R L D , i e r r ) if ( i r a n k == i r o o t ) t h e n do i =1 , nn w r i t e (6 , ’ ( a3 , i3 ,2 x , a6 , f5 .1) ’ ) ’ i : ␣ ’ ,i , ’ z ( i ): ␣ ’ , z ( i ) end doend if
d e a l l o c a t e ( z , s t a t = i s t a t ) c a l l m p i _ f i n a l i z e ( i e r r ) s t o p
end p r o g r a m M P I g a t h T e s t
call mpi gather(x,n,MPI REAL8,z,n,MPI REAL8,iroot,. . . 各プロセスのデータを rank が iroot のプロセスに集める．
• 1 つ目の引数 x: 各プロセスから送信されるデータの入った配列． • 2 つ目の引数 n: 送信されるデータ配列の長さ．
• 3 つ目の引数 MPI REAL8: 送信されるデータの型が 8 バイト実数であることを 示す．(環境によっては MPI DOUBLE PRECISION にする必要があるかもしれな い．) 値が MPI REAL4 なら 4 バイト実数，MPI INTEGER なら整数型のデータを 送ることになる． • 4 つ目の引数 z: データを受信するための配列．配列の長さは 2 つ目の引数の 値にプロセス数を掛けた数． • 5 つ目の引数 n: 1 プロセスから送られてくるデータの長さ．普通は 2 つ目の引 数と同じにする． • 6 つ目の引数 MPI REAL8: 受信するデータの型．普通は 3 つ目の引数と同じに する． • 7 つ目の引数 iroot: 受信プロセス (データを集める側) の rank の番号を入れる．
4.3
データの同期
プロセスのある変数 (配列) の値を特定のプロセスにおける値で同期させる mpi bcast という subroutine も使われる．以下は，mpi bcast を使った例である．プログラム 5: mpi bcast 使用例
p r o g r a m M P I b r o a d c a s t i m p l i c i t n o n e
i n c l u d e ’ m p i f . h ’ i n t e g e r :: nprocs , irank , i r o o t i n t e g e r :: i e r r i n t e g e r :: i integer , p a r a m e t e r :: n =5 r e a l (8) :: x ( n ) i r o o t =0 c a l l m p i _ i n i t ( i e r r ) c a l l m p i _ c o m m _ s i z e ( M P I _ C O M M _ W O R L D , nprocs , i e r r ) c a l l m p i _ c o m m _ r a n k ( M P I _ C O M M _ W O R L D , irank , i e r r ) if ( i r a n k == i r o o t ) t h e n do i =1 , n x ( i )= i end do e l s e x ( : ) = - 1 . 0 end if
w r i t e (6 , ’ ( a16 , i3 ,2 x , a3 ,5 f6 .1) ’ ) ’ B e f o r e ... ␣ R a n k : ␣ ’ , irank ,& ’ x : ␣ ’ ,( x ( i ) , i =1 , n )
c a l l m p i _ b c a s t ( x , n , M P I _ R E A L 8 , iroot , M P I _ C O M M _ W O R L D , i e r r ) w r i t e (6 , ’ ( a16 , i3 ,2 x , a3 ,5 f6 .1) ’ ) ’ A f t e r ... ␣ ␣ R a n k : ␣ ’ , irank ,&
’ x : ␣ ’ ,( x ( i ) , i =1 , n ) c a l l m p i _ f i n a l i z e ( i e r r ) s t o p
end p r o g r a m M P I b r o a d c a s t
call mpi bcast(x,n,MPI REAL8,iroot,MPI COMM WORLD,ierr)
ある変数について，rank が iroot のプロセスでの値を全プロセスに送信し，全プロ セスで rank が iroot のプロセスでの値と同一の値をとるようにする． • 1 つ目の引数 x: 送信元では送信するデータの入った配列．受信側ではデータ を受信するための配列． • 2 つ目の引数 n: 送信するデータ配列の長さ． • 4 つ目の引数 z: 送信元の rank の番号． 上のプログラムでをプロセス数 4 で実行した場合，出力は以下のようになる．
B e f o r e ... R a n k : 0 x : 1.0 2.0 3.0 4.0 5.0 B e f o r e ... R a n k : 1 x : -1.0 -1.0 -1.0 -1.0 -1.0 B e f o r e ... R a n k : 2 x : -1.0 -1.0 -1.0 -1.0 -1.0 B e f o r e ... R a n k : 3 x : -1.0 -1.0 -1.0 -1.0 -1.0 A f t e r ... R a n k : 0 x : 1.0 2.0 3.0 4.0 5.0 A f t e r ... R a n k : 1 x : 1.0 2.0 3.0 4.0 5.0 A f t e r ... R a n k : 2 x : 1.0 2.0 3.0 4.0 5.0 A f t e r ... R a n k : 3 x : 1.0 2.0 3.0 4.0 5.0
4.4
一対一での通信
ある特定のプロセスから別のある特定のプロセスへデータを送る場合には，mpi send と mpi recv という subroutine を，例えば以下のようにして用いる．プログラム 6: MPI による一対一での送受信 p r o g r a m M P I c o m m i m p l i c i t n o n e i n c l u d e ’ m p i f . h ’ i n t e g e r :: m p i s t a t ( M P I _ S T A T U S _ S I Z E ) i n t e g e r :: nprocs , irank , i r o o t i n t e g e r :: ierr , istat , i t a g i n t e g e r :: isrc , i d e s t i n t e g e r :: nn , i integer , p a r a m e t e r :: n =2 r e a l (8) :: x ( n ) i s r c =0 i d e s t =1 x ( : ) = 0 . 0 c a l l m p i _ i n i t ( i e r r ) c a l l m p i _ c o m m _ s i z e ( M P I _ C O M M _ W O R L D , nprocs , i e r r ) c a l l m p i _ c o m m _ r a n k ( M P I _ C O M M _ W O R L D , irank , i e r r ) if ( i r a n k == i s r c ) t h e n do i =1 , n x ( i )= d b l e ( i ) end do end if if ( i r a n k == i d e s t ) t h e n w r i t e (6 , ’ ( a10 ,2 f6 .1) ’ ) ’ B e f o r e ... ␣ ’ ,( x ( i ) , i =1 , n ) end if
i t a g =0
if ( i r a n k == i s r c ) t h e n
ca l l m p i _ s e n d ( x , n , M P I _ R E A L 8 , idest , itag , M P I _ C O M M _ W O R L D ,& i e r r )
end if
if ( i r a n k == i d e s t ) t h e n
ca l l m p i _ r e c v ( x , n , M P I _ R E A L 8 , isrc , itag , M P I _ C O M M _ W O R L D ,& mpistat , i e r r ) end if if ( i r a n k == i d e s t ) t h e n w r i t e (6 , ’ ( a10 ,2 f6 .1) ’ ) ’ A f t e r ... ␣ ␣ ’ ,( x ( i ) , i =1 , n ) end if c a l l m p i _ f i n a l i z e ( i e r r ) s t o p end p r o g r a m M P I c o m m
call mpi send(x,n,MPI REAL8,idest,itag,MPI COMM WORLD,ierr) データを rank が idest のプロセスへ送信する．
• 1 つ目の引数 x: 送信するデータの入った配列． • 2 つ目の引数 n: 送信するデータ配列の長さ． • 4 つ目の引数 idest: 送信先の rank．
• 5 つ目の引数 itag: 送信データに付けるタグ (識別番号)．
call mpi recv(x,n,MPI REAL8,isrc,itag,MPI COMM WORLD,mpistat,ierr) Rankが isrc のプロセスから送信されたデータを受信する．
• 1 つ目の引数 x: 受信したデータを格納する配列． • 2 つ目の引数 n: 受信するデータ配列の長さ． • 4 つ目の引数 idest: 送信元の rank．
• 5 つ目の引数 itag: 送信データに付いた識別番号で，mpi send で指定したの と同じ値を入れる必要がある．
• 7つ目の引数 mpistat: 長さがMPI STATUS SIZEの整数配列を与える．MPI STATUS SIZE は ‘mpif.h’ で定義されている．
なお，mpi send，mpi recv はデータの送信，受信手続きを終えるまではその後の処理 を実行せずに待ち状態になる．プログラムの書き方がまずいと，全プロセスが送信終了 待ちあるいは受信終了待ちになってプログラムが固まってしまう場合があるので，送受 信の順序には注意する必要がある．
送信，受信手続きを完了したかどうかに関わらず次の処理に進む mpi isend，mpi irecv という subroutine もあるが，これを使った場合，送信手続きが終わる前に送信すべきデータ が消されてしまったり，受信手続きが終わる前に受信データを読みにいってしまったりとい った問題が起こるおそれがある．このようなことを防ぐために，mpi isend や mpi irecv を 使う場合は，送受信に関わるメモリ領域にアクセスがある前で mpi wait という subroutine を call しておき，それ以上は送受信が完了するまで処理が進まないようにしておく必要 がある．
4.5
構造体
各粒子の状態変数を構造体 (derived type) にまとめた場合，derived type の各要素の値 をまとめて別のプロセスに送信できた方が便利である．mpi type struct は，mpi send， mpi recvを 1 回 (ずつ)call すればデータをまとめて送れるようにするために，MPI 上の 構造体を定義する subroutine である．
ただし，Fortran 構造体の各要素がどのようにメモリ上に配置されるのかは，処理系に 依存するので，「MPI 構造体」を定義する方法もそれによって変わってくる可能性がある． 以下で示す例は，Intel Fortran compiler 10.1 の環境では動いているが，処理系によって は MPI 構造体の定義の仕方を変える必要が生じる可能性があるので，留意されたい．
筆者がよく扱っているモデルでは，状態変数を以下のような Fortran 構造体にまとめて いる．
プログラム 7: Fortran 構造体 (derived type) の定義例
m o d u l e r a m p a r a m s use g r i d _ p a r a m e t e r s p u b l i c type , p u b l i c :: d r i f t _ p a r a m s s e q u e n c e r e a l :: ddp
r e a l :: p o t e n t ( nr , np ) r e a l :: vcl ( nr , np ) , p1 ( nr , np ) r e a l :: p o t c o e f ( 0 : 2 4 , 0 : m a x t i m e s e q ) end t y p e d r i f t _ p a r a m s type , p u b l i c :: b o u n d _ p a r a m s s e q u e n c e r e a l :: fb ( np , nw , nk ) r e a l :: d e n s _ p a r (8) , t e m p e _ p a r (8) i n t e g e r :: ib0 ( np ) end t y p e b o u n d _ p a r a m s type , p u b l i c :: r a m _ d a t a s e q u e n c e r e a l :: f2 ( nr , np , nw , nk ) r e a l :: t t y p e ( d r i f t _ p a r a m s ) :: d _ p a r s r e a l :: h d e n s _ f a c t (0: m a x t i m e s e q ) t y p e ( b o u n d _ p a r a m s ) :: b _ p a r s i n t e g e r :: nid end t y p e r a m _ d a t a end m o d u l e r a m p a r a m s sequenceは，Fortran 構造体の各要素をメモリ上に書いた順番に並べることを指示する ものである．
上のような Fortran 構造体を MPI で送信するための MPI 構造体は，Intel compiler の 環境だと，以下のような subroutine で多分生成できるのではないかと思われる． プログラム 8: MPI 構造体の定義例 s u b r o u t i n e g e n _ t y p e ( m o d e l _ t y p e ) integer , p a r a m e t e r :: n b l o c k =13 , n c o m b l o c k =12 ! In o r d e r to p r e s e r v e the id for e a c h i n s t a n c e , ! the l a s t b l o c k is not c o p i e d in r e s a m p l i n g . integer , d i m e n s i o n ( n b l o c k ) :: b l o c k _ l e n , b l o c k _ d i s p l , t y p e s i n t e g e r :: m o d e l _ t y p e integer , p a r a m e t e r :: b y t e 1 r d =8 , b y t e 1 i d =4 ! B l o c k l e n g t h for e a c h e l e m e n t b l o c k _ l e n ( 1 ) = nr * np * iw * ik b l o c k _ l e n ( 2 ) = 1 b l o c k _ l e n ( 3 ) = 1 b l o c k _ l e n ( 4 ) = nr * np b l o c k _ l e n ( 5 ) = nr * np b l o c k _ l e n ( 6 ) = nr * np b l o c k _ l e n ( 7 ) = 2 5 * ( m a x t i m e s e q +1) b l o c k _ l e n ( 8 ) = m a x t i m e s e q +1
b l o c k _ l e n ( 9 ) = np * iw * ik b l o c k _ l e n ( 1 0 ) = 8 b l o c k _ l e n ( 1 1 ) = 8 b l o c k _ l e n ( 1 2 ) = np b l o c k _ l e n ( 1 3 ) = 1 ! T y p e for e a c h e l e m e n t do i =1 ,11 t y p e s ( i )= M P I _ R E A L 8 end do do i =12 ,13 t y p e s ( i )= M P I _ I N T E G E R end do ! B l o c k d i s p l a c e m e n t b l o c k _ d i s p l ( 1 ) = 0 * b y t e 1 r d do i =2 , n b l o c k s e l e c t c a s e ( t y p e s ( i - 1 ) ) c a s e ( M P I _ R E A L 8 ) b l o c k _ d i s p l ( i )= b l o c k _ d i s p l ( i - 1 ) + b y t e 1 r d * b l o c k _ l e n ( i -1) c a s e ( M P I _ I N T E G E R ) b l o c k _ d i s p l ( i )= b l o c k _ d i s p l ( i - 1 ) + b y t e 1 i d * b l o c k _ l e n ( i -1) c a s e ( M P I _ C H A R A C T E R ) b l o c k _ d i s p l ( i )= b l o c k _ d i s p l ( i - 1 ) + b l o c k _ l e n ( i -1) end s e l e c t end do ca l l m p i _ t y p e _ s t r u c t ( n c o m b l o c k , b l o c k _ l e n , b l o c k _ d i s p l , types ,& m o d e l _ t y p e , i e r r ) ca l l m p i _ t y p e _ c o m m i t ( m o d e l _ t y p e , i e r r ) ! c a l l m p i _ t y p e _ e x t e n t ( m o d e l _ t y p e , iext , i e r r ) r e t u r n end s u b r o u t i n e g e n _ t y p e
call mpi type struct 新しい型を作る．
• 1 つ目の引数 ncomblock: 新しい型に含まれる要素の個数．
• 2 つ目の引数 block len: 長さ ncomblock の配列．各要素の長さを指定． • 3 つ目の引数 block displ: 長さ ncomblock の配列．各要素が 1 つめの要素か
• 4つ目の引数 types: 長さncomblockの配列．各要素の型をMPI REAL8, MPI INTEGER などの値で指定．
• 5 つ目の引数 model type: 新たに作られる型を参照するための変数． call mpi type commit
mpi type structで新しい型 model type を作った後，これを call することで，新 しい型が MPI で使えるようになる．
このうち，3 つ目の引数 block displ は，処理系によって与え方を変える必要があるか もしれない．Intel compiler では，derived type を定義するときに sequence を付けておく と，各要素が隙間なく sequential に並べられる仕様になっているのでこれで問題はないが， 処理系によっては隙間が入ったりするので，そのときは block displ を適宜調整する．
構造体がうまく定義できていたら，以下のようにして mpi send，mpi recv を用いて， データのやりとりができるようになるはずである．
プログラム 9: MPI による一対一での送受信
t y p e ( r a m _ d a t a ) :: ram
! ! 中略()
if ( i r a n k == i s r c ) t h e n
ca l l m p i _ s e n d ( ram ,1 , m o d e l _ t y p e , idest , itag , M P I _ C O M M _ W O R L D , i e r r ) end if
if ( i r a n k == i d e s t ) t h e n
ca l l m p i _ r e c v ( ram ,1 , m o d e l _ t y p e , isrc , itag , M P I _ C O M M _ W O R L D ,& mpistat , i e r r ) end if ただし上述のように，構造体の定義の仕方は処理系に依存する部分もあるので，構造体 を定義したときは，それがちゃんとうまく送受信できるのかきちんと確認をしておくべ きである． わざわざ MPI 構造体を定義するのが面倒であれば，一旦別の配列に押し込んでおいて その配列を送信する，あるいは各要素を逐一送信する，といった方法を取ってもよいだろ う．ただ，MPI の通信 subroutine は，subroutine を呼ぶだけでも多少時間がかかってし まうので，同じデータ量をやりとりする場合でも，データを小分けにして mpi send を何 度も call するよりは，ある程度データをまとめて送った方が，計算時間が少なく済むこと
が多い．したがって，少なくとも，各要素を逐一送信するという手段を取るよりは何らか の形でデータをまとめて送ることを考えた方がよいと思われる．
5
並列計算による粒子フィルタ
粒子フィルタでは，システムノイズを与えてシミュレーションを進めるという予測のス テップについては，完全に並列処理が可能である．一方，フィルタの部分では，各プロセ スで得られる粒子 x(i)t|t−1 の情報を交換する必要が生じるので，MPI のプロセス間通信の 出番となる． フィルタリングの手続きは，上述のように 1. 各粒子の尤度を計算する． 2. 尤度に基づいて，各粒子の重みを計算する． 3. 各粒子の重みにしたがって，どの粒子の複製をどの粒子と置き換えるかについての 対応テーブルを作る． 4. 対応テーブルに沿って，実際に粒子の分配操作を行う． の 4 つの段階に分けられるが，順を追ってどのように MPI を使って実装するかを紹介す ることにする．5.1
尤度の計算
尤度の計算には観測データが必要となるが，データはファイルから入力することが多 く，MPI を使わずに各プロセスが直接データを読み込むようにするのが普通であると思 われる．この場合，尤度の計算に必要な観測データと x(i)t|t−1の情報は，どちらも個々のプ ロセスが自分で取得できることになるので，MPI を使うことなく，個々のプロセスで尤 度を計算させればよい．5.2
重みの計算
重みを計算するには，尤度 λ(i)t の全粒子についての総和∑Ni=1λ(i)t が必要となるので，一 旦，一つのプロセスに尤度の値を集める必要がある．そこで，mpi gather を使う1)． 全粒子数を nptotal 個とし，その nptotal 個の粒子を 1 プロセスあたりに m 個ずつ割 り振っていたものとする．このような場合に，rank が iroot のプロセスに全粒子の尤度 の情報を集めるには，以下のようにするとよいだろう． プログラム 10: 尤度の値の回収 i n t e g e r :: p r e c _ r e a l = M P I _ R E A L 8 r e a l (8) :: x l a m b d a ( m ) r e a l (8) :: x l a m a r r ( n p t o t a l ) ! ! 中略()
c a l l m p i _ g a t h e r ( xlambda , m , prec , xlamarr , m , prec , iroot ,& M P I _ C O M M _ W O R L D , i e r r ) なお，xlambda はそれぞれのプロセスで割り当てられた粒子についての尤度を保持した 配列である．
5.3
リサンプリングの対応テーブル
対応テーブルの作成は，重みの情報を特定のプロセスに集めてしまえば，プログラム 2 のやり方と同様の方法でできる．その次のステップのために，対応テーブルを全プロセ スに送信する必要があるが，それは mpi bcast で実現すればよい．5.4
粒子の再分配
前のステップで得たテーブルにしたがって破棄する粒子の場所に増やすべき粒子の複 製を配置するという手続きを行っていく．並列計算の場合，新たに配置する粒子 (の複製) が同じプロセスのものとは限らない—つまりプロセス内で新たに配置すべき粒子の情報1)単に総和を取るだけなら，mpi reduce という subroutine を使う方法もあるのだが，ここでは一つのプ
を保持しているとは限らない—ので，そのような場合には，mpi send と mpi recv を使っ て粒子の情報を他プロセスから運んでくる必要がある．粒子フィルタでは，この粒子の 情報の受け渡しが様々なプロセスでばらばらに発生するので，あるプロセスでは送信用 の routine mpi send が call され，また別のプロセスでは受信用の routine mpi recv が call されるといったような形で，ややプログラムは複雑な構造になってしまう．
以下のサンプルでは，各プロセスがそれぞれリサンプリングテーブル isrc を逐一調べ ていく．送信 rank と受信 rank が同じでかつそれが自分の rank であったら粒子のコピー を行う．送信 rank と受信 rank が異なった場合には，自分が送信 rank なら mpi send を call，自分が受信 rank なら mpi recv を call するという流れで粒子の再分配を行う．一 見，複雑な構造はしているが，どのプロセスも同じリサンプリングテーブルを同じ順序 で参照するので，一応，データの入れ違いなどはたぶん起こらないはずである． プログラム 11: 粒子の再分配手続き c a l l m p i _ c o m m _ s i z e ( M P I _ C O M M _ W O R L D , nprocs , i e r r ) c a l l m p i _ c o m m _ r a n k ( M P I _ C O M M _ W O R L D , irank , i e r r ) i t a g =0 do k =0 , nprocs -1 ! 受信側のプロセスの ra n k do m =1 , n p e a c h ! 受信プロセス上での粒子の席番号 i d e s t = m + k * n p e a c h ! 粒子の席の全体での通し番号 ! k s r c : i d e s t の席に座る粒子の全体での通し番号 ! m s r c : i d e s t の席に座る粒子が存在するプロセスの ra n k k s r c =( i s r c ( i d e s t ) -1)/ n p e a c h m s r c = i s r c ( i d e s t ) - k s r c * n p e a c h t r a n s =. t r u e . ! 送信側と受信側とで r a n k が等しければ if ( k == k s r c ) t h e n t r a n s =. f a l s e . ! R a n k が i r a n k ならコピー if ( k == i r a n k ) t h e n xp ( m )= xp ( m s r c ) end if e l s e ! m > 1 ならすでに同じ粒子を受けていないか調べる if ( m > 1) t h e n do mm =1 , m -1 if ( i s r c ( i d e s t )== i s r c ( mm + k * n p e a c h )) t h e n t r a n s =. f a l s e . if ( i r a n k == k ) xp ( m )= xp ( mm )
ex i t end if end do end if ! プロセス間で粒子を移動する必要が生じた場合 if ( t r a n s ) t h e n ! R a n k が k s r c なら送信 if ( i r a n k == k s r c ) t h e n c a l l m p i _ s e n d ( xp ( m s r c ) ,1 , m o d e l _ t y p e , k , itag ,& M P I _ C O M M _ W O R L D , i e r r ) ! R a n k が k なら受信 el s e if ( i r a n k . eq . k ) t h e n c a l l m p i _ r e c v ( xp ( m ) ,1 , m o d e l _ t y p e , ksrc , itag ,& M P I _ C O M M _ W O R L D , mpistat , i e r r ) end if ! 関係ない r a n k なら何もしない it a g = i t a g +1 if ( i t a g > M P I _ T A G _ U B ) i t a g =0 end if end if end do end do
