Merging particle ﬁlter for high-dimensional nonlinear problems
S. Nakano, G. Ueno, and T. Higuchi(The Institute of Statistical Mathematics, ROIS, Tokyo, JAPAN)

1Introduction

Sequential data assimilation is an on line approach in which a state at each time
step is estimated as the observation is obtained. Most of sequential data assimilation techniques are based on Bayesian approach. The information of observations
is incorporated into the system model by considering a posterior probability density function (PDF) derived from a prior PDF, where the prior PDF is obtainedusing past data and a system model. The particle ﬁlter (PF) (Gordon et al., 1993;
Kitagawa, 1993, 1996) is one of such sequential data assimilation algorithms. The
PF approximates PDFs at each time step by an ensemble of a large number ofparticles. An estimation of a posterior PDF is obtained by resampling with replacement from a prior ensemble. Since the PF does not require assumptions oflinearity or Gaussianness, it is applicable to general nonlinear problems including cases with nonlinear observations which other algorithms such as the ensemble
Kalman ﬁlter (Evensen, 1994; Burgers et al., 1998) do not provide good estimation.
However, the PF often encounters a problem called ‘degeneration’, which does
not occur in the EnKF. As resampling procedures are applied recursively, most of
the particles are replaced by particles that ﬁt the observed data better, and theposterior PDF is eventually represented by only a few of the particles among themembers of the initial ensemble. This reduces the validity of ensemble approximation especially in applying to high-dimensional systems. Although this problem
could be avoided by increasing the number of particles in the ensemble, the increase of the number of particles requires a prohibitive computational cost at each
forecast step.
In order to overcome this problem, we have devised another technique, themerging particle ﬁlter (MPF). In the MPF, a ﬁltering procedure is performed bymerging several particles of a prior ensemble, which is rather similar to the genetic
algorithm (e.g., Goldberg, 1989). This merging procedure is performed such thatthe ﬁrst and second moments of a posterior PDF is preserved. The MPF allows
the degeneration problem to be avoided without unduly increasing the number ofparticles and thus the MPF is applicable to higher-dimensional systems. Moreover,1the primary advantage of the PF over the EnKF is inherited; that is, the MPF isapplicable even to cases in which the relationship between a state and observed
data is nonlinear. In our previous paper (Nakano et al., 2007), the basic idea ofthe MPF is given. The purpose of this paper is to brieﬂy review the algorithm of
the MPF and to discuss some properties of the MPF.

2

Particle ﬁlter
In the following, we consider a state space model as:
xk = Fk (xk−1 , v k )(1a)y k = Hk (xk ) + w k(1b)

where the vectors xk and y k indicate the state of a system and observed data at a
discrete time T = tk (k = 1, . . .), respectively, and the vectors v k and w k denotesystem noise and observation noise, respectively. The operator Fk represents the
temporal evolution from time tk−1 to time tk according to the simulation. The
operator Hk projects the state vector xk to the observation space.
The PF considers a PDF of a state xk . The PDF is approximated by an ensemble
consisting of a large number of discrete samples called ‘particles’. Suppose that aﬁltered distribution at time T = tk−1 , p(xk−1 |y 1:k−1 ), is approximated by particles
(1)
(2)
(N )
{xk−1|k−1, xk−1|k−1, · · · , xk−1|k−1} as

1  
(i)
δ xk−1 − xk−1|k−1p(xk−1 |y 1:k−1 ) ≈
N i=1
N

(2)

where δ is Dirac’s delta function, and N is the number of particles in the ensemble.
Here we expressed p(xk−1 |y 1 , · · · , y k−1 ) as p(xk−1 |y 1:k−1 ). We obtain an ensemble
approximation of the forecast distribution of the state at the next observation time
T = tk as
N

1  
(i)
p(xk |y 1:k−1 ) ≈δ xk − xk|k−1
(3)
N i=1
(i)

(i)

(i)

where xk|k−1 is given by Fk (xk−1|k−1 , vk ) for each particle i. This procedure iscalled the forecast step.
From the forecast distribution p(xk |y 1:k−1) and observed data y k , we obtain an2

approximation of the ﬁltered PDF p(xk |y 1:k ):p(xk |y 1:k−1 ) p(y k |xk )p(xk |y 1:k−1 ) p(y k |xk )dxk
N
 



1
(i)
(i)

≈ 
p y k |xk|k−1 δ xk − xk|k−1(j)i=1
j p y k |xk|k−1

p(xk |y 1:k ) = 

=

N

i=1

(4)



(i)
wi δ xk − xk|k−1

(i)

(i)

where p(y k |xk|k−1) is the likelihood of xk|k−1 given the data y k and the weight wi
is deﬁned as(i)p(y k |xk|k−1).
(5)
wi = (j)j p(y k |xk|k−1 )This is called the ﬁltering step.
(1)(N )We then obtain a new ensemble {xk|k , · · · , xk|k } by resampling the forecast
(1)

(N )

ensemble {xk|k−1, · · · , xk|k−1} with replacement. The new ensemble may containmultiple duplicates originating from the same particle in the forecast ensemble.
(i)
The resampling is performed such that the number of copies for a particle xk|k−1
is nearly proportional to the weight wi . Denoting the number of copies for a(i)particle xk|k−1 as mi , mi satisﬁes(6)
mi = N; mi ≥ 0
mi ≈ Nwi(i)for each xk|k−1. Using Eq. (6), Eq. (4) can be further approximated asp(xk |y 1:k ) ≈

N


mi 
(i)
δ xk − xk|k−1
N
i=1


1  
(i)
=
δ xk − xk|k .
N i=1
N

(7)

Thus, the newly generated ensemble approximates the ﬁltered PDF p(xk |y 1:k ).
Equation (7) has the same form as Eq. (2). We can then recursively repeat theabove procedure from Eq. (2) to Eq. (7) to incorporate a sequence of observed
data into the system model.

3

Merging particle ﬁlter

The MPF is a variant of the PF. In the MPF, each particle of a ﬁltered ensemble
is generated by combining multiple particles drawn from the forecast ensemble,

3

which is rather similar to the genetic algorithm. When the number of particles tobe combined is denoted by n, it is necessary to draw n×N samples from the forecast
(1,1)
(n,1)
(1,N )(n,N )
ensemble. to obtain an ensemble: {x̂k|k , · · · , x̂k|k , · · · , x̂k|k , · · · , x̂k|k }. The(i)sampling of the n × N samples is performed such that a particle xk|k−1 is drawn
with a probability wi like the PF described in the previous section. It should be(j,1)(j,N )
noted that a subset {x̂k|k , · · · , x̂k|k } aﬀords an ensemble approximation of theﬁltered PDF as
N

1  
(j,i)
δ xk − x̂k|kp(xk |y 1:k ) ≈
(8)
N i=1
(i)

because copies of each particle xk|k−1 is contained with a ratio proportional to wi in
(1)

(N )

the subset. We then make a new ensemble consisting of N particles {xk|k , · · · , xk|k }.
Each particle in this new ensemble is generated as a weighted sum of n samples
from the n × N sample set as:(i)
xk|k =

n

j=1

(j,i)

αj x̂k|k .

(9)In order that the newly generated ensemble approximately preserves the mean andcovariances of the ﬁltered PDF, the weights αj are set to satisfynαj = 1(10a)

αj2 = 1

(10b)

j=1
n

j=1

where each αj is a real number. When the merging weights satisfy Eq. (10a), the(1)(N )mean of the PDF approximated by the new ensemble {xk|k , · · · , xk|k } becomes


N
N
N
n

1  
1  (i)
1 
(i)
(j,i)
δ xk − xk|k dxk =x =αj x̂k|kxk
N i=1
N i=1 k|k N i=1 j=1
 
n
N


1  
(j,i)
=
δ xk − x̂k|k dxk
αj xk
N i=1
j=1

n

≈
αj xk p(xk |y 1:k ) dxk
j=1


=

xk p(xk |y 1:k ) dxk = µk|k

4

(11)

where µk|k is the mean of the ﬁltered PDF p(xk |y 1:k ). In addition, if the merging
weights αj satisfy Eq. (10b), the covariances given by the new ensemble become



1  
(i)
(xk − µk|k )(xk − µk|k )δ xk − xk|k dxk
N i=1

N



(i)
T 1
= xk xk
δ xk − xk|k dxk − µk|k µTk|k
N i=1
N

T

1  (i) (i) Tx x− µk|k µTk|k
N i=1 k|k k|k
N

=

(12)

N 
n


1
(j,i) (j,i) T
αj2 x̂k|k x̂k|k − µk|k µTk|k
N i=1 j=1

n
N




(j,i)
2
T 1
=
αj xk xk
δ xk − x̂k|k dxk − µk|k µTk|k
N i=1
j=1

≈ xk xTk p(xk |y 1:k ) dxk − µk|k µTk|k = Σk|k=where Σk|k is the covariance matrix of p(xk |y 1:k ). Here, we used an approximationasN
1  (j1 ,i)(j ,i)(x̂− µk|k )(x̂k 2 − µk|k )T ≈ 0 (if j1 = j2 ),N i=1 k(j ,1)(j ,N )

(j ,1)(j ,N )

which is justiﬁed because the two sets of samples {x̂k|k1 , · · · , x̂k|k1 } and {x̂k|k2 , · · · , x̂k|k2
are obtained through independent random sampling and would not correlate witheach other. The ensemble obtained using Eq. (9) therefore gives an approximationof p(xk |y 1:k ) preserving the mean and covariances as

1  
(i)
δ xk − xk|k .
N i=1
Np(xk |y 1:k ) ≈

4

(13)

Discussion

The MPF approximately preserves the mean and covariance (i.e., the ﬁrst andsecond moments) of the ﬁltered PDF, and it asymptotically preserves the mean andcovariance as the number of particles approaches inﬁnity. However, the momentsof higher order than the second moment are not preserved and thus the MPF does
not ensure that the shape of the ﬁltered PDF is preserved. It should be noted thatthe moments of higher order than the second moment are not preserved even if the5}number of merged particles n is increased. Here we assume that the dimension of
xk is one for simplicity. The m-th moment of xk is then given as(xk − μk|k )m p(xk |y1:k ) dxk
≈(xk − μk|k )m


1  
(i)
δ xk − xk|k dxk
N i=1
N

N
N m
1  (i)1 m=(x − μk|k ) =N i=1 k|k
N i=1 l=1

n

jl =1

(j ,i)

αjl x̂k|kl − μk|k


N

m



1
(j,i)
(j,i)
m
m
m 1
≈
αj x̂k|k − μk|k=
αj (xk − μk|k )δ xk − x̂k|k dxk
N i=1 j=1
N i=1
j=1

n

m
≈
αj (xk − μk|k )m p(xk |y1:k ) dxk
N 
n


n


j=1

(14)
where we use an approximation as
1 (j ,i)
x̂k|kl − μk|k ≈ 0 (unless j1 = j2 = · · · = jm ).
N i=1 l=1
Nm

In order to preserve the m-th moment even after taking the weighted sum as Eq.
(9), a set of weights {αj }nj=1 must satisfyn

αjm = 1.

(15)

j=1

On the other hand, if the set of real numbers {αj }nj=1 satisfy Eq. (10b), 0 < |αj | < 1
are satisﬁed for all j and thus
|αj |2 > |αj |m(16)

if an integer m is larger than 2. Since Eq. (16) is satisﬁed for all j,
n


αj2 >

j=1

n


αjm .

(17)

j=1

If Eq. (10b) is satisﬁed and m > 2,n

αjm < 1.

j=1

6

(18)

Therefore, if Eq. (10b) is satisﬁed, Eq. (15) can not be satisﬁed when m > 2. Thismeans that there does not exist a set of weights {αj }nj=1 which preserves the second
moment and moments of higher order than the second moment simultaneously.
In implementing the MPF, the number of merged particles n can be chosen
almost arbitrarily. However, n must be equal to or greater than 3. If n = 1, theweight α1 must be 1 in order to satisfy both Eq. (10a) and Eq. (10b), which is
obviously equivalent to the normal PF. If n = 2, then one of merging weights
must be 1, and the other must be 0, so as to satisfy both Eq. (10a) and Eq. (10b).
This setting is also equivalent to the normal PF, which means that the merging
procedure does not make sense. Although there is no upper limit for n, it is notnecessary to set n to be large.
Although the ﬁltered ensemble in the MPF may contain duplicate particles, thenumber of the duplicate particles is much smaller than that for the PF. If theweights are set such that no two weights are equal to each other and that noneof the weights are zero, the number of the duplicate particles can be remarkably
reduced in comparison with the PF. In the MPF, two duplicate particles in the(1)(N )ﬁltered ensemble {xk|k , · · · , xk|k } can be generated only from two identical setsof n particles in the forecast ensemble, if duplicate particles are not contained in(1)(N )the forecast ensemble {xk|k−1, · · · , xk|k−1}. When the probability that a particle(i)xk|k−1 is drawn from the forecast ensemble is wi (0 ≤ wi < 1), the probability
(i )

(i )

n
1
, · · · , xk|k−1} is drawn is nj=1 wij . Thus,
that a sequence of n particles {xk|k−1
the number of duplicate particles contained in the ﬁltered ensemble is, at most,approximately N × (max wi )n for the MPF, while it is N × max wi for the PF.
When n is equal to or greater than 3, there are inﬁnite allowable sets of merging
weights: {α1 , · · · , αn }. Although there is no deﬁnitive way to determine the valuesof the weights, it would be better to set them such that no two weights are equal toeach other and that none of the weights are zero in order to maintain the diversityof the ensemble. However, characteristics of an ensemble obtained as a result of
the MPF depend to some extent on how the set of merging weights is given. For
example, if one of the weights is nearly equal to 1 and the other weights are nearly
equal to zero, the MPF shows a similar behavior to the PF; that is, the shape of aﬁltered PDF is mostly preserved but the degeneration problem tend to be serious(Nakano et al., 2008).

5

Summary

The merging particle ﬁlter (MPF) provides an ensemble approximation of theﬁltered PDF. Each member of a ﬁltered ensemble is generated from a weighted sum

7

of multiple samples from the forecast ensemble such that the mean and covarianceof the ﬁltered distribution are approximately preserved. The MPF can not preserve
the shape of the ﬁltered PDF because the moments of higher order than the second
moment are not preserved. However, when one merging weights is set to be near
1 and the other weights are set to be small, it shows s similar behavior to the PF;
that is, the shape of a ﬁltered PDF is mostly preserved although the degeneration
problem could be serious.

References
Burgers, G., P. J. van Leeuwen, and G. Evensen, 1998., Analysis scheme in the ensemble
Kalman filter, Mon. Wea. Rev., 126 , 1719.
Evensen, G., 1994., Sequential data assimilation with a nonlinear quasi-geostrophic
model using Monte Carlo methods to forecast error statistics, J. Geophys. Res.,
99(C5), 10143.
Goldberg, D. E., Genetic algorithms in search, optimization and machine learning,
Addison-Wesley, Reading, 1989.
Gordon, N. J., D. J. Salmond, and A. F. M. Smith, 1993., Novel approach tononlinear/non-Gaussian Bayesian state estimation, IEE Proceedings F , 140 , 107.
Kitagawa, G., 1993., Monte Carlo filtering and smoothing method for non-Gaussiannonlinear state space model, Inst. Statist. Math. Res. Memo..
Kitagawa, G., 1996., Monte Carlo filter and smoother for non-Gaussian nonlinear state
space models, J. Comp. Graph. Statist., 5 , 1.
Nakano, S., G. Ueno, and T. Higuchi, 2007., Merging particle filter for sequential data
assimilation, Nonlin. Process. Geophys., 14 , 395–408.
Nakano, S., G. Ueno, K. Nakamura, and T. Higuchi, 2008., The merging particle filter
and its characteristics (in Japanese), Proc. of Inst. of Statist. Math., 56 , 225–234.

8