Annecy Del Moral MTT

7/27/2019 Annecy Del Moral MTT

1/48

Particle approximation of multiple object filtering problems

P. Del Moral

INRIA Bordeaux & IMB & CMAP Polytechnique Univ. New South Wales, Sydney (Jan. 2014)

IMACS Monte Carlo, July 15-19, 2013Some hyper-refs

Mean field simulation for Monte Carlo integration. Chapman & Hall, Series : Maths and Stat. (2013).

Particle approximations of a class of branching distribution flows arising in multi-target tracking.SIAM Control. & Opt. (2011). (joint work with Caron, Doucet, Pace)

On the Conditional Distributions of Spatial Point Processes.

Advances in Applied Probability (2011). (joint work with Caron, Doucet, Pace).

On the Stability & the Approximation of Branching Distribution Flows, with Applications to NonlinearMultiple Target Filtering.Stochastic Analysis and Applications (2011). (joint work with Caron, Pace, Vo).

Comparison of implementations of Gaussian mixture PHD filters. FUSION (2010).(joint work with Caron, Pace, Vo).

More references on the website http://www.math.u-bordeaux1.fr/delmoral/index.html [+ Links]
http://hal.archives-ouvertes.fr/docs/00/46/41/30/PDF/RR-7233.pdfhttp://hal.archives-ouvertes.fr/docs/00/46/41/27/PDF/RR-7232.pdfhttp://hal.archives-ouvertes.fr/docs/00/46/41/27/PDF/RR-7232.pdfhttp://hal.inria.fr/docs/00/51/65/07/PDF/RR-7376.pdfhttp://hal.inria.fr/docs/00/51/65/07/PDF/RR-7376.pdfhttp://ieeexplore.ieee.org/xpl/freeabs_all.jsp?reload=true&arnumber=5711953&tag=1http://www.math.u-bordeaux1.fr/~delmoral/index.htmlhttp://www.math.u-bordeaux1.fr/~delmoral/index.htmlhttp://www.math.u-bordeaux1.fr/~delmoral/index.htmlhttp://www.math.u-bordeaux1.fr/~delmoral/index.htmlhttp://ieeexplore.ieee.org/xpl/freeabs_all.jsp?reload=true&arnumber=5711953&tag=1http://hal.inria.fr/docs/00/51/65/07/PDF/RR-7376.pdfhttp://hal.archives-ouvertes.fr/docs/00/46/41/27/PDF/RR-7232.pdfhttp://hal.archives-ouvertes.fr/docs/00/46/41/30/PDF/RR-7233.pdfhttp://find/


2/48

Contents

Introduction/notation

Multiple objects branching signals

Multiple targets filtering models

General measure valued equations

Particle association measures
http://goforward/http://find/http://goback/


3/48

Introduction/notationDefense Industrial Research projectSome basic notationSpatial Branching modelsFirst moments recursion




http://find/http://goback/


4/48

2 Industrial research project

1. Defense industrial Contract : ALEA INRIA team & DCNS (2009)

2. National Research project : ANR PROPAGATION[2,3Me] (2009-2012):

Passive radar tracking and optronics liabilities for the protection of

coastal infrastructures

ALEA INRIA team DCNS SIS, THALES, ECOMER, EXAVISION

Project members : + D. Arrivault, Fr. Caron, M. Pace. Visiting researchers :

D. Clark, A. Doucet, J. Houssineau, S.S. Sing, B.N. Vo.
http://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://www.polemerpaca.com/Securite-et-surete-maritime/Protection-maritime-etendue-et-rapprochee/PROPAGATIONhttp://find/


5/48

Lebesgue integral on a measurable state space E

(, f) =(measure, function) (M(E) Bb(E))

(f) =

(dx) f(x)

Delta-Dirac Measure at a E

= a

a

(f) = f(x) a(dx) = f(a)Normalization of a positive measure M+(E) (when (1) = 0)

(dx) := (dx)/(1) = probability P(E)
http://find/


6/48

Lebesgue integral on a measurable state space E

(, f) =(measure, function) (M(E) Bb(E))

(f) =

(dx) f(x)

Delta-Dirac Measure at a E

= a

a

(f) = f(x) a(dx) = f(a)Normalization of a positive measure M+(E) (when (1) = 0)

(dx) := (dx)/(1) = probability P(E)

Example = 10 Law(X)

(1) = 10 & = Law(X) & (f) = E(f(X))
http://find/


7/48

Q(x, dy) integral operator from E into E

Two operator actions :

f Bb(E) Q(f) Bb(E) and M(E) Q M(E)

with

Q(f)(x) = Q(x, dx)f(x)

[Q](dx) =

(dx)Q(x, dx) ( [Q](f) := [Q(f)] )

and the composition

(Q1Q2)(x, dx) = Q1(x, dx)Q2(x, dx)Q Markov operator Q(1) = 1 =unit function notation : M or K (Markov transitions-kernel/Stochastic matrices)
http://find/


8/48


9/48

Boltzmann-Gibbs transformation : G 0 s.t. (G) > 0

G()(dx) =1

(G)

G(x) (dx)

Bayes rule (with fixed observation y)

(dx) = p(x)dx and G(x) = p(y|x)

G()(dx) =

1

p(y)p(y|x) p(x)dx = p(x|y) dx


10/48


G()(dx) =1

(G)

G(x) (dx)

Bayes rule (with fixed observation y)

(dx) = p(x)dx and G(x) = p(y|x)

G()(dx) =

1

p(y)p(y|x) p(x)dx = p(x|y) dx

Restriction

(dx) = P(X dx) = p(x)dx and G(x) = 1A(x)

G()(dx) =1

P(X A) 1A(x) p(x)dx = P(X dx | X A)
http://find/


11/48


G()(dx) =1

(G)G(x) (dx)

(non unique) Markov transport equation

G()(dy) =

(dx)S(x, dy) G() = S

Example 1 : (G 1) accept/reject/recycling/interacting jumps

S(x, dy) = G(x)x(dy) + (1 G(x)) G()(dy)
http://find/


12/48


G()(dx) =1

(G)G(x) (dx)

(non unique) Markov transport equation

G()(dy) =

(dx)S(x, dy) G() = S

Example 1 : (G 1) accept/reject/recycling/interacting jumps

S(x, dy) = G(x)x(dy) + (1 G(x)) G()(dy)

Note :

S1N

1jN Xj

(Xi, dy)

= G(Xi)Xi(dy) + (1 G(Xi))

1jNG(Xj)

1kN G(Xk)

Xj(dy)
http://find/


13/48

Other examples of transport equations G() = SExample 2 : s.t. G 1 a.e.

S(x, dy) = G(x) x(dy) + (1 G(x)) G()(dy)

Example 3 : a s.t. G > a

S

(x, dy) =a

(G)

x(dy) + 1

a

(G) Ga()(dy)Example 4 : G

S(x, dy) = (x) x(dy) + (1 (x)) [GG(x)]+ ()(dy)

with the acceptance rate

(x) = [G G(x)]/(G)
http://find/


14/48

Spatial Branching models (time index n N, state spaces En) 3 ingredients : Gn(x)

1, n(dx)

0, and Mn(xn1, dxn) Markov.

Branching rule (spawning) :

Random mapping x gn(x) N offsprings, with E(gn(x)) = Gn(x)

survival rates en(x) + cemetery states : Gn en(x)Gn(x)

Spontaneous births: Spatial Poisson with intensity n(dx) Free motion between branching times : Mn-evolutions

Random occupation measure (after the n-th evolution step)

Xn =Nn

i=1

Xin

En := {types, locations, labels, excursions, paths,. . . }

S i l B hi d l ( i i d N E )
http://find/


15/48

Spatial Branching models (time index n N, state spaces En) First moment recursion = branching intensity distribution

n+1(f) := E (X

n+1(f)) = n(Qn+1(f)) + n+1(f)

withQn+1(x, dy) = Gn(x)Mn+1(x, dy)

S ti l B hi d l (ti i d N t t E )
http://find/


16/48

Spatial Branching models (time index n N, state spaces En) First moment recursion = branching intensity distribution

n+1(f) := E (X

n+1(f)) = n(Qn+1(f)) + n+1(f)

withQn+1(x, dy) = Gn(x)Mn+1(x, dy)

Sketched proof (n = 0): Xn+1 =Nn+1i=1

Xin+1 =

Nni=1

gin(Xin)

j=1

Xi,jn+1

E (Xn+1(f) | Xn, gn(Xn)) =Nn

i=1

gi

n(Xi

n) Mn+1(f)(Xi

n)

E (

Xn+1(f)

| Xn) =

Nn

i=1Gn(X

in) Mn+1(f)(X

in) =

Xn(Qn+1(f))
http://find/


17/48

Continuous time models

Geometric clocks exponential rates time mesh(tn tn1) 0

Xn = XtnGn = survival

[spawning

mean offsprings + (1

spawning)

1]
http://find/


18/48

Continuous time models

Geometric clocks exponential rates time mesh(tn tn1) 0

Xn = XtnGn = survival

[spawning

mean offsprings + (1

spawning)

1]

(tn tn1) 0 G = 1 + V dt and M = Id + L dt and tn t

d

dtt(f) = t(L

V(f)) + t(f) with LV = L + V

Schrodinger operator
http://find/


19/48

n = 0 Classical Feynman-Kac models

Feynman-Kac representation ( Application domains)

n+1(f) = 0(1) E0

f(Xn+1)

0pnGp(Xp)

Particle approximations = Genetic type algo = Particle filters = ...

Qn+1(x, dy) = Gn(x) Selection potential

Mn+1(x, y)

Mutation transition
http://find/


20/48


Multiple objects branching signalsEvolution equationsStability propertiesThree typical scenariosAn extended Feynman-Kac modelMean field particle interpretationsSome convergence results



http://find/


21/48

More general spatial Branching models (hyp. 0 = 0)

n = n1Qn + n and n := n/n(1)

(n(1), n) := p,n (p(1), p)

Some problems Problem 1: Mass n(1) unstable n(1) or n(1) 0 as n Problem 2: Xn =

Nni=1 Xin generally NOT POISSON random field.

Problem 3:

non degenerate numerical sampling method?

Problem 4: non degenerate approximation ofn?

Th i M M G G [ ]


22/48

Three scenarios Mn = M, Gn = G [g, g+], n = 1. G = 1 := M (independent of)

n(1) = 0(1) + n(1) and n tv = O(1/n)2. g+ < 1 := /(1) with given by

:=n0

Qn Poisson equation (Id Q) =

and|n(f) (f)| |n(f) (f)| c gn+ f

Th i M M G G [ ]
http://find/


23/48

Three scenarios Mn = M, Gn = G [g, g+], n = 1. G = 1 := M (independent of)

n(1) = 0(1) + n(1) and n tv = O(1/n)2. g+ < 1 := /(1) with given by

:=n0

Qn Poisson equation (Id Q) =

and|n(f) (f)| |n(f) (f)| c gn+ f

Continuous time models G = eVt

& M = L t + Id

t(f) =

t0

E

f(Xs) exp

s

0

V(Xr)dr

ds

t

Poisson equation L

V = , with LV = L + V

The 3 d sce a io (M M G G [ ] )
http://find/


24/48

The 3-rd scenario (Mn = M, Gn = G [g, g+], n = )

g > 1

(f) := Q(f)/Q(1) (independent of)

limn

1

nlog n(1) = log (G) and n tv c en

= [quasi-invariant meas., Yaglom meas., ground states, Feynman-Kacsg fixed points, infinite population stationary measure, . . . ]

Hyper-refs :

On the stability of interacting processes with applications to filtering and geneticalgorithms. (joint work with A. Guionnet) Annales IHP (2001).

Particle approximations of Lyapunov exponents connected to Schrodinger operators andFeynman-Kac semigroups. (joint work with L. Miclo) ESAIM: P&S (2003).

Particle Motions in Absorbing Medium with Hard and Soft Obstacles.

(joint work with A. Doucet) Stochastic Analysis and Applications (2004).

Nonlinear equations
http://www.math.u-bordeaux1.fr/~pdelmora/ihp.pshttp://www.math.u-bordeaux1.fr/~pdelmora/ihp.pshttp://www.math.u-bordeaux1.fr/~pdelmora/exponent.pshttp://www.math.u-bordeaux1.fr/~pdelmora/exponent.pshttp://www.math.u-bordeaux1.fr/~pdelmora/exponent.pshttp://www.math.u-bordeaux1.fr/~pdelmora/obstacle.pshttp://www.math.u-bordeaux1.fr/~pdelmora/obstacle.pshttp://www.math.u-bordeaux1.fr/~pdelmora/exponent.pshttp://www.math.u-bordeaux1.fr/~pdelmora/exponent.pshttp://www.math.u-bordeaux1.fr/~pdelmora/ihp.pshttp://www.math.u-bordeaux1.fr/~pdelmora/ihp.pshttp://find/


25/48

Nonlinear equations

n+1 n(1) nQn+1 + n+1(1) n+1

Nonlinear & interacting mass + proba measures equations

n+1(1) = n(1) n(Gn) + n+1(1)

n+1 = Gn (n)Mn+1,(n(1),n)

with the Markov transitions:

Mn+1,(m,)(x, dy) := n (m, ) Mn+1(x, dy) + (1 n (m, )) n+1(dy)

and the collection of [0, 1]-parameters

n (m, ) =m (Gn)

m (Gn) + n+1(1)

An extended Feynman Kac model
http://find/


26/48

An extended Feynman-Kac model

nupdating n := Gn (n) = nSn,n prediction n+1 := nMn+1,(n(1),n)

A couple of equations:

The total mass evolution

n+1(1) = n(1) n(Gn) + n+1(1)

The nonlinear filtering/Feynman-Kac type conservative equations

n+1 = nSn,n Mn+1,(n(1),n) := n Kn+1,(n(1),n) Markov transition
http://find/


27/48

Mean field interacting particle models

Nn =1

N 1iNin N n and Nn (1) N n(1)

the total mass evolution [deterministic]

Nn+1(1) := Nn (1)

Nn (Gn) + n+1(1)

Mean field particle model

in+1 = r.v. with distribution Kn+1,(Nn (1),Nn )(in, dxn+1)

(Local) Stochastic perturbation model:

Nn+1 := Nn Kn+1,(Nn (1),Nn ) +

1N

WNn+1

Theoretical convergence results
http://find/


28/48

g

Independent local sampling error fluctuations

(WNn )n0

N iid centered Gaussian fields (Wn)n0

http://find/


29/48


(WNn )n0


Functional CLT(s) (withNn :=

Nn (1) Nn

)

V,Nn :=

N (Nnn) & V

,Nn :=

N (Nn

n)

N V

n & V

n

http://find/


30/48


(WNn )n0



Nn (1) Nn

)

V,Nn :=

N (Nnn) & V

,Nn :=

N (Nn

n)

N V

n & V

n

Uniform cv results (under some mixing conditions on Mn)

supn0

ENn

n (f)

p

c(p)/Np/2 (

uniform concentration)

http://find/


31/48


(WNn )n0



Nn (1) Nn

)

V,Nn :=

N (Nnn) & V

,Nn :=

N (Nn

n)

N V

n & V

n

Uniform cv results (under some mixing conditions on Mn)

supn0

ENn

n (f)

p

c(p)/Np/2 (

uniform concentration)

Unbiased particle total mass with variance (N n)

E1 Nn (1)/n(1)

2

c n/N
http://find/


32/48



Multiple targets filtering modelsConditioning principlesPHD filtering equationStability properties



Conditioning principles for marked point processes
http://find/


33/48

Poisson point process X with intensity (dx1) Q(x1, dx2) onE = (E1 E2)

X:= mN(X1,X2) = 1iN

(Xi1,Xi2) and Xj := mN(Xj) = 1iN

Xij

http://find/


34/48


X:= mN(X1,X2) = 1iN


Xij

2 Bayes rules: Normalization p(x2|x1) Markov operator p(x1|x2)

Q(x1, dx2) =Q(x1, dx2)

Q(x1,E2)and (dx1) Q(x1, dx2) = (Q) (dx2) Q(x2, dx1)

http://find/


35/48


X:= mN(X1,X2) = 1iN


Xij

2 Bayes rules: Normalization p(x2|x1) Markov operator p(x1|x2)

Q(x1, dx2) =Q(x1, dx2)Q(x1,E2)

and (dx1) Q(x1, dx2) = (Q) (dx2) Q(x2, dx1)

2 conditional distributions formulae:

E (F1(X1) | X2 ) = F1 (mN(x1)) 1iN

Q(Xi2, dx

i1)

E (F2(X2) | X1 ) =

F2 (mN(x2))

1iNQ(Xi1, dx

i2)

http://find/


36/48

g p p p p

(X1,X2) = (X,Y), X Poisson Signal (dx) YPoisson Obs.

(Xi = x) (Yi = y) (x) g(x, y) (dy) + (1 (x)) c(dy)

Clutter Y Poisson with intensity (dy) = h(y) (dy)

http://find/


37/48

(X1,X2) = (X,Y), X Poisson Signal (dx) YPoisson Obs.

(Xi = x) (Yi = y) (x) g(x, y) (dy) + (1 (x)) c(dy)

Clutter Y Poisson with intensity (dy) = h(y) (dy)

Observables Y0 = Y 1=c ( = detection rate)

(f) := E (X(f) | Yo)= ((1 )f) +

Yo(dy) (1 (y)) g(y,.)()(f)

with the conditional clutter probability density

(y) = h(y)/[h(y) + (g(y, .))]

Ex.: full detect and no clutter = 1 & h = 0 Yo = Y
http://find/


38/48

Conditional mean number of targets and their distributions

(1) = Y(1)and

(f) := (f)

(1) = Y(dy) =Y/Y(1)g(y,.)()(f) Bayes rule

with := /(1)

Single target Yo

= Y Classical filtering updating equations = g(Y,.)()

PHD filt i ti [Si l b hi d l (Q )]


39/48

PHD filtering equation [Signal branching model (Qn, n)]

Hyp.: Xn+1 Poisson n+1 = nQn + n with obs. Y0n+1 as before

= PHD filtering equations:

n+1 := nQn + nn(f) := n((1 n)f) + Yon (dy) (1 n (y)) ngn (y,.)(n)(f)

A class of measure valued equations PHD; Bernoulli filters, etc.

n+1 = nQn+1,n
http://find/


40/48

Stability properties of meas. valued equations

n = n/n(1) Nonlinear semigroup (n(1), n) = p,n(p(1), p)

Stability Theorem :

p,n(m

, )

p,n(m, )

c e(np)

Regularity prop. 3 natural conditions on the PHD filter/model

1. small clutter intensities

2. high detection probability

3. high spontaneous birth rates
http://find/


41/48




General measure valued equationsNonlinear evolution equationsMean field particle approximation

http://find/


42/48

Mean field particle models


43/48


Nn =1


with

Nn+1(1) = Nn (1) Nn (Gn,Nn (1)Nn )

in+1 = random var. with law Kn+1,(Nn (1)Nn )(in, dx)

Same theorems as before with uniform convergence estimates

http://find/


44/48


Nn =1


with

Nn+1(1) = Nn (1) Nn (Gn,Nn (1)Nn )

in+1 = random var. with law Kn+1,(Nn (1)Nn )(in, dx)

Same theorems as before with uniform convergence estimates

Abstract general models

numerical scheme with local errors Interacting Kalman type filters particle associations measures ( GM-PHD)
http://find/


45/48






Association measures [ = 1 & h = 0 & Qn = Mn]
http://find/


46/48

Ex. : Computable (exact or approximate) filters

The mappings ynn+1() := gn(yn,.

)()Mn+1

{Kalman, EKF, Ensemble Kalman filters, particle filters,. . . }

Initial association measure

1 :=

Y0(dy0) y01 (0) N1 :=

YN0 (dy0) y01 (0)

for instance

YN0 =1

N

1iN

Yi0 i.i.d. samples from Y0 or (if possible) YN

0 = Y0

Particle association measures [ = 1 & h = 0 & Qn = Mn]
http://find/


47/48

2 Y1(dy1)

y12 (

N1 )

=

Y1(dy1) YN0 (dy0)

y01 (0)(g1(y1, .))YN0 (dy0) y01 (0)(g1(y1, .)) [y12 y01 ] (0)

YN

0,1(d(y0, y1)) [y1

2 y0

1

] (0)

for instance

YN0,1 =1

N

1iN

(Yi0,1,Yi1,1)

i.i.d. samples from the N Y1(1) supported measures

Y1(dy1) YN0 (dy0)y01 (0)(g1(y1, .))YN0 (dy0) y01 (0)(g1(y1, .)) (y0,y1)

and so on ...

Particle association measures - Track management
http://find/


48/48

Association particle tree genealogies

Nn+1

:= YN0,n(d(y0, . . . , yn)) ynn+1 . . . y01 (0)with

YN0,n :=1

N

1iN

(Yi0,n,Yi1,n,...,Yin,n)

Stochastic models and cv analysis :

General case :

(miss-detect, survival, spontaneous birth) = as before virtual obs.

Abstract models of the form n+1 = nQn+1,n . Mean field particle models Association particle measures. Lp-bounds Concentration sub-Gaussian inequalities.
http://find/

Documents

Annecy Del Moral MTT