Theory And Implementation Of Particle Filters

Theory and Implementation of Particle Filters Miodrag Bolic Assistant Professor School of Information Technology and Engineering University of Ottawa [email protected] 12 Nov 2004

1

Big picture Observed signal 1

sensor

t Observed signal 2

Particle Filter

Estimation

t

‰ ‰

Goal: Estimate a stochastic process given some noisy observations Concepts: ‰ ‰

12 Nov 2004

t

Bayesian filtering Monte Carlo sampling 2

Particle filtering operations „

„

„

Particle filter is a technique for implementing recursive Bayesian filter by Monte Carlo sampling The idea: represent the posterior density by a set of random particles with associated weights. Compute estimates based on these samples and weights Posterior density

Sample space 12 Nov 2004

3

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

4

Motivation „

„

„ „

The trend of addressing complex problems continues Large number of applications require evaluation of integrals Non-linear models Non-Gaussian noise

12 Nov 2004

5

Sequential Monte Carlo Techniques „ „ „ „ „

Bootstrap filtering The condensation algorithm Particle filtering Interacting particle approximations Survival of the fittest

12 Nov 2004

6

History First attempts – simulations of growing polymers

„

M. N. Rosenbluth and A.W. Rosenbluth, “Monte Carlo calculation of the average extension of molecular chains,” Journal of Chemical Physics, vol. 23, no. 2, pp. 356–359, 1956.

‰

First application in signal processing - 1993

„

N. J. Gordon, D. J. Salmond, and A. F. M. Smith, “Novel approach to nonlinear/nonGaussian Bayesian state estimation,” IEE Proceedings-F, vol. 140, no. 2, pp. 107–113, 1993.

‰

Books

„

A. Doucet, N. de Freitas, and N. Gordon, Eds., Sequential Monte Carlo Methods in Practice, Springer, 2001. B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House Publishers, 2004.

‰

‰

Tutorials

„ ‰

M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-gaussian Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002.

12 Nov 2004

7

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

8

Applications „

Signal processing ‰

‰ ‰

„

Image processing and segmentation Model selection Tracking and navigation

Communications ‰ ‰ ‰

Channel estimation Blind equalization Positioning in wireless networks

„

Other applications1) ‰ ‰ ‰ ‰ ‰ ‰ ‰

‰ ‰

1)

Biology & Biochemistry Chemistry Economics & Business Geosciences Immunology Materials Science Pharmacology & Toxicology Psychiatry/Psychology Social Sciences

A. Doucet, S.J. Godsill, C. Andrieu, "On Sequential Monte Carlo Sampling Methods for Bayesian Filtering", Statistics and Computing, vol. 10, no. 3, pp. 197-208, 2000

12 Nov 2004

9

Bearings-only tracking „

„

„ „ „

The aim is to find the position and velocity of the tracked object. The measurements taken by the sensor are the bearings or angles with respect to the sensor. Initial position and velocity are approximately known. System and observation noises are Gaussian. Usually used with a passive sonar.

12 Nov 2004

10

Bearings-only tracking ‰ States: position and velocity xk=[xk, Vxk, yk, Vyk]T ‰ Observations: angle

zk

‰ Observation equation:

zk=atan(yk/ xk)+vk

‰ State equation:

xk=Fxk-1+ Guk y

yk yk+1 Trajectory zk

zk+1 xk xk+1

x

12 Nov 2004

11

Bearings-only tracking

12 Nov 2004

‰

Blue – True trajectory

‰

Red – Estimates

12

Car positioning „

„

„

„

„

Observations are the velocity and turn information1) A car is equipped with an electronic roadmap The initial position of a car is available with 1km accuracy In the beginning, the particles are spread evenly on the roads As the car is moving the particles concentrate at one place

1) Gustafsson et al., “Particle Filters for Positioning, Navigation, and Tracking,” IEEE Transactions on SP, 2002

12 Nov 2004

13

Detection over flat-fading channels „

„

„

Detection of data transmitted over unknown Rayleigh fading channel The temporal correlation in the channel is modeled using AR(r) process At any instant of time t, the unknowns are , and , and our main objective is to detect the transmitted symbol sequentially h(t) st

g(t)

y(t)

s(t) Channel

12 Nov 2004

v(t) yt Sampling

14

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

15

Fundamental concepts „

„ „ „

State space representation Bayesian filtering Monte-Carlo sampling Importance sampling

State space model Solution

Problem

Estimate posterior

Integrals are not tractable

Monte Carlo Sampling

Difficult to draw samples

Importance Sampling

12 Nov 2004

16

Representation of dynamic systems ‰

The state sequence is a Markov random process

State equation: xk=fx(xk-1, uk) ‰ ‰ ‰

xk state vector at time instant k fx state transition function uk process noise with known distribution

Observation equation: zk=fz(xk, vk) ‰ ‰ ‰

zk observations at time instant k fx observation function vk observation noise with known distribution

12 Nov 2004

17

Representation of dynamic systems The alternative representation of dynamic system is by densities.

„

State equation: p(xk|xk-1) Observation equation: p(zk|xk)

„

The form of densities depends on:

„

‰ ‰

Functions fx(·) and fz(·) Densities of uk and vk

12 Nov 2004

18

Bayesian Filtering „

The objective is to estimate unknown state xk, based on a sequence of observations zk, k=0,1,… . Objective in Bayesian approach ↓ Find posterior distribution p(x0:k|z1:k)

„

By knowing posterior distribution all kinds of estimates can be computed:

12 Nov 2004

19

Update and propagate steps „ „

k=0 Bayes theorem Filtering density: Predictive density:

p (x1 | z 0 ) = ∫ p (x1 | x 0 ) p (x 0 | z 0 )dx 0

z1

z0 p(x0)

p(z 0 | x 0 ) p(x 0 ) p(z 0 )

p(x 0 | z 0 ) =

Update

Propagate p(x0|z0)

12 Nov 2004

p(x1|z0)

Update

z2 Propagate p(x1|z1)

…

p(x2|z1)

Update p(xk|zk-1)

Propagate p(xk|zk)

p(xk+1|zk)

20

Update and propagate steps „ „

k>0 Derivation is based on Bayes theorem and Markov property Filtering density:

Predictive density:

p (x k | z1:k ) =

p (z k | x k ) p (x k | z1:k −1 ) p (z k | z1:k −1 )

p(x k +1 | z1:k ) = ∫ p (x k +1 | x k ) p (x k | z1:k )dx k

12 Nov 2004

21

Meaning of the densities Bearings-only tracking problem „ p(xk|z1:k) posterior ‰

„

p(xk|xk-1) prior ‰

„

What is the probability that the object is at the location xk for all possible locations xk if the history of measurements is z1:k? The motion model – where will the object be at time instant k given that it was previously at xk-1?

p(zk|xk) likelihood ‰

The likelihood of making the observation zk given that the object is at the location xk.

12 Nov 2004

22

Bayesian filtering - problems „ „

„

„

Optimal solution in the sense of computing posterior The solution is conceptual because integrals are not tractable Closed form solutions are possible in a small number of situations Gaussian noise process and linear state space model ↓ Optimal estimation using the Kalman filter Idea: use Monte Carlo techniques

12 Nov 2004

23

Monte Carlo method „

Example: Estimate the variance of a zero mean Gaussian process +∞ v = ∫ x 2 p ( x)dx −∞

„

1.

Monte Carlo approach: Simulate M random variables from a Gaussian distribution

x ( m ) ~ N (0, σ 2 ) 2.

Compute the average

v= 12 Nov 2004

1 M

∑

M

(m) 2 ( x ) m =1

24

Importance sampling Classical Monte Carlo integration – Difficult to draw samples from the desired distribution Importance sampling solution:

„

„

Draw samples from another (proposal) distribution Weight them according to how they fit the original distribution

1. 2.

Free to choose the proposal density Important:

„ „

It should be easy to sample from the proposal density Proposal density should resemble the original density as closely as possible

‰ ‰

12 Nov 2004

25

Importance sampling „

Evaluation of integrals E ( f ( X )) = ∫ f ( x) p ( x)dx = ∫ f ( x) X

X

p( x) π ( x)dx π ( x)

Monte Carlo approach: 1. Simulate M random variables from proposal density π(x) „

x ( m ) ~ π ( x) 2. Compute the average 1 E ( f ( x)) ≈ M

∑

M m =1

f (x

(m)

p( x ( m) ) ) π ( x (m) ) 1 424 3 w( m )

12 Nov 2004

26

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

27

Sequential importance sampling Idea: „ Update filtering density using Bayesian filtering „ Compute integrals using importance sampling „

The filtering density p(xk|z1:k) is represented using

{x

particles and their weights „

Compute weights using:

(m) k

w

( m) k

, wk( m )

}

M

m =1

p ( xk( m ) , z1:k ) = π ( xk( m ) , z1:k )

Posterior

x 12 Nov 2004

28

Sequential importance sampling „ „

Let the proposal density be equal to the prior Particle filtering steps for m=1,…,M:

1. Particle generation

xk( m ) ~ p ( xk | xk −1 )

2a. Weight computation

wk*( m ) = wk*(−m1 ) p ( z k | xk( m ) ) (m) k

w

2b. Weight normalization

=

wk*( m ) M

∑w m =1

*( m ) k M

E ( g ( xk | z1:k )) = ∑ g ( xk( m ) ) wk( m )

3. Estimate computation

m =1

12 Nov 2004

29

Resampling Problems: ‰

Weight Degeneration

‰

Wastage of computational resources Solution ⇒ RESAMPLING

‰

Replicate particles in proportion to their weights

‰

Done again by random sampling M

⎧ 1⎫ (m) (m) x , k ⎨ ⎬ ~ xk , wk M ⎭ m =1 ⎩ ~ (m)

12 Nov 2004

{

}

M

m =1

30

Resampling M

⎧ (m) 1 ⎫ ⎨ xk + 2 , ⎬ M ⎭ m =1 ⎩

M

⎧ ~ (m) 1 ⎫ ⎨ x k +1 , ⎬ M ⎭m =1 ⎩

{x

( m) k +1

, wk( m+1)

}

M

m =1

M

⎧ (m) 1 ⎫ ⎨ xk +1 , ⎬ M ⎭ m =1 ⎩

M

⎧ ~ (m) 1 ⎫ ⎨xk , ⎬ M ⎭m =1 ⎩

{x

(m) k

, wk( m )

}

M

m =1

M

⎧ (m) 1 ⎫ ⎨ xk −1 , ⎬ M ⎭ m =1 ⎩ x 12 Nov 2004

31

Particle filtering algorithm Initialize particles New observation Particle generation 1

2

... M

1

2

... M

Weigth computation Normalize weights

Output estimates Resampling Output yes

More observations? no Exit

12 Nov 2004

32

Bearings-only tracking example MODEL „ States: xk=[xk, Vxk, yk, Vyk]T „

„

Observations: zk Noise u k ~N( 0 ,σ u2 ), vk ~N( 0 ,σ v2 )

„

„

ALGORITHM „ Particle generation

State equation: xk=Fxk-1+ Guk Observation equation: zk=atan(yk/ xk)+vk

‰

2 Generate M random numbers u(m) k ~N( 0 ,σ u )

‰

Particle computation

x

(m) k

~ (m)

= F x k + Gu(m) k

„

Weight computation yk( m ) 2 *( m ) wk = Ν ( z k − atan ( m ) , σ v ) xk Weight normalization

„

Resampling

„

„

M

⎧ ~ (m) 1 ⎫ (m) (m) ⎨x k , ⎬ ~ x k , wk M ⎭m =1 ⎩

{

}

M

m =1

Computation of the estimates

12 Nov 2004

33

Bearings-Only Tracking Example

12 Nov 2004

34

Bearings-Only Tracking Example

12 Nov 2004

35

Bearings-Only Tracking Example

12 Nov 2004

36

General particle filter „

„

If the proposal is a prior density, then there can be a poor overlap between the prior and posterior Idea: include the observations into the proposal density π ( xk | xk −1 , zk ) = p( xk | x ( m ) , z k ) k −1

„

This proposal density minimize Var ( wk*(m ) )

12 Nov 2004

37

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

38

Advantages of particle filters „ „

„ „

Ability to represent arbitrary densities Adaptive focusing on probable regions of state-space Dealing with non-Gaussian noise The framework allows for including multiple models (tracking maneuvering targets)

12 Nov 2004

39

Disadvantages of particle filters „ „

„

„

„

High computational complexity It is difficult to determine optimal number of particles Number of particles increase with increasing model dimension Potential problems: degeneracy and loss of diversity The choice of importance density is crucial

12 Nov 2004

40

Variations Rao-Blackwellization: „

„

Some components of the model may have linear dynamics and can be well estimated using a conventional Kalman filter. The Kalman filter is combined with a particle filter to reduce the number of particles needed to obtain a given level of performance.

12 Nov 2004

41

Variations Gaussian particle filters „

„

„ „

Approximate the predictive and filtering density with Gaussians Moments of these densities are computed from the particles Advantage: there is no need for resampling Restriction: filtering and predictive densities are unimodal

12 Nov 2004

42

Outline „ „ „ „ „ „

Motivation Applications Fundamental concepts Sample importance resampling Advantages and disadvantages Implementation of particle filters in hardware

12 Nov 2004

43

Challenges and results Challenges ‰ ‰ ‰

Reducing computational complexity Randomness – difficult to exploit regular structures in VLSI Exploiting temporal and spatial concurrency

Results ‰

‰ ‰

New resampling algorithms suitable for hardware implementation Fast particle filtering algorithms that do not use memories First distributed algorithms and architectures for particle filters

12 Nov 2004

44

Complexity Complexity

Initialize particles New observation Particle generation

1

2

... M

1

2

... M

4M random number generations

M exponential and arctangent functions

Weigth computation Normalize weights

Output estimates Resampling

Propagation of the particles Output

More observations?

yes

Bearings-only tracking problem Number of particles M=1000

no Exit

12 Nov 2004

45

Mapping to the parallel architecture Start New observation Particle generation 1

2

... M

1

2

... M

Processing Element 1

Central Unit

Weight computation Resampling Propagation of particles

Processing Element 3

Processing elements (PE) 9 Particle generation 9 Weight Calculation

Exit 12 Nov 2004

Processing Element 2

Processing Element 4

Central Unit 9 Algorithm for particle propagation 9 Resampling 46

Propagation of particles p PE 1

PE 2

PE 3

PE 4

Particles after resampling

Disadvantages of the particle propagation step ‰ Random communication pattern ‰ Decision about connections is not known before the run time Processing Element 1

Processing Element 2

‰ Requires dynamic type of a network ‰ Speed-up is significantly affected

Central Unit Processing Element 3

Processing Element 4

12 Nov 2004

47

Parallel resampling N=0 4

1

4

‰

‰

N=0

2

1

4

N=8

N=4

2

1

1

N=4 1

2

1

3

4

3

N=0

N=3

N=0

4

4

3

N=8

N=4

1 1

4 N=4

Solution ‰

‰

N=13

The way in which Monte Carlo sampling is performed is modified

Advantages ‰

Propagation is only local

‰

Propagation is controlled in advance by a designer

‰

Performances are the same as in the sequential applications

Result ‰

Speed-up is almost equal to the number of PEs (up to 8 PEs)

12 Nov 2004

48

Architectures for parallel resampling ‰ ‰

Controlled particle propagation after resampling Architecture that allows adaptive connection among the processing elements PE1

PE3

Central Unit

PE2

PE4

12 Nov 2004

49

Space exploration ‰ ‰ ‰

Hardware platform is Xilinx Virtex-II Pro Clock period is 10ns PFs are applied to the bearings-only tracking problem Limit: Available memory

1000

Sample period (us)

Limit: Logic blocks 1

100

Number of particles M

Virtex II Pro design space

500 2

1000 5000

4

10000 50000

8 10 16 32

1 1

10

K=14

100

Number of PEs

12 Nov 2004

50

Summary „

„

Very powerful framework for estimating parameters of non-linear and non-Gaussian models Main research directions ‰ ‰

‰

„

Finding new applications for particle filters Developing variations of particle filters which have reduced complexity Finding the optimal parameters of the algorithms (number of particles, divergence tests)

Challenge ‰

Popularize the particle filter so that it becomes a standard tool for solving many problems in industry

12 Nov 2004

51