中 慎也
中野
情報・システム研究機構 統計数理研究所

A li ti off thti l filt
Application
the particlefilter
 The particle filter tends to be computationallyexpensive in applying to high-dimensional systems.
 An enormous number of particles is required in order

to appropriately approximate a PDF.
 Computation
pfor temporalpevolution is required
qfor

each of the particles.
 Computation
pfor temporalpevolution for each p
particle isusually expensive but it can be performed in parallel.

C t tContents
 粒子フィルタの概要
 超並列計算機への実装に向けての問題点
 超並列計算機への実装に向けたアルゴリズム
 仮想データによる実験
仮想デ タ実験

P ti l filt
ParticlefilterIn the particle filter, a probability
density function (PDF) isapproximated by a large number
pdrawn from the PDFof samples
as1 Np ( x)    ( x  x (i ) ).
N i 1Each of the samples is calleda ‘particle’.

1 Np( x)    ( x  x (i ) ) N i 1

1 N(i )
(x)
p(x)
dx


(x)

N i 1

P ti l filt
Particlefilter(e.g., Gordon et al. 1993) If we define an operator describing the time transition of

the system state as f, a posterior particle at the previous
step evolves into a predictive particle at the current step
(i )
(i )
(i )
as xk |k 1  f ( xk 1|k 1 )  vk for each i,where1 N
p ( xk 1 | yk 1 )    ( xk 1  xk(i)1|k 1 ).
N i 1

P ti l filt
Particlefilter
An approximation of theposterior distribution is
obtained by resampling theparticles
ti l with
ith replacementltsuch that each particle aredrawn with the probability
proportional to the likelihood.

p ( xk | y1:k )  p( yk | xk ) p ( xk | y1:k 1 )1 Np ( xk | y1:k 1 )    ( xk  xk(i|k) 1 )N i 1

The number of duplicates  The likelihood for each prior particle

D i ti off th
ti l filt
Derivation
the particle
filter1 N(i )
 If we approximate the prior PDF as p ( xk | y1:k 1 )    ( xk  xk |k 1 )N i 1the posterior PDF can be approximated asp ( xk | y1:k ) p ( yk | xk ) p ( xk | y1:k 1 )

 p( y | x ) p( x | yk

 Nk1:k 1

k1

N(i )
(i )
p(y|x)
(xx k k|k 1 k k|k 1 )

(i )
i 1
p
(y|x)
 k k|k 1
j 1N  wk(i ) ( xk  xk(i|k) 1 ).
i 1

(i )
where wk p( yk | xk(i|k) 1 )Ni 1) dxk(i )k |k 1p ( yk | x).

D i ti off thti l filt
Derivationthe particlefilter A new population {xk |k , , xk |k } is obtained by resampling
(1)

(N )

the predictive population{xk(1)|k 1 , , xk( |Nk ) 1} with replacement
such that each particle from the predictive population are
drawn with probability proportional to the weight wk(i.)(i ) The number of duplicates mk satisfies

mk(i )  Nwk(i )

(i )
m

0;
m
 kk 0(i )for each k. Then we can approximate the posterior PDF asmk(i )p( xk | y1:k )   ( xk  xk(i|k) 1 )i 1 N1 N
   ( xk  xk(i|k) ).
N i 1N

高次元問題 の粒子フ ルタの適用
高次元問題への粒子フィルタの適用
 できれば，数十，数百次元以上の非線型問題を解きたい．
できれば 数十 数百次元以上の非線型問題を解きたい 一方，所要粒子(=サンプル)数は，次元の増大に対して指数関数的に増大する．
関数的に増大する高次元問題への適用(=粒子数の増大)に対応するための高次元問題適用( 粒子数 増大) 対応するため手段の一つとして，超並列計算機の利用が挙げられる．
昨今のスカラー超並列計算機の発展は，比較的規模の大きな問題に対しても粒子フィルタを適用できる可能性を切り開く可能性がある．

分散メモリ型並列計算機における粒子フィルタ

予測: 容易に並列化可能

フィルタ: ランダムにノード間の通信が発生分散メモリ型並列計算機への粒子フィルタの実装
 一期先予測期先予測 PE (Processing Element)数に比例する形で粒子数を増やしても計算時間には影響しないと思われる．

 フィルタ(=リサンプリング) 大量の粒子のどれをどこに移動させるか，エラーを起こさずにスケジューリングするのはかなり難しい．
 任意の2 nodes間でほぼat randomに通信が発生し，並列化しにくい．
 PE数(ノード数)が増えるほど，大量の通信が発生．各ノードにとっての通信先が増えることになるので，計算ノドにと ての通信先が増えることになるので 計算効率が著しく低下しうる．

GGroupingi
We can divide the whole population into multiple
sub-pupulations (groups).
(k ) : the weight for each groupk(  , ): the local weight for each particle in a group

The global weight in the whole population is given by the productof the two weights as  k k( )(  , )for each particle.

D i ti off thti l filt
Derivation
the particle
filter1 N(i )
 If we approximate the prior PDF as p ( xk | y1:k 1 )    ( xk  xk |k 1 )N i 1the posterior PDF can be approximated asp ( xk | y1:k ) p ( yk | xk ) p ( xk | y1:k 1 )

 p( y | x ) p( x | yk

 Nk1:k 1

k1

N(i )
(i )
p(y|x)
(xx k k|k 1 k k|k 1 )(i )i 1p(y|x)
 k k|k 1
j 1N  wk(i ) ( xk  xk(i|k) 1 ).
i 1(i )where wk p( yk | xk(i|k) 1 )

Ni 1) dxk(i )k |k 1p ( yk | x).

D i ti off th
ti l filt
Derivation
the particlefilter A new population {xk |k , , xk |k } is obtained by resampling
(1)

(N )

the predictive population{xk(1)|k 1 , , xk( |Nk ) 1} with replacement
such that each particle from the predictive population are
drawn with probability proportional to the weight wk(i.)(i ) The number of duplicates mk satisfies

mk(i )  Nwk(i )(i )m0;m kk 0(i )for each k. Then we can approximate the posterior PDF as
mk(i )p( xk | y1:k )  
 ( xk  xk(i|k) 1 )i 1 N1 N   ( xk  xk(i|k) ).
N i 1NGiGrouping If the constraints

(k )k(  , )  wk([  1]  ) ,(  , ) k  1,

 1

( )
 k  1.

 1

are satisfied, we can approximate the posterior PDF as:
Np( xk | y1:k )   w  ( xk  xi 1

(i )
k(i )
k |k 1)  1( )k(  , )([  1]  )(xx).
 kkk |k 1

 1

 We can then consider the local weights for each particle k(  , )and the weights for each group (k ) separatelyseparately.

GGroupingi リサンプリングをグループごとに行っても，事後分布の近似リサンプリングをグル プごとに行っても 事後分布の近似となる．
 グループ内で退化(=アンサンブルの多様性が失われる)が起こりやすくなるが，手続きは並列化しやすくなる．
p ( xk | y1:k )    1  1( )
k( )
k



(  , )([  1]  )


(xx)
 kkk |k 1

 1



1



(k )  1

1



mk(  , )

 ( xk  xk([|k11]  ) ) (x  x    ) k([ 1]  )k |kGGroupingi ここでは，2次元トここでは，2次元トーラスモデルを考ラスモデルを考え，各ノードに均等に粒子を割り当てるものとする．
 さらに，2x2個のノードに割り当てられた粒子を一つのグループに属するものと見なす．

Gi strategy
t tGrouping グル
グループ分けを毎ステップ固定する
プ分けを毎ステップ固定する必要はない. (Bolić et al. 2005)
 そこで，毎ステップグループ分けを切り替えることを考える．

Alt
ti llattice
tti patterntt
it hi
Alternativeswitching 2つのグループ分けを交互に用いる．
 グループ分けの切り替えを何ステップか繰り返せば，あるノードのグプ替
何プ 繰
ば情報が全体に伝搬する．
 2次元トーラスモデルの各ノードを実際のPEに対応づけることで，リ次元ト ラデ の各
ド 実際の に対応 ける と ，リサンプリングを各グループで並列に実行することが可能．

EExperimentsit Lorenz 96 model (Lorenz and Emanuel 1998)
1998).

dx jddt

 ( x j 1  x j  2 ) x j 1  x j  8for j  1,,, 40( x1  x39 , x0  x40 , x41  x1.)
 A simple 1-dimensional nonlinear model. (40 variables) 1 time step Δt is set to be 0.005.
 We assume that half of the 40 state variables can beobserved.

EExperimentsit An artificial data set is generated

every 20 time step by adding
errors with a standard deviation of1.5 to a ‘true’ run.
 The system noise was assumed

to be a Gaussian with a diagonalcovariance as diag(0.25, 0.25, 0.25).

EExperimentsit The following
follo ing obser
observation
ation

model p(yk|xk) is assumed: || yk  Hxk |k 1 ||2 
1p ( yk | xk |k 1 ) exp  
.
22
2


 The observation noise wasassumed to be a Gaussian witha diagonal covariance asdiag(σ2, σ2, σ2) where σ= 3.0.

Rlt
ResultPEs

Particles

Global PF

Alternative switching

RMSE Elapsed times

RMSE Elapsed times

4x465536

0 920.9201m34s

0 94
0.9400m54s

8x8262144

0.8504m56s

0.8600m59s

16x1610485760.81

18m41s

0.8101m02s

32x32

4194304

0.7973m54s

0.8001m04s

 TheTh experiments
it were performed
fd on CRAY XT6XT6m (CPU(CPU: OOpteront
6174
12 core x2, Memory: 32GB; for each node).
 A virtual two-dimensional torus network is assumed. Each core in thereal architecture is assigned to each node in the virtual torus network.

Rlt
ResultPEs

Particles

RMSE

Elapsed times Effective
sample size

4x4

2621440.86

06m14s75660

8x8262144

0.8600m59s

72403

16x16262144

0.8900m12s

55381

32x32

2621440 95
0.9500m04s

18564

 粒子数を固定してノード数を増やすと，制度が落ちていく．
 有効粒子数 (ここでは，エントロピー S  義) も落ちていく．

N w log w に基づいて定i 1ii

Eff ti numberb off particles
ti l
Effective
PEs

Particles RMSE Elapsed
times

Effective
Effective
Effective
sample size number of sample size
nodes
for a node

4x4

2621440.86

06m14s

75660

15.8

4803.08x82621440.8600m59s

72403

58.5

1237.716x16

2621440.8900m12s55381

166.2

333.6

32x32

2621440.9500m04s

18564

201.3

92.2





S   (k )k(  , ) log (k )k(  , )
 1  1   1   1( )
k( )
k



(  , )( )(  , )
log


log



 kkk

 1log ( )
k



 
 1( )k(  , )(  , )
log

 kk

 1

E llik lih d dit ib tiEcology
off llog-likelihood
distribution

8x832x32

 粒子数を262144個に固定し，グループあたりのノード数も2x2個に固定した場合．
 ノード数を増やすと，情報が全体に伝わるのに必要なステップ数が増えて，ノードごとの尤度のばらつきが大きくなり，それが性能に影響する．

Rlt
ResultPEsParticles

RMSE

Elapsed times Effective
sample size

4x465536

0.9400m54s

18751

8x8

2621440.8600m59s

72403

16x16

10485760.8101m02s

25474032x324194304

0 80
0.8001m04s

839029

 但し，少なくとも今回実験した条件の範囲においては，ノードあたりの但し 少なくとも今回実験した条件の範囲においては ノ ドあたりの粒子数を固定すれば有効粒子数も増大するという結果が得られた．したがって，ノード数を増やすこと自体には意味があると考えられる．

Eff ti numberb off particles
ti lEffective
PEsParticles RMSE Elapsed
times

Effective
Effective
Effective
sample size number of sample size
nodes
for a node

4x4

65536

0.9400m54s

18751

15.5

1207.18x82621440.8600m59s

72403

58.5

1237.716x161048576

0.81

01m02s254740201.7

1261.432x324194304

0.80

01m04s

839029667.8

1256.4S   (k )k(  , ) log (k )k(  , )
 1  1
  
 1
   1( )k( )k(  , )( )(  , )
log
log

 kkk

 1log ( )k  1( )k(  , )(  , )
log

 kk

 1

Summary 超並列計算機へ実装するに際しては，リサンプリング手続きに伴う大量のノード間通信，およびリサンプリングの並列化の困難さが懸念される．
 粒子のグループ分けは，単純なだけに並列化・高速化の手段としてある程度有効と見られる．
る程度有効と見られる
 グループあたりのノード数を固定した場合，ノード数を増やすに連れて退化が起 りやすくなる退化が起こりやすくなる．
 しかし，各ノードに割り当てる粒子の数を一定とすれば，ノード数が増えるに従って有効粒子数も増やせるので，速く解くという方向ではなく，大規模な問題を解くという方向であれば，ノード数を増やす(=並列度を上げる)ことには意味があると考えられる．