11 Feedback Particle Filter
Tao Yang, Prashant G. Mehta and Sean P. Meyn
Abstract
A new formulation of the particle filter for nonlinear filtering is presented, based on concepts fromoptimal control, and from the mean-field game theory. The optimal control is chosen so that the posteriordistribution of a particle matches as closely as possible the posterior distribution of the true state giventhe observations. This is achieved by introducing a cost function, defined by the Kullback-Leibler (K-L)divergence between the actual posterior, and the posterior of any particle.The optimal control input is characterized by a certain Euler-Lagrange (E-L) equation, and is shownto admit an innovation error-based feedback structure. For diffusions with continuous observations, thevalue of the optimal control solution is ideal. The two posteriors match exactly, provided they areinitialized with identical priors. The feedback particle filter is defined by a family of stochastic systems,each evolving under this optimal control law.A numerical algorithm is introduced and implemented in two general examples, and a neuroscienceapplication involving coupled oscillators. Some preliminary numerical comparisons between the feed-back particle filter and the bootstrap particle filter are described.
I. I
NTRODUCTION
We consider a scalar filtering problem:d X t = a ( X t ) d t + σ B d B t , (1a)d Z t = h ( X t ) d t + σ W d W t , (1b)where X t ∈ R is the state at time t , Z t ∈ R is the observation process, a ( · ) , h ( · ) are C functions,and { B t } , { W t } are mutually independent standard Wiener processes. Unless otherwise noted,the stochastic differential equations (SDEs) are expressed in Itˆo form.The objective of the filtering problem is to compute or approximate the posterior distributionof X t given the history Z t : = σ ( Z s : s ≤ t ) . The posterior p ∗ is defined so that, for any measurableset A ⊂ R , (cid:90) x ∈ A p ∗ ( x , t ) d x = P { X t ∈ A | Z t } . (2) T. Yang and P. G. Mehta are with the Coordinated Science Laboratory and the Department of Mechanical Sci-ence and Engineering at the University of Illinois at Urbana-Champaign (UIUC) [email protected];[email protected]
S. P. Meyn is with the Department of Electrical and Computer Engineering at the University of Florida [email protected]
Financial support from the AFOSR grant FA9550-09-1-0190 and the NSF grant EECS-0925534 is gratefully acknowledged.The conference version of this paper appeared in [25], [24].
February 27, 2013 DRAFT a r X i v : . [ m a t h . NA ] F e b The filter is infinite-dimensional since it defines the evolution, in the space of probabilitymeasures, of { p ∗ ( · , t ) : t ≥ } . If a ( · ) , h ( · ) are linear functions, the solution is given by thefinite-dimensional Kalman filter. The theory of nonlinear filtering is described in the classicmonograph [15].The article [3] surveys numerical methods to approximate the nonlinear filter. One approachdescribed in this survey is particle filtering.The particle filter is a simulation-based algorithm to approximate the filtering task [13], [11],[8]. The key step is the construction of N stochastic processes { X it : 1 ≤ i ≤ N } . The value X it ∈ R is the state for the i th particle at time t . For each time t , the empirical distribution formed by,the “particle population” is used to approximate the conditional distribution. Recall that this isdefined for any measurable set A ⊂ R by, p ( N ) ( A , t ) = N N ∑ i = { X it ∈ A } . (3)A common approach in particle filtering is called sequential importance sampling , whereparticles are generated according to their importance weight at every time stage [8], [3]. Bychoosing the sampling mechanism properly, particle filtering can approximately propagate theposterior distribution, with the accuracy improving as N increases [5].The objective of this paper is to introduce an alternative approach to the construction of aparticle filter for (1a)-(1b) inspired by mean-field optimal control techniques; cf., [14], [26]. Inthis approach, the model for the i th particle is defined by a controlled system,d X it = a ( X it ) d t + σ B d B it + d U it , (4)where X it ∈ R is the state for the i th particle at time t , U it is its control input, and { B it } are mutually independent standard Wiener processes. Certain additional assumptions are maderegarding admissible forms of control input.Throughout the paper we denote conditional distribution of a particle X it given Z t by p where,just as in the definition of p ∗ : (cid:90) x ∈ A p ( x , t ) d x = P { X it ∈ A | Z t } . (5)The initial conditions { X i } Ni = are assumed to be i.i.d., and drawn from initial distribution p ∗ ( x , ) of X (i.e., p ( x , ) = p ∗ ( x , ) ).The control problem is to choose the control input U it so that p approximates p ∗ , andconsequently p ( N ) (defined in (3)) approximates p ∗ for large N . The synthesis of the controlinput is cast as an optimal control problem, with the Kullback-Leibler metric serving as the costfunction. The optimal control input is obtained via analysis of the first variation. February 27, 2013 DRAFT
The main result of this paper is to derive an explicit formula for the optimal control input,and demonstrate that under general conditions we obtain an exact match: p = p ∗ under optimalcontrol. The optimally controlled dynamics of the i th particle have the following Itˆo form,d X it = a ( X it ) d t + σ B d B it + K ( X it , t ) d I it + Ω ( X it , t ) d t (cid:124) (cid:123)(cid:122) (cid:125) optimal control, d U i ∗ t , (6)in which Ω ( x , t ) : = σ W K ( x , t ) K (cid:48) ( x , t ) , K (cid:48) ( x , t ) = ∂ K ∂ x ( x , t ) , and I i is similar to the innovationprocess that appears in the nonlinear filter,d I it : = d Z t − ( h ( X it ) + ˆ h t ) d t , (7)where ˆ h t : = E [ h ( X it ) | Z t ] = (cid:82) h ( x ) p ( x , t ) d x . In a numerical implementation, we approximateˆ h t ≈ ˆ h ( N ) t : = N N ∑ i = h ( X it ) . (8)The gain function K is shown to be the solution to the following Euler-Lagrange boundaryvalue problem (E-L BVP): − ∂∂ x (cid:18) p ( x , t ) ∂∂ x { p ( x , t ) K ( x , t ) } (cid:19) = σ W h (cid:48) ( x ) , (9)with boundary conditions lim x →± ∞ p ( x , t ) K ( x , t ) =
0, where h (cid:48) ( x ) = dd x h ( x ) .Note that the gain function needs to be obtained for each value of time t . If the right handside of (9) is non-negative valued, it then follows from the minimum principle for elliptic BVPsthat the gain function K is non-negative valued [10].The contributions of this paper are as follows: • Variational Problem.
The construction of the feedback particle filter is based on a variationalproblem, where the cost function is the Kullback-Leibler (K-L) divergence between p ∗ ( x , t ) and p ( x , t ) . The feedback particle filter (6)-(9), including the formula (7) for the innovation error andthe E-L BVP (9), is obtained via analysis of the first variation. • Consistency.
The particle filter model (6) is consistent with nonlinear filter in the followingsense: Suppose the gain function K ( x , t ) is obtained as the solution to (9), and the priors areconsistent, p ( x , ) = p ∗ ( x , ) . Then, for all t ≥ x , p ( x , t ) = p ∗ ( x , t ) . • Algorithms.
Numerical techniques are proposed for synthesis of the gain function K ( x , t ) .If a ( · ) and h ( · ) are linear and the density p ∗ is Gaussian, then the gain function is simplythe Kalman gain. At time t , it is a constant given in terms of variance alone. The variance isapproximated empirically as a sample covariance. February 27, 2013 DRAFT
In the nonlinear case, numerical approximation techniques are described. Other approachesusing sum of Gaussian approximation also exist but are omitted on account of space. Details forthe latter can be found in [24].In recent decades, there have been many important advances in importance sampling basedapproaches for particle filtering; cf., [8], [3], [22]. A crucial distinction here is that there is noresampling of particles.We believe that the introduction of control in the feedback particle filter has several usefulfeatures/advantages:
Does not require sampling.
There is no re-sampling required as in the conventional particle filter.This property allows the feedback particle filter to be flexible with regards to implementationand does not suffer from sampling-related issues.
Innovation error.
The innovation error-based feedback structure is a key feature of the feedbackparticle filter (6). The innovation error in (6) is based on the average value of the prediction h ( X it ) of the i th -particle and the prediction ˆ h t due to the entire population.The feedback structure is easier to see when the filter is expressed in its Stratonovich form:d X it = a ( X it ) d t + d B it + K ( X i , t ) ◦ (cid:18) d Z t − ( h ( X it ) + ˆ h t ) d t (cid:19) . (10)Given that the Stratonovich form provides a mathematical interpretation of the (formal) ODEmodel [20, Section 3.3], we also obtain the ODE model of the filter. Denoting Y t . = d Z t d t andwhite noise process ˙ B it . = d B it d t , the ODE model of the filter is given by,d X it d t = a ( X it ) + ˙ B it + K ( X i , t ) · (cid:18) Y t − ( h ( X it ) + ˆ h t ) (cid:19) . The feedback particle filter thus provides for a generalization of the Kalman filter to nonlinearsystems, where the innovation error-based feedback structure of the control is preserved (seeFig. 1). For the linear case, the optimal gain function is the Kalman gain. For the nonlinear case,the Kalman gain is replaced by a nonlinear function of the state.
Feedback structure.
Feedback is important on account of the issue of robustness . A filter is basedon an idealized model of the underlying dynamic process that is often nonlinear, uncertain andtime-varying. The self-correcting property of the feedback provides robustness, allowing one totolerate a degree of uncertainty inherent in any model.In contrast, a conventional particle filter is based upon importance sampling. Although theinnovation error is central to the Kushner-Stratonovich’s stochastic partial differential equation(SPDE) of nonlinear filtering, it is conspicuous by its absence in a conventional particle filter.Arguably, the structural aspects of the Kalman filter have been as important as the algorithmitself in design, integration, testing and operation of the overall system. Without such structuralfeatures, it is a challenge to create scalable cost-effective solutions.
February 27, 2013 DRAFT
Fig. 1. Innovations error-based feedback structure for (a) Kalman filter and (b) nonlinear feedback particle filter.
The “innovation” of the feedback particle filter lies in the (modified) definition of innovationerror for a particle filter. Moreover, the feedback control structure that existed thusfar only forKalman filter now also exists for particle filters (compare parts (a) and (b) of Fig. 1).
Variance reduction.
Feedback can help reduce the high variance that is sometimes observed inthe conventional particle filter. Numerical results in Sec. V support this claim — See Fig. 3 fora comparison of the feedback particle filter and the bootstrap filter.
Ease of design, testing and operation.
On account of structural features, feedback particlefilter-based solutions are expected to be more robust, cost-effective, and easier to debug andimplement.
Applications.
Bayesian inference is an important paradigm used to model functions of certainneural circuits in the brain [9]. Compared to techniques that rely on importance sampling, afeedback particle filter may provide a more neuro-biologically plausible model to implementfiltering and inference functions [25]. This is illustrated here with the aid of a filtering probleminvolving nonlinear oscillators. Another application appears in [21].
A. Comparison with Relevant Literature
Our work is motivated by recent development in mean-field games, but the focus there hasbeen primarily on optimal control [14], [26].In nonlinear filtering, there are two directly related works: Crisan and Xiong [6], and Mitterand Newton [19]. In each of these papers, a controlled system is introduced, of the formd X it = (cid:0) a ( X it ) + u ( X it , t ) (cid:1) d t + σ B d B it . The objective is to choose the control input to obtain a solution of the nonlinear filtering problem.The approach in [19] is based on consideration of a finite-horizon optimal control problem.It leads to an HJB equation whose solution yields the optimal control input.The work of Crisan and Xiong is closer to our paper in terms of both goals and approaches.Although we were not aware of their work prior to submission of our original conferencepapers [25], [24], Crisan and Xiong provide an explicit expression for a control law that is similarto the feedback particle filter, with some important differences. One, the considerations of Crisan
February 27, 2013 DRAFT and Xiong (and also of Newton and Mitter) require introduction of a smooth approximation ofthe process “ d Z d t − ˆ h t ,” which we avoid with our formulation. Two, the filter derived in Crisanand Xiong has a structure based on a gain feedback with respect to the smooth approximation,while the feedback particle filter is based on the formula for innovation error I it as given in (7).This formula is fundamental to construction of particle filters in continuous time settings. Weclarify here that the formula for innovation error is not assumed, comes about as a result of theanalysis of the variational problem.Remarkably, both the feedback particle filter and Crisan and Xiong’s filter require solution ofthe same boundary value problem, and as such have the same computational complexity. TheBVP is solved to obtain the gain function. However, the particular solution described in Crisanand Xiong for the BVP may not work in all cases, including the linear Gaussian case. Additionaldiscussion appears in Sec. III-E.Apart from these two works, Daum and Huang have introduced the information flow filter forthe continuous-discrete time filtering problem [7]. Although an explicit formula for the filter isdifficult to obtain, a closely related form of the boundary value problem appears in their work.There is also an important discussion of both the limitations of the conventional particle filter,and the need to incorporate feedback to ameliorate these issues. Several numerical experimentsare presented that describe high variance and robustness issues, especially where signal modelsare unstable. These results provide significant motivation to the work described here. B. Outline
The variational setup is described in Sec. II: It begins with a discussion of the continuous-discrete filtering problem: the equation for dynamics is defined by (1a), but the observations aremade only at discrete times. The continuous-time filtering problem (for (1a)-(1b)) is obtained asa limiting case of the continuous-discrete problem.The feedback particle filter is introduced in Sec. III. Extension to the multivariable caseis briefly described in Sec. III-D, followed by a comparison with Crisan and Xiong’s filterin Sec. III-E.Algorithms are discussed in Sec. IV, and numerical examples are described in Sec. V, includingthe neuroscience application involving coupled oscillator models. These models (also consideredin our earlier mean-field control paper [26]) provided some of the initial motivation for thepresent work. II. V
ARIATIONAL P ROBLEM
The control problem posed by any one of the i th particles can be cast as a partially observedoptimal control problem. The observations are given by { X it , Z t } , and the state process is two-dimensional, { X it , X t } . In partially observed optimal control problems, it is typical to take the February 27, 2013 DRAFT “belief state” p ∗ t as the state process, which is known to serve as a sufficient statistic for optimalcontrol under general conditions. Since our cost function is taken as the KL divergence between p ∗ t and p t (defined in (2) and (5), respectively), a natural state process for the purposes of optimalcontrol is the triple { X it , p t , p ∗ t } .The precise formulation of the optimal control problem begins with the continuous time model,with sampled observations. The equation for dynamics is given by (1a), and the observationsare made only at discrete times { t n } : Y t n = h ( X t n ) + W (cid:52) t n , (11)where (cid:52) : = t n + − t n and { W (cid:52) t n } is i.i.d and drawn from N ( , σ W (cid:52) ) .The particle model in this case is a hybrid dynamical system: For t ∈ [ t n − , t n ) , the i th particleevolves according to the stochastic differential equation,d X it = a ( X it ) d t + σ B d B it , t n − ≤ t < t n , (12)where the initial condition X it n − is given. At time t = t n there is a potential jump that is determinedby the input U it n : X it n = X it − n + U it n , (13)where X it − n denotes the right limit of { X it : t n − ≤ t < t n } . The specification (13) defines the initialcondition for the process on the next interval [ t n , t n + ) .The filtering problem is to construct a control law that defines { U it n : n ≥ } such that p ( · , t n ) approximates p ∗ ( · , t n ) for each n ≥
1. To solve this problem we first define “belief maps” thatpropagate the conditional distributions of X and X i . A. Belief Maps
The observation history is denoted Y n : = σ { Y t i : i ≤ n , i ∈ N } . For each n , various conditionaldistributions are considered:1) p ∗ n and p ∗− n : The conditional distribution of X t n given Y n and Y n − , respectively.2) p n and p − n : The conditional distribution of X it n given Y n and Y n − , respectively.These densities evolve according to recursions of the form, p ∗ n = P ∗ ( p ∗ n − , Y t n ) , p n = P ( p n − , Y t n ) . (14)The mappings P ∗ and P can be decomposed into two parts. The first part is identical for eachof these mappings: the transformation that takes p n − to p − n coincides with the mapping from p ∗ n − to p ∗− n . In each case it is defined by the Kolmogorov forward equation associated with thediffusion on [ t n − , t n ) . February 27, 2013 DRAFT
The second part of the mapping is the transformation that takes p ∗− n to p ∗ n , which is obtainedfrom Bayes’ rule: Given the observation Y t n made at time t = t n , p ∗ n ( s ) = p ∗− n ( s ) · p Y | X ( Y t n | s ) p Y ( Y t n ) , s ∈ R , (15)where p Y denotes the pdf for Y t n , and p Y | X ( · | s ) denotes the conditional distribution of Y t n given X t n = s . Applying (11) gives, p Y | X ( Y t n | s ) = (cid:113) πσ W / (cid:52) exp (cid:18) − ( Y t n − h ( s )) σ W / (cid:52) (cid:19) . Combining (15) with the forward equation defines P ∗ .The transformation that takes p − n to p n depends upon the choice of control U it n in (13). Attime t = t n , we seek a control input U it n that is admissible . Definition 1 (Admissible Input):
The control sequence { U it n : n ≥ } is admissible if there isa sequence of maps { v n ( x ; y n ) } such that U it n = v n ( X it − n , Y t , . . . , Y t n ) for each n , and moreover,(i) E [ | U it n | ] < ∞ , and with probability one,lim x →± ∞ v n ( x , Y t , . . . , Y t n ) p − n ( x ) = . (ii) v n is twice continuously differentiable as a function of x .(iii) 1 + v (cid:48) n ( x ) is non-zero for all x , where v (cid:48) n ( x ) = dd x v n ( x ) .We will suppress the dependency of v n on the observations (and often the time-index n ),writing U it n = v ( x ) when X it − n = x . Under the assumption that 1 + v (cid:48) ( x ) is non-zero for all x , wecan write, p n ( x + ) = p − n ( x ) | + v (cid:48) ( x ) | , where x + = x + v ( x ) . (16) B. Variational Problem
Our goal is to choose an admissible input so that the mapping P approximates the mapping P ∗ in (14). More specifically, given the pdf p n − we have already defined the mapping P sothat p n = P ( p n − , Y t n ) . We denote ˆ p ∗ n = P ∗ ( p n − , Y t n ) , and choose v n so that these pdfs are asclose as possible. We approach this goal through the formulation of an optimization problemwith respect to the KL divergence metric. That is, at time t = t n , the function v n is the solutionto the following optimization problem, v n ( x ) = arg min v KL ( p n (cid:107) ˆ p ∗ n ) . (17)Based on the definitions, for any v the KL divergence can be expressed,KL ( p n (cid:107) ˆ p ∗ n ) = − (cid:90) R p − n ( x ) (cid:110) ln | + v (cid:48) ( x ) | + ln (cid:0) p − n ( x + v ( x )) p Y | X ( Y t n | x + v ( x )) (cid:1)(cid:111) d x + C , (18) February 27, 2013 DRAFT where C = (cid:82) R p − n ( x ) ln ( p − n ( x ) p Y ( Y t n )) d x is a constant that does not depend on v ; cf., App. VII-Afor the calculation.The solution to (17) is described in the following proposition, whose proof appears in App. VII-B. Proposition 2.1:
Suppose that the admissible input is obtained as the solution to the sequenceof optimization problems (17). Then for each n , the function v = v n is a solution of the followingEuler-Lagrange (E-L) BVP:dd x (cid:18) p − n ( x ) | + v (cid:48) ( x ) | (cid:19) = p − n ( x ) ∂∂ v (cid:0) ln ( p − n ( x + v ) p Y | X ( Y t n | x + v )) (cid:1) , (19)with boundary condition lim x →± ∞ v ( x ) p − n ( x ) = optimal control function . Additional details on the continuous-discrete time filter appear in our conference paper [25].III. F EEDBACK P ARTICLE F ILTER
We now consider the continuous time filtering problem (1a, 1b) introduced in Sec. I.
A. Belief State Dynamics and Control Architecture
The model for the particle filter is given by the Itˆo diffusion,d X it = a ( X it ) d t + σ B d B it + u ( X it , t ) d t + K ( X it , t ) d Z t (cid:124) (cid:123)(cid:122) (cid:125) d U it , (20)where X it ∈ R is the state for the i th particle at time t , and { B it } are mutually independent standardWiener processes. We assume the initial conditions { X i } Ni = are i.i.d., independent of { B it } , anddrawn from the initial distribution p ∗ ( x , ) of X . Both { B it } and { X i } are also assumed to beindependent of X t , Z t .As in Sec. II, we impose admissibility requirements on the control input U it in (20): Definition 2 (Admissible Input):
The control input U it is admissible if the random variables u ( x , t ) and K ( x , t ) are Z t = σ ( Z s : s ≤ t ) measurable for each t . Moreover, each t ,(i) E [ | u ( X it , t ) | + | K ( X it , t ) | ] < ∞ , and with probability one,lim x →± ∞ u ( x , t ) p ( x , t ) = , (21a)lim x →± ∞ K ( x , t ) p ( x , t ) = . (21b)where p is the posterior distribution of X it given Z t , defined in (5).(ii) u : R → R , K : R → R are twice continuously differentiable in their first arguments. February 27, 2013 DRAFT0
The functions { u ( x , t ) , K ( x , t ) } represent the continuous-time counterparts of the optimal con-trol function v n ( x ) (see (17)). We say that these functions are optimal if p ≡ p ∗ , where recall p ∗ is the posterior distribution of X t given Z t as defined in (2). Given p ∗ ( · , ) = p ( · , ) , ourgoal is to choose { u , K } in the feedback particle filter so that the evolution equations of theseconditional distributions coincide.The evolution of p ∗ ( x , t ) is described by the Kushner-Stratonovich (K-S) equation:d p ∗ = L † p ∗ d t + σ W ( h − ˆ h t )( d Z t − ˆ h t d t ) p ∗ , (22)where ˆ h t = (cid:82) h ( x ) p ∗ ( x , t ) d x , and L † p ∗ = − ∂ ( p ∗ a ) ∂ x + σ B ∂ p ∗ ∂ x .The evolution equation of p ( x , t ) is described next. The proof appears in App. VII-C. Proposition 3.1:
Consider the process X it that evolves according to the particle filter model (20).The conditional distribution of X it given the filtration Z t , p ( x , t ) , satisfies the forward equationd p = L † p d t − ∂∂ x ( K p ) d Z t − ∂∂ x ( up ) d t + σ W ∂ ∂ x (cid:0) p K (cid:1) d t . (23) B. Consistency with the Nonlinear Filter
The main result of this section is the construction of an optimal pair { u , K } under the followingassumption: Assumption A1
The conditional distributions ( p ∗ , p ) are C , with p ∗ ( x , t ) > p ( x , t ) > x ∈ R , t > { u , K } as the solution to a certain E-L BVP based on p : the function K is the solution to − ∂∂ x (cid:18) p ( x , t ) ∂∂ x { p ( x , t ) K ( x , t ) } (cid:19) = σ W h (cid:48) ( x ) , (24)with boundary condition (21b). The function u ( · , t ) : R → R is obtained as: u ( x , t ) = K ( x , t ) (cid:18) − ( h ( x ) + ˆ h t ) + σ W K (cid:48) ( x , t ) (cid:19) , (25)where ˆ h t = (cid:82) h ( x ) p ( x , t ) d x . We assume moreover that the control input obtained using { u , K } isadmissible. The particular form of u given in (25) and the BVP (24) is motivated by consideringthe continuous-time limit of (19), obtained on letting (cid:52) : = t n + − t n go to zero; the calculationsappear in App. VII-D.Existence and uniqueness of { u , K } is obtained in the following proposition — Its proof isgiven in App. VII-E. Proposition 3.2:
Consider the BVP (24), subject to Assumption A1. Then,
February 27, 2013 DRAFT1
1) There exists a unique solution K , subject to the boundary condition (21b).2) The solution satisfies K ( x , t ) ≥ x , t , provided h (cid:48) ( x ) ≥ x .The following theorem shows that the two evolution equations (22) and (23) are identical.The proof appears in App. VII-F. Theorem 3.3:
Consider the two evolution equations for p and p ∗ , defined according to thesolution of the forward equation (23) and the K-S equation (22), respectively. Suppose that thecontrol functions u ( x , t ) and K ( x , t ) are obtained according to (24) and (25), respectively. Then,provided p ( x , ) = p ∗ ( x , ) , we have for all t ≥ p ( x , t ) = p ∗ ( x , t ) Remark 1:
Thm. 3.3 is based on the ideal setting in which the gain K ( X it , t ) is obtained as afunction of the posterior p = p ∗ for X it . In practice the algorithm is applied with p replaced bythe empirical distribution of the N particles.In this ideal setting, the empirical distribution of the particle system will approximate theposterior distribution p ∗ ( x , t ) as N → ∞ . The convergence is in the weak sense in general. Toobtain almost sure convergence, it is necessary to obtain sample path representations of thesolution to the stochastic differential equation for each i (see e.g. [16]). Under these conditionsthe solution to the SDE (4) for each i has a functional representation, X it = F ( X i , B i [ , t ] ; Z [ , t ] ) , where the notation Z [ , t ] signifies the entire sample path { Z s : 0 ≤ s ≤ t } for a stochastic process Z ; F is a continuous functional (in the uniform topology) of the sample paths { B i [ , t ] , Z [ , t ] } along with the initial condition X i . It follows that the empirical distribution has a functionalrepresentation, p ( N ) ( A , t ) = N N ∑ i = { F ( X i , B i [ , t ] ; Z [ , t ] ) ∈ A } The sequence { ( X i , B i [ , t ] ) : i = , ... } is i.i.d. and independent of Z . It follows that the summand { { F ( X i , B i [ , t ] ; Z [ , t ] ) : i = , . . . } is also i.i.d. given Z [ , t ] . Almost sure convergence follows fromthe Law of Large Numbers for scalar i.i.d. sequences.In current research we are considering the more difficult problem of performance bounds forthe approximate implementations described in Sec. IV. Remark 2:
On integrating (24) once, we obtain an equivalent characterization of the E-L BVP: ∂∂ x ( p K ) = − σ W ( h − ˆ h t ) p , (26) February 27, 2013 DRAFT2 now with a single boundary condition lim x →− ∞ p K ( x , t ) =
0. The resulting gain function can bereadily shown to yield admissible control input under certain additional technical assumptionson density p and the function h .Given the scope of this paper, and the fact that the same apriori bounds apply also to themultivariable case, we defer additional discussion to a future publication. Remark 3:
Although the methodology and the filter is presented for Gaussian process andobservation noise, the case of non-Gaussian process noise is easily handled – simply replace thenoise model in the filter with the appropriate model of the process noise.For other types of observation noise, one would modify the conditional distribution p Y | X inthe optimization problem (17). The derivation of filter would then proceed by consideration ofthe first variation (see App. VII-D). C. Example: Linear Model
It is helpful to consider the feedback particle filter in the following simple linear setting,d X t = α X t d t + σ B d B t , (27a)d Z t = γ X t d t + σ W d W t , (27b)where α , γ are real numbers. The initial distribution p ∗ ( x , ) is assumed to be Gaussian withmean µ and variance Σ .The following lemma provides the solution of the gain function K ( x , t ) in the linear Gaussiancase. Lemma 3.4:
Consider the linear observation equation (27b). If p ( x , t ) is assumed to be Gaus-sian with mean µ t and variance Σ t , then the solution of E-L BVP (9) is given by: K ( x , t ) = Σ t γσ W . (28)The formula (28) is verified by direct substitution in the ODE (9) where the distribution p isGaussian.The optimal control yields the following form for the particle filter in this linear Gaussianmodel: d X it = α X it d t + σ B d B it + Σ t γσ W (cid:18) d Z t − γ X it + µ t t (cid:19) . (29)Now we show that p = p ∗ in this case. That is, the conditional distributions of X and X i coincide, and are defined by the well-known dynamic equations that characterize the mean andthe variance of the continuous-time Kalman filter. The proof appears in App. VII-G. February 27, 2013 DRAFT3
Theorem 3.5:
Consider the linear Gaussian filtering problem defined by the state-observationequations (27a,27b). In this case the posterior distributions of X and X i are Gaussian, whoseconditional mean and covariance are given by the respective SDE and the ODE,d µ t = α µ t d t + Σ t γσ W (cid:16) d Z t − γ µ t d t (cid:17) (30)dd t Σ t = α Σ t + σ B − ( γ ) Σ t σ W (31)Notice that the particle system (29) is not practical since it requires computation of theconditional mean and variance { µ t , Σ t } . If we are to compute these quantities, then there isno reason to run a particle filter!In practice { µ t , Σ t } are approximated as sample means and sample covariances from theensemble { X it } Ni = : µ t ≈ µ ( N ) t : = N N ∑ i = X it , Σ t ≈ Σ ( N ) t : = N − N ∑ i = ( X it − µ ( N ) t ) . (32)The resulting equation (29) for the i th particle is given byd X it = α X it d t + σ B d B it + Σ ( N ) t γσ W (cid:32) d Z t − γ X it + µ ( N ) t t (cid:33) . (33)It is very similar to the mean-field “synchronization-type” control laws and oblivious equilibriaconstructions as in [14], [26]. The model (29) represents the mean-field approximation obtainedby letting N → ∞ . D. Feedback Particle Filter for the Multivariable Model
Consider the model (1a)-(1b) in which the state X t is d -dimensional, with d ≥
2, so that a ( · ) is a vector-field on R d . For ease of presentation σ B is assumed to be scalar, and the observationprocess Z t ∈ R real-valued.To aid comparison with Crisan and Xiong’s work, we express the feedback particle filter inits Stratonovich form: d X it = a ( X it ) d t + σ B d B it + K ( X it , t ) ◦ d I it (34)where the innovation error is as before,d I it : = d Z t − ( h ( X it ) + ˆ h t ) d t , (35) February 27, 2013 DRAFT4 and the gain function K ( x , t ) = ( K , K , ..., K d ) T is now vector-valued. It is given by the solutionof a BVP, the multivariable counterpart of (26): ∇ · ( p K ) = − σ W ( h − ˆ h t ) p , (36)where ∇ · denotes the divergence operator.It is straightforward to prove consistency by repeating the steps in the proof of Thm. 3.3, nowwith the Kolmogorov forward operator:d p = L † p d t − ∇ · ( K p ) d Z t − ∇ · ( up ) d t + σ W d ∑ i , j = ∂ [( KK T ) i j p ] ∂ x i ∂ x j d t (37)where L † p = − ∇ · ( pa ) + σ B ∆ p , ∆ is the Laplacian, and u is the multivariable counterpartof (25): u = − K ( x , t ) h ( x ) + ˆ h t + Ω ( x , t ) , where Ω = ( Ω , Ω , ..., Ω d ) T is the Wong-Zakai correction term: Ω l ( x , t ) : = σ W d ∑ k = K k ( x , t ) ∂ K l ∂ x k ( x , t ) . As with the scalar case, the multivariable feedback particle filter requires solution of aBVP (36) at each time step.Following the work of Crisan and Xiong [6], we might assume the following representationin an attempt to solve (36), p K = ∇ φ . (38)where φ is assumed to be sufficiently smooth. Substituting (38) in (36) yields the Poissonequation, ∆ φ = − σ W ( h − ˆ h t ) p . (39)A solution to Poisson’s equation with d ≥ G ( r ) = π ln ( r ) for d = d ( − d ) ω d r − d for d > ω d is the volume of the unit ball in R d . A solution to (39) is then given by, φ ( x ) = − σ W (cid:90) R d G ( | y − x | )( h ( y ) − ˆ h t ) p ( y , t ) d y , where | y − x | : = (cid:16) ∑ dj = ( y j − x j ) (cid:17) is the Euclidean distance. February 27, 2013 DRAFT5
On taking the gradient and using (38), one obtains an explicit formula for the gain function: K ( x , t ) = σ W Γ ( p , h )( x ) p ( x , t ) = : K g ( x , t ) , (40)where Γ ( p , g )( x ) : = d ω d (cid:82) y − x | y − x | d ( g ( y ) − ˆ g ) p ( y ) d y .While this leads to a solution to (36), it may not lead to an admissible control law. Thisdifficulty arises in the prior work of Crisan and Xiong. E. Comparison with Crisan and Xiong’s Filter
In [6] and in Sec. 4 of [22], a particle filter of the following form is presented:d X it = a ( X it ) d t + σ B d B it + K g ( X it , t ) (cid:16) dd t ˜ I t (cid:17) d t , (41)where ˜ I t is a certain smooth approximation obtained from the standard form of the innovationerror I t : = Z t − (cid:82) t ˆ h t d t . A consistency result is described for this filter.We make the following comparisons:1) Without taking a smooth approximation, the filter (41) is formally equivalent to the followingSDE expressed here in its Stratonovich form:d X it = a ( X it ) d t + σ B d B it + K g ( X it , t ) ◦ (cid:16) d Z t − ˆ h t d t (cid:17) . (42)In this case, using (37), it is straightforward to show that the consistency result does not hold .In particular, there is an extra second order term that is not present in the K-S equation forevolution of the true posterior p ∗ ( x , t ) .2) The feedback particle filter introduced in this paper does not require a smooth approximationand yet achieves consistency. The key breakthrough is the modified definition of the innovationerror (compare (42) with (34)). Note that the innovation error (35) is not assumed apriori butcomes about via analysis of the variational problem. This is one utility of introducing thevariational formulation. Once the feedback particle filter has been derived, it is straightforwardto prove consistency (see the Proof of Thm. 3.3).3) The computational overhead for the feedback particle filter and the filter of Crisan and Xiongare equal. Both require the approximation of the integral (40), and division by (a suitableregularized approximation of) p ( x , t ) . Numerically, the Poisson equation formulation (39) ofthe E-L BVP (36) is convenient. There exist efficient numerical algorithms to approximate theintegral solution (40) for a system of N particles in arbitrary dimension; cf., [12].However, while appealing, the function K g is not the correct form of the gain function in themultivariable case, even for linear models: It is straightforward to verify that the Kalman gainis a solution of the boundary value problem (36). Using the Kalman gain for the gain function February 27, 2013 DRAFT6 in (34) yields the feedback particle filter for the multivariable linear Gaussian case. The filter isa straightforward extension of (33) in Sec. III-C.However, the Kalman gain solution is not of the form (38). Thus, the integral solution (40)does not equal the Kalman gain in the linear Gaussian case (for d ≥ | K g ( x , t ) | → ∞ as | x | → ∞ , and E [ | K g | ] = E [ | K g | ] = ∞ . A proof is given in App. VII-H.It follows that the control input obtained using K g is not admissible, and hence the Kolmogorovforward operator is no longer valid. Filter implementations using K g suffer from numerical issueson account of large unbounded gains. In contrast, the feedback particle filter using Kalman gainworks both in theory and in practice.The choice of gain function in the multivariable case requires careful consideration of theuniqueness of the solutions of the BVP: The solution of (36) is not unique, even thoughuniqueness holds when p K is assumed to be of a gradient-form.Before closing this section, we note that [6, Proposition 2.4] concerns another filter that doesnot rely on smooth approximation,d X it = a ( X it ) d t + σ B d B it + K g ( X it , t ) ◦ d Z t − σ W Γ ( p , | h | )( X it ) p ( X it , t ) d t . Our calculations indicate that consistency is also an issue for this filter. The issue with K g alsoapplies to this filter. A more complete comparison needs further investigation.IV. S YNTHESIS OF THE G AIN F UNCTION
Implementation of the nonlinear filter (6) requires solution of the E-L BVP (9) to obtain thegain function K ( x , t ) for each fixed t . A. Direct Numerical Approximation of the BVP solution
The explicit closed-form formula (76) for the solution of the BVP (9) can be used to constructa direct numerical approximation of the solution. Using (76), we have K ( x , t ) = p ( x , t ) σ W (cid:90) x − ∞ ( ˆ h t − h ( y )) p ( y , t ) d y . The approximation involves three steps:1) Approximation of ˆ h t by using a sample mean:ˆ h t ≈ N N ∑ j = h ( X jt ) = : ˆ h ( N ) t .
2) Approximation of the integrand: ( ˆ h t − h ( y )) p ( y , t ) ≈ N N ∑ j = ( ˆ h ( N ) t − h ( X jt )) δ ( y − X jt ) , February 27, 2013 DRAFT7 where δ ( · ) is the Dirac delta function.3) Approximation of the density p ( x , t ) in the denominator, e.g., as a sum of Gaussian: p ( x , t ) ≈ N N ∑ j = q jt ( x ) = : ˜ p ( x , t ) , (43)where q jt ( x ) = q ( x ; X jt , (cid:15) ) = √ π (cid:15) exp (cid:16) − (cid:15) ( x − X jt ) (cid:17) . The appropriate value of (cid:15) depends uponthe problem. As a function of N , (cid:15) can be made smaller as N grows; As N → ∞ , (cid:15) → K ( x , t ) = p ( x , t ) σ W N N ∑ j = ( ˆ h ( N ) t − h ( X jt )) H ( x − X jt ) , (44)where H ( · ) is the Heaviside function.Note that the gain function needs to be evaluated only at the particle locations X it . An efficient O ( N ) algorithm is easily constructed to do the same: K ( X it , t ) = p ( X it , t ) σ W N ∑ j : X jt < X it ( ˆ h ( N ) t − h ( X jt )) + ( ˆ h ( N ) t − h ( X it )) K (cid:48) ( X it , t ) = σ W ( ˆ h ( N ) t − h ( X it )) − ˜ b ( X it ) K ( X it , t ) , (45)where ˜ b ( x ) : = ∂∂ x ( ln ˜ p )( x , t ) . For ˜ p defined using the sum of Gaussian approximation (43), aclosed-form formula for ˜ b ( x ) is easily obtained. B. Algorithm
For implementation purposes, we use the Stratonovich form of the filter (see (10)) together withan Euler discretization. The resulting discrete-time algorithm appears in Algorithm 1. At eachtime step, the algorithm requires approximation of the gain function. A DNS-based algorithmfor the same is summarized in Algorithm 2.In practice, one can use a less computationally intensive algorithm to approximate the gainfunction. An algorithm based on sum-of-Gaussian approximation of density appears in ourconference paper [24]. In the application example presented in Sec. V-C, the gain functionis approximated by using Fourier series.
C. Further Remarks on the BVP
Recall that the solution of the nonlinear filtering problem is given by the Kushner-Stratonovichnonlinear evolution PDE. The feedback particle filter instead requires, at each time t , a solutionof the linear BVP (36) to obtain the gain function K : ∇ · ( p K ) = − σ W ( h − ˆ h t ) p . February 27, 2013 DRAFT8
Algorithm 1
Implementation of feedback particle filter Initialization for i : = N do Sample X i from p ( x , ) end for Assign value t : = Iteration [from t to t + (cid:52) t ] Calculate ˆ h ( N ) t : = N ∑ Ni = h ( X it ) for i : = N do Generate a sample, (cid:52) V , from N ( , ) Calculate (cid:52) I it : = (cid:52) Z t − (cid:16) h ( X it ) + ˆ h ( N ) t (cid:17) (cid:52) t Calculate the gain function K ( X it , t ) (e.g., by using Alg. 2) X it + (cid:52) t : = X it + a ( X it ) (cid:52) t + σ B √ (cid:52) t (cid:52) V + K ( X it , t ) (cid:52) I it end for t : = t + (cid:52) t Algorithm 2
Synthesis of gain function K ( x , t ) Calculate ˆ h t ≈ ˆ h ( N ) t ; Approximate p ( x , t ) as a sum of Gaussian: p ( x , t ) ≈ ˜ p ( x , t ) : = N N ∑ j = q jt ( x ) , where q jt ( x ) = √ π (cid:15) exp (cid:18) − ( x − X jt ) (cid:15) (cid:19) . Calculate the gain function K ( x , t ) : = p ( x , t ) σ W N N ∑ j = (cid:16) ˆ h ( N ) t − h ( X jt ) (cid:17) H ( x − X jt ) , where H ( · ) is the Heaviside function.We make the following remarks:1) There are close parallels between the proposed algorithm and the vortex element method(VEM) developed by Chorin and others for solution of the Navier-Stokes evolution PDE;cf., [4], [18]. In VEM, as in the feedback particle filter, one obtains the solution of a nonlinearevolution PDE by flowing a large number of particles. The vector-field for the particles isobtained by solving a linear BVP at each time. February 27, 2013 DRAFT9 N = 10 t N = 10 N = 10 X ( t )¯ µ N ( t ) StateCond.mean M.C. est. Cond. cov. Cond. cov. M.C. est. Σ( t )¯Σ N ( t ) Cond. cov. Cond. cov. M.C. est. Σ( t )¯Σ N ( t ) Fig. 2. (a) Comparison of the true state { X t } and the conditional mean { ¯ µ ( N ) t } . (b) and (c) Plots of estimated conditionalcovariance with N = ,
000 and N =
100 particles, respectively. For comparison, the true conditional covariance obtained usingKalman filtering equations is also shown.
Algorithms based on VEM are popular in the large Reynolds number regime when the domainis not too complicated. The latter requirement is necessary to obtain solution of the linear BVPin tractable fashion [12].2) One may ask what is the benefit, in terms of accuracy and computational cost, of the feedbackparticle filter-based solution when compared to a direct solution of the nonlinear PDE (Kushner-Stratonovich equation) or the linear PDE (Zakai equation)?The key point, we believe, is robustness on account of the feedback control structure. Specifi-cally, the self-correcting property of the feedback provides robustness, allowing one to toleratea degree of uncertainty inherent in any model or approximation scheme. This is expected toyield accurate solutions in a computationally efficient manner. A complete answer will requirefurther analysis, and as such reflects an important future direction.3) The biggest computational cost of our approach is the need to solve the BVP at each time-step, that additionally requires one to approximate the density. We are encouraged howeverby the extensive set of tools in feedback control: after all, one rarely needs to solve the HJBequations in closed-form to obtain a reasonable feedback control law. Moreover, there aremany approaches in nonlinear and adaptive control to both approximate control laws as wellas learn/adapt these in online fashion; cf., [2].V. N
UMERICS
A. Linear Gaussian Case
Consider the linear system:d X t = α X t d t + d B t , (46a)d Z t = γ X t d t + σ W d W t , X ∼ N ( , ) , (46b) February 27, 2013 DRAFT0 where { B t } , { W t } are mutually independent standard Wiener process, and parameters α = − . γ = σ W = . N particles is described by the linear SDE,d X it = α X it d t + d B it + γ ¯ Σ ( N ) t σ W [ d Z t − γ X it + ¯ µ ( N ) t t ] , (47)where { B it } are mutually independent standard Wiener process; the particle system is initializedby drawing initial conditions { X i } Ni = from the distribution N ( , ) , and the parameter values arechosen according to the model.In the simulation discussed next, the mean ¯ µ ( N ) t and the variance ¯ Σ ( N ) t are obtained from theensemble { X it } Ni = according to (32).Fig. 2 summarizes some of the results of the numerical experiments: Part (a) depicts a samplepath of the state { X t } and the mean { ¯ µ ( N ) t } obtained using a particle filter with N = , Σ ( N ) t and the true errorvariance Σ t that one would obtain by using the Kalman filtering equations. The accuracy of theresults is sensitive to the number of particles. For example, part (c) of the figure provides acomparison of the variance with N =
100 particles.
Comparison with the bootstrap filter : We next provide a performance comparison betweenthe feedback particle filter and the bootstrap particle filter for the linear problem (46a, 46b) inregard to both error and running time.For the linear filtering problem, the optimal solution is given by the Kalman filter. We usethis solution to define the relative mean-squared error: mse = T (cid:90) T (cid:32) Σ ( N ) t − Σ t Σ t (cid:33) d t , (48)where Σ t is the error covariance using the Kalman filter, and Σ ( N ) t is its approximation using theparticle filter.Fig. 3(a) depicts a comparison between mse obtained using the feedback particle filter (47) andthe bootstrap filter. The latter implementation is based on an algorithm taken from Ch. 9 of [1].For simulation purposes, we used a range of values of α ∈ {− . , , . } , γ = σ B = σ W = . (cid:52) t = .
01, and T =
50. The plot is generated using simulations with N = , , , , , α > α >
0, our numerical simulations with the bootstrap filter “blew-up” (similar conclusions
February 27, 2013 DRAFT1 (b) Time
BPFFPF (a) Error
FPF −0.500.5BPF −0.50
Bootstrap (BPF)Feedback (FPF)
Fig. 3. Comparison of (a) the mse , (b) the computational time using feedback particle filter and the bootstrap particle filter. were also arrived independently in [7]) while the feedback particle filter is provably stable basedon observability of the model (46a)-(46b) (see Fig. 3 mse plot with α = . p ( x , t ) . Detailed comparisons between the feedbackparticle filter and the bootstrap particle filter will appear elsewhere.In general, the main computational burden of the feedback particle filter is to obtain gainfunction which can be made efficient by using various approximation approaches. B. Nonlinear example
This nonlinear SDE is chosen to illustrate the tracking capability of the filter in highly nonlinearsettings, d X t = X t ( − X t ) d t + σ B d B t , (49a)d Z t = X t d t + σ W d W t . (49b) February 27, 2013 DRAFT2 X ( t ) E ( t ) X ( t ) Fig. 4. Comparison of the true state X ( t ) and the conditional mean ¯ X ( t ) by using feedback particle filter. The error E ( t ) = X ( t ) − ¯ X ( t ) remains small even during a transition of the state. When σ B =
0, the ODE (49a) has two stable equilibria at ±
1. With σ B >
0, the state of the SDE“transitions” between these two “equilibria”.Fig. 4 depicts the simulation results obtained using the nonlinear feedback particle filter (6),with σ B = . σ W = .
2. The implementation is based on an algorithm described in Sec. IVof [24], and the details are omitted here on account of space. We initialize the simulation withtwo Gaussian clusters. After a brief period of transients, these clusters merge into a single cluster,which adequately tracks the true state including the transition events.
C. Application: nonlinear oscillators
We consider the filtering problem for a nonlinear oscillator:d θ t = ω d t + σ B d B t mod 2 π , (50)d Z t = h ( θ t ) d t + σ W d W t , (51)where ω is the frequency, h ( θ ) = [ + cos ( θ )] , and { B t } and { W t } are mutually independentstandard Wiener process. For numerical simulations, we pick ω = σ B = . σ W = .
4. We consider oscillator models because of their significanceto applications including neuroscience; cf., [26].The feedback particle filter is given by:d θ it = ω d t + σ B d B it + K ( θ it , t ) ◦ [ d Z t − ( h ( θ it ) + ˆ h t ) d t ] mod 2 π , (52) i = , ..., N , where the function K ( θ , t ) is obtained via the solution of the E-L equation: − ∂∂ θ (cid:18) p ( θ , t ) ∂∂ θ { p ( θ , t ) K ( θ , t ) } (cid:19) = − sin θ σ W . (53) February 27, 2013 DRAFT3
Although the equation (53) can be solved numerically to obtain the optimal control function K ( θ , t ) , here we investigate a solution based on perturbation method. Suppose, at some time t , p ( θ , t ) = π = : p , the uniform density. In this case, the E-L equation is given by: ∂ θ θ K = sin θ σ W . A straightforward calculation shows that the solution in this case is given by K ( θ , t ) = − sin θ σ W = : K ( θ ) . (54)To obtain the solution of the E-L equation (53), we assume that the density p ( θ , t ) is a smallharmonic perturbation of the uniform density. In particular, we express p ( θ , t ) as: p ( θ , t ) = p + ε ˜ p ( θ , t ) , (55)where ε is a small perturbation parameter. Since p ( θ , t ) is a density, (cid:82) π ˜ p ( θ , t ) d θ = K ( θ , t ) = K ( θ ) + ε ˜ K ( θ , t ) . (56)On substituting the ansatz (55) and (56) in (53), and retaining only O ( ε ) term, we obtain thefollowing linearized equation: ∂ θ θ ˜ K = − π∂ θ [( ∂ θ ˜ p ) K ] . (57)The linearized E-L equation (57) can be solved easily by considering a Fourier series expansionof ε ˜ p ( θ , t ) : ε ˜ p ( θ , t ) = P c ( t ) cos θ + P s ( t ) sin θ + h.o.h , (58)where “h.o.h” denotes the terms due to higher order harmonics. The Fourier coefficients aregiven by, P c ( t ) = π (cid:90) π p ( θ , t ) cos θ d θ , P s ( t ) = π (cid:90) π p ( θ , t ) sin θ d θ . For a harmonic perturbation, the solution of the linearized E-L equation (57) is given by: ε ˜ K ( θ , t ) = π σ W ( P c ( t ) sin 2 θ − P s ( t ) cos 2 θ ) = : K ( θ ; P c ( t ) , P s ( t )) (59)For “h.o.h” terms in the Fourier series expansion (58) of the density in p ( θ , t ) , the linearizedE-L equation (57) can be solved in a similar manner. In numerical simulation provided here,we ignore the higher order harmonics, and use a control input as summarized in the followingproposition: Proposition 5.1:
Consider the E-L equation (53) where the density p ( θ , t ) is assumed to be asmall harmonic perturbation of the uniform density π , as defined by (55) and (58). As ε → K ( θ , t ) = K ( θ ) + K ( θ ; P c ( t ) , P s ( t )) + o ( ε ) , (60) February 27, 2013 DRAFT4 t t θ π π π π (a) (b) (c) BoundsStateConditional mean est. θ ( t )¯ θ N ( t ) Second harmonicFirst harmonic PDF estimate x -2 x -1 ˜ p ( θ, t ) Fig. 5. Summary of the numerical experiments with the nonlinear oscillator filter: (a) Comparison of the true state { θ t } andthe conditional mean { ¯ θ Nt } . (b) The mean-squared estimate of the first and second harmonics of the density p ( θ , t ) and (c) aplot of a typical empirical distribution. where P c ( t ) , P s ( t ) denote the harmonic coefficients of density p ( θ , t ) . For large N , these areapproximated by using the formulae: P c ( t ) ≈ π N N ∑ j = cos θ j ( t ) , P s ( t ) ≈ π N N ∑ j = sin θ j ( t ) . (61)We next discuss the result of numerical experiments. The particle filter model is given by (52)with gain function K ( θ it , t ) , obtained using formula (60). The number of particles N = , { θ i } Ni = was sampled from a uniform distribution on circle [ , π ] .Fig. 5 summarizes some of the results of the numerical simulation. For illustration purposes,we depict only a single cycle from a time-window after transients due to initial condition haveconverged. Part (a) of the figure compares the sample path of the actual state { θ t } (as a dashedline) with the estimated mean { ¯ θ ( N ) t } (as a solid line). The shaded area indicates ± one standarddeviation bounds. Part (b) of the figure provides a comparison of the magnitude of the first andthe second harmonics (as dashed and solid lines, respectively) of the density p ( θ , t ) . The densityat any time instant during the time-window is approximately harmonic (see also part (c) wherethe density at one typical time instant is shown).Note that at each time instant t , the estimated mean, the bounds and the density p ( θ , t ) shownhere are all approximated from the ensemble { θ it } Ni = . For the sake of illustration, we have useda Gaussian mixture approximation to construct a smooth approximation of the density.VI. C ONCLUSIONS
In this paper, we introduced a new formulation of the nonlinear filter, referred to as the feedback particle filter . The feedback particle filter provides for a generalization of the Kalman
February 27, 2013 DRAFT5 filter to a general class of nonlinear non-Gaussian problems. Feedback particle filter inheritsmany of the properties that has made the Kalman filter so widely applicable over the past fivedecades, including innovation error and the feedback structure (see Fig. 1).Feedback is important on account of the issue of robustness . In particular, feedback can helpreduce the high variance that is sometimes observed in the conventional particle filter. Numericalresults are presented to support this claim (see Fig. 3).Even more significantly, the structural aspects of the Kalman filter have been as important asthe algorithm itself in design, integration, testing and operation of a larger system involvingfiltering problems (e.g., navigation systems). We expect feedback particle filter to similarlyprovide for an integrated framework, now for nonlinear non-Gaussian problems. We refer thereader to our paper [23] where feedback particle filter-based algorithms for nonlinear filteringwith data association uncertainty are described.VII. A
PPENDIX
A. Calculation of KL divergence
Recall the definition of K-L divergence for densities,KL ( p n (cid:107) ˆ p ∗ n ) = (cid:90) R p n ( s ) ln (cid:16) p n ( s ) ˆ p ∗ n ( s ) (cid:17) d s . We make a co-ordinate transformation s = x + v ( x ) and use (16) to express the K-L divergenceas: KL ( p n (cid:107) ˆ p ∗ n ) = (cid:90) R p − n ( x ) | + v (cid:48) ( x ) | ln ( p − n ( x ) | + v (cid:48) ( x ) | ˆ p ∗ n ( x + v ( x )) ) | + v (cid:48) ( x ) | d x The expression for K-L divergence given in (18) follows on using (15).
B. Solution of the optimization problem
Denote: L ( x , v , v (cid:48) ) = − p − n ( x ) (cid:16) ln | + v (cid:48) | + ln ( p − n ( x + v ) p Y | X ( Y t n | x + v )) (cid:17) . (62)The optimization problem (17) is a calculus of variation problem:min v (cid:90) L ( x , v , v (cid:48) ) d x . The minimizer is obtained via the analysis of first variation given by the well-known Euler-Lagrange equation: ∂ L ∂ v = dd x (cid:18) ∂ L ∂ v (cid:48) (cid:19) , Explicitly substituting the expression (62) for L , we obtain (19). February 27, 2013 DRAFT6
C. Derivation of the Forward Equation
We denote the filtration B t = σ ( X i , B is : s ≤ t ) , and we recall that Z t = σ ( Z s : s ≤ t ) for t ≥ independent by construction.On denoting ˜ a ( x , t ) = a ( x ) + u ( x , t ) , the particle evolution (20) is expressed, X it = X i + (cid:90) t ˜ a ( X is , s ) d s + (cid:90) t K ( X is , s ) d Z ( s ) + σ B B it . (63)By assumption on Lipschitz continuity of ˜ a and K , there exists a unique solution that is adaptedto the larger filtration B t ∨ Z t = σ ( X i , B is , Z s : s ≤ t ) . In fact, there is a functional F t such that, X it = F t ( X i , B it , Z t ) , (64)where Z t : = { Z s : 0 ≤ s ≤ t } denotes the trajectory.The conditional distribution of X it given Z t = σ ( Z s : s ≤ t ) was introduced in Sec. II-A: Itsdensity is denoted p ( x , t ) , defined by any bounded and measurable function f : R → R via, E [ f ( X it ) | Z t ] = (cid:90) R p ( x , t ) f ( x ) d x = : (cid:104) p t , f (cid:105) . We begin with a result that is the key to proving Prop. 3.1. The proof of Lemma 7.1 is omittedon account of space.
Lemma 7.1:
Suppose that f is an B t ∨ Z t -adapted process satisfying E (cid:82) t | f ( s ) | d s < ∞ . Then, E (cid:104) (cid:90) t f ( s ) d s | Z t (cid:105) = (cid:90) t E [ f ( s ) | Z s ] d s , (65) E (cid:104) (cid:90) t f ( s ) d Z s | Z t (cid:105) = (cid:90) t E [ f ( s ) | Z s ] d Z s . (66)We now provide a proof of the Proposition 3.1. Proof of Proposition 3.1
Applying Itˆo’s formula to equation (20) gives, for any smooth andbounded function f ,d f ( X it ) = L f ( X it ) d t + K ( X it , t ) ∂ f ∂ x ( X it ) d Z t + σ B ∂ f ∂ x ( X it ) d B it , where L f : = ( a + u ) ∂ f ∂ x + ( σ W K + σ B ) ∂ f ∂ x . Therefore, f ( X it ) = f ( X i ) + (cid:90) t L f ( X is ) d s + (cid:90) t K ( X is , s ) ∂ f ∂ x ( X is ) d Z s + σ B (cid:90) t ∂ f ∂ x ( X is ) d B is . Taking conditional expectations on both sides, (cid:104) p t , f (cid:105) = E [ f ( X i ) | Z t ] + E (cid:104) (cid:90) t L f ( X is ) d s | Z t (cid:105) + E (cid:104) (cid:90) t K ( X is , s ) ∂ f ∂ x ( X is ) d Z s | Z t (cid:105) + σ B E (cid:104) (cid:90) t ∂ f ∂ x ( X is ) d B is | Z t (cid:105) On applying Lemma 7.1, and the fact that B it is a Wiener process, we conclude that (cid:104) p t , f (cid:105) = (cid:104) p , f (cid:105) + (cid:90) t (cid:104) p s , L f (cid:105) d s + (cid:90) t (cid:104) p s , K ∂ f ∂ x (cid:105) d Z s . The forward equation (23) follows using integration by parts.
February 27, 2013 DRAFT7
D. Euler-Lagrange equation for the continuous-time filter
In this section we describe, formally, the continuous-time limit of the discrete-time E-LBVP (19). In the continuous-time case, the control and the observation models are of the form(see (20) and (1b)): d U it = u ( X it , t ) d t + K ( X it , t ) d Z t , d Z t = h ( X t ) d t + σ W d W t . In discrete-time, these are approximated as (cid:52) U it = u ( X it , t ) (cid:52) t + K ( X it , t ) (cid:52) Z t , (67) (cid:52) Z t = h ( X t ) (cid:52) t + σ W (cid:52) W t , where (cid:52) t is the small time-increment at t . It follows that the conditional distribution of Y t . = (cid:52) Z t (cid:52) t given X t is the density, p Y | X ( Y t |· ) = (cid:113) πσ W / (cid:52) t exp (cid:18) − ( (cid:52) Z t − h ( · ) (cid:52) t ) σ W (cid:52) t (cid:19) . (68)Substituting (67)-(68) in the E-L BVP (19) for the continuous-discrete time case, we arriveat the formal equation: ∂∂ x (cid:18) p ( x , t ) + u (cid:48) (cid:52) t + K (cid:48) (cid:52) Z t (cid:19) = p ( x , t ) ∂∂ v (cid:16) ln p ( x + v , t ) + ln p Y | X ( Y t | x + v ) (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) v = u (cid:52) t + K (cid:52) Z t . (69)For notational ease, we use primes to denote partial derivatives with respect to x : p is usedto denote p ( x , t ) , p (cid:48) : = ∂ p ∂ x ( x , t ) , p (cid:48)(cid:48) : = ∂ p ∂ x ( x , t ) , u (cid:48) : = ∂ u ∂ x ( x , t ) , K (cid:48) : = ∂ K ∂ x ( x , t ) etc. Note that thetime t is fixed.A sketch of calculations to obtain (24) and (25) starting from (69) appears in the followingthree steps: Step 1:
The three terms in (69) are simplified as: ∂∂ x (cid:18) p + u (cid:48) (cid:52) t + K (cid:48) (cid:52) Z t (cid:19) = p (cid:48) − f (cid:52) t − ( p (cid:48) K (cid:48) + p K (cid:48)(cid:48) ) (cid:52) Z t p ∂∂ v ln p ( x + v ) (cid:12)(cid:12)(cid:12)(cid:12) v = u (cid:52) t + K (cid:52) Z t = p (cid:48) + f (cid:52) t + ( p (cid:48)(cid:48) K − p (cid:48) K p ) (cid:52) Z t p ∂∂ v ln p Y | X ( Y t | x + v ) (cid:12)(cid:12)(cid:12)(cid:12) v = u (cid:52) t + K (cid:52) Z t = p σ W ( h (cid:48) (cid:52) Z t − hh (cid:48) (cid:52) t ) + ph (cid:48)(cid:48) K (cid:52) t where we have used Itˆo’s rules ( (cid:52) Z t ) = σ W (cid:52) t , (cid:52) Z t (cid:52) t = f = ( p (cid:48) u (cid:48) + pu (cid:48)(cid:48) ) − σ W ( p (cid:48) K (cid:48) + p K (cid:48) K (cid:48)(cid:48) ) , f = ( p (cid:48)(cid:48) u − p (cid:48) up ) + σ W K (cid:18) p (cid:48)(cid:48)(cid:48) − p (cid:48) p (cid:48)(cid:48) p + p (cid:48) p (cid:19) . February 27, 2013 DRAFT8
Collecting terms in O ( (cid:52) Z t ) and O ( (cid:52) t ) , after some simplification, leads to the following ODEs: E ( K ) = σ W h (cid:48) ( x ) (70) E ( u ) = − σ W h ( x ) h (cid:48) ( x ) + h (cid:48)(cid:48) ( x ) K + σ W G ( x , t ) (71)where E ( K ) = − ∂∂ x (cid:16) p ( x , t ) ∂∂ x { p ( x , t ) K ( x , t ) } (cid:17) , and G = − K (cid:48) K (cid:48)(cid:48) − ( K (cid:48) ) ( ln p ) (cid:48) + K ( ln p ) (cid:48)(cid:48)(cid:48) . Step 2.
Suppose ( u , K ) are admissible solutions of the E-L BVP (70)-(71). Then it is claimedthat − ( p K ) (cid:48) = h − ˆ h σ W p (72) − ( pu ) (cid:48) = − ( h − ˆ h ) ˆ h σ W p − σ W ( p K ) (cid:48)(cid:48) . (73)Recall that admissible here meanslim x →± ∞ p ( x , t ) u ( x , t ) = , lim x →± ∞ p ( x , t ) K ( x , t ) = . (74)To show (72), integrate (70) once to obtain − ( p K ) (cid:48) = σ W hp + C p , where the constant of integration C = − ˆ h σ W is obtained by integrating once again between − ∞ to ∞ and using the boundary conditions for K (74). This gives (72).To show (73), we denote its right hand side as R and claim (cid:18) R p (cid:19) (cid:48) = − hh (cid:48) σ W + h (cid:48)(cid:48) K + σ W G . (75)The equation (73) then follows by using the ODE (71) together with the boundary conditions for u (74). The verification of the claim involves a straightforward calculation, where we use (70)to obtain expressions for h (cid:48) and K (cid:48)(cid:48) . The details of this calculation are omitted on account ofspace. Step 3.
The E-L equation for K is given by (70) which is the same as (24). The proof of (25)involves a short calculation starting from (73), which is simplified to the form (25) by using (72). Remark 4:
The derivation of Euler-Lagrange equation, as presented above, is a heuristic onaccount of Step 1. A similar heuristic also appears in the original paper of Kushner [17]. There,the Kushner-Stratonovich PDE (22) is derived by considering a continuous-time limit of theBayes formula (15). The It ˆo’s rules are used to obtain the limit. Rigorous justification of thecalculation in Step 1, or its replacement by an alternate argument is the subject of future work.The calculation in Steps 2 and 3 require additional regularity assumptions on density p andfunction h : p is C and h is C . February 27, 2013 DRAFT9
E. Proof of Proposition 3.2.
Consider the ODE (24). It is a linear ODE whose unique solution is given by K ( x , t ) = p ( x , t ) (cid:18) C + C (cid:90) x − ∞ p ( y , t ) d y − σ W (cid:90) x − ∞ h ( y ) p ( y , t ) d y (cid:19) , (76)where the constant of integrations C = C = ˆ h t σ W because of the boundary conditions for K . Part 2 is an easy consequence of the minimum principle for elliptic PDEs [10]. F. Proof of Thm. 3.3
It is only necessary to show that with this choice of { u , K } , we have d p ( x , t ) = d p ∗ ( x , t ) ,for all x and t , in the sense that they are defined by identical stochastic differential equations.Recall d p ∗ is defined according to the K-S equation (22), and d p according to the forwardequation (23).If K solves the E-L BVP (24) then using (76), ∂∂ x ( p K ) = − σ W ( h − ˆ h ) p . (77)On multiplying both sides of (25) by − p , we have − up = ( h − ˆ h ) p K − σ W ( p K ) ∂ K ∂ x + ˆ hp K = − σ W ∂ ( p K ) ∂ x K − σ W ( p K ) ∂ K ∂ x + ˆ hp K = − σ W ∂∂ x ( p K ) + ˆ hp K , where we have used (77) to obtain the second equality. Differentiating once with respect to x and using (77) once again, − ∂∂ x ( up ) + σ W ∂ ∂ x ( p K ) = − ˆ h σ W ( h − ˆ h ) p . (78)Using (77)-(78) in the forward equation (23), we haved p = L † p + σ W ( h − ˆ h )( d Z t − ˆ h d t ) p . This is precisely the SDE (22), as desired.
February 27, 2013 DRAFT0
G. Proof of Thm. 3.5
The Gaussian density is given by: p ( x , t ) = √ π Σ t exp ( − ( x − µ t ) Σ t ) , (79)The density (79) is a function of the stochastic process µ t . Using Itˆo’s formula,d p ( x , t ) = ∂ p ∂ µ d µ t + ∂ p ∂ Σ d Σ t + ∂ p ∂ µ d µ t , where ∂ p ∂ µ = x − µ t Σ t p , ∂ p ∂ Σ = Σ t (cid:16) ( x − µ t ) Σ t − (cid:17) p , and ∂ p ∂ µ = Σ t (cid:16) ( x − µ t ) Σ t − (cid:17) p . Substituting these intothe forward equation (23), we obtain a quadratic equation Ax + Bx =
0, where A = d Σ t − (cid:18) α Σ t + σ B − γ Σ t σ W (cid:19) d t , B = d µ t − (cid:18) α µ t d t + γ Σ t σ W ( d Z t − γ µ t d t ) (cid:19) . This leads to the model (30) and (31).
H. BVP for Multivariable Feedback Particle Filter
Consider the multivariable linear system,d X t = α X t d t + σ B d B t (80a)d Z t = γ T X t d t + σ W d W t (80b)where X t ∈ R d , Z t ∈ R , α is an d × d matrix, γ is an d × { B t } is an d − dimensionalWiener process, { W t } is a scalar Wiener process, and { B t } , { W t } are assumed to be mutuallyindependent. We assume the initial distribution p ∗ ( x , ) is Gaussian with mean vector µ andvariance matrix Σ .The following proposition shows that the Kalman gain is a solution of the multivariableBVP (36), the Kalman gain solution does not equal the solution K g (see (40)), and that thesolution given by K g is not integrable with respect to p : Proposition 7.2:
Consider the d-dimensional linear system (80a)-(80b), where d ≥
2. Suppose p ( x , t ) is assumed to be Gaussian: p ( x , t ) = ( π ) d | Σ t | exp (cid:0) − ( x − µ t ) T Σ − t ( x − µ t ) (cid:1) , where x =( x , x , ..., x d ) T , µ t is the mean, Σ t is the covariance matrix, and | Σ t | > K ( x , t ) = σ W Σ t γ (81)2) Suppose that the Kalman gain K given in (81) is non-zero. For this solution to (36), theredoes not exist a function φ such that p K = ∇ φ . February 27, 2013 DRAFT1
3) Consider the solution K g ( x , t ) of the BVP (36) as given by (40). This gain function isunbounded: | K g ( x , t ) | → ∞ as | x | → ∞ , and moreover (cid:90) R d | K g ( x , t ) | p ( x , t ) d x = ∞ , (cid:90) R d | K g ( x , t ) | p ( x , t ) d x = ∞ . Proof:
The Kalman gain solution (81) is verified by direct substitution in the BVP (36)where the distribution p is Gaussian.The proof of claim 2 follows by contradiction. Suppose a function φ exists such that p K = ∇ φ ,then we have ∂ ( K p ) i ∂ x j = ∂ φ∂ x j ∂ x i = ∂ φ∂ x i ∂ x j = ∂ ( K p ) j ∂ x i , ∀ i , j (82)where ( K p ) i is the i th entry of vector K p . By direct evaluation, we have ∂ ( K p ) i ∂ x j = K i · (cid:0) Σ − t ( x − µ t ) (cid:1) j p . Using (82), we obtain K i ( Σ − t ) jk = K j ( Σ − t ) ik , ∀ i , j , k (83)Setting k = i , summing over the index i and using (81), we arrive attr ( Σ − t ) K = Σ − t K , where tr ( Σ − t ) denotes the trace of the matrix Σ − t . This provides a contradiction because K (cid:54)≡ Σ t is a positive definite symmetric matrix with | Σ t | > K g as given by (40): p K g ( x , t ) = σ W d ω d (cid:90) y − x | y − x | d ( h ( y ) − ˆ h ) p ( y , t ) d y For this integral, with h ( y ) ≡ γ T y , we have the following asymptotic formula for | x | ∼ ∞ , p K g ( x , t ) ∼ C | x | d + o ( | x | d ) , where C does not vary as a function of | x | (its value depends only upon the angular coordinates).For example, in dimension d = C is given by C ( x , x ) = C ( | x | cos ( θ ) , | x | sin ( θ )) = − d ω d (cid:32) cos ( θ ) sin ( θ ) sin ( θ ) cos ( θ ) (cid:33) Σ t γσ W , where Σ t γσ W is the Kalman gain vector.The result follows because p | x | d → ∞ and | x | d is not integrable in R d . Using the Cauchy-Schwarz inequality, (cid:90) | K g ( x , t ) | p ( x , t ) d x ≤ (cid:18) (cid:90) | K g ( x , t ) | p ( x , t ) d x (cid:19) , which shows that K g is not square-integrable. February 27, 2013 DRAFT2 R EFERENCES [1] A. Bain and D. Crisan.
Fundamentals of Stochastic Filtering . Springer, Cambridge, Mass, 2010.[2] D. P. Bertsekas and J. N. Tsitsiklis.
Neuro-Dynamic Programming . Atena Scientific, Cambridge, Mass, 1996.[3] A. Budhiraja, L. Chen, and C. Lee. A survey of numerical methods for nonlinear filtering problems.
Physica D: NonlinearPhenomena , 230(1-2):27 – 36, 2007.[4] A. J. Chorin. Numerical study of slightly viscous flow.
J. Fluid Mech. , 57:785–796, 1973.[5] D. Crisan and A. Doucet. A survey of convergence results on particle filtering methods for practitioners.
IEEE Trans.Signal Process. , 50(3):736–746, 2002.[6] D. Crisan and J. Xiong. Approximate McKean-Vlasov representations for a class of SPDEs.
Stochastics: An InternationalJournal of Probability and Stochastic Processes , pages 1–16, 2009.[7] Fred Daum and Jim Huang. Generalized particle flow for nonlinear filters. In
Proc. SPIE , pages 76980I–76980I–12, 2010.[8] A. Doucet, N. de Freitas, and N. Gordon.
Sequential Monte-Carlo Methods in Practice . Springer-Verlag, April 2001.[9] K. Doya, S. Ishii, A. Pouget, and R. P. N. Rao.
Bayesian Brain . Comput. Neurosci. MIT Press, Cambridge, MA, 2007.[10] L. C. Evans.
Partial Differential Equations . American Mathematical Society, 1998.[11] N. J. Gordon, D. J. Salmond, and A. F. M. Smith. Novel approach to nonlinear/non-Gaussian Bayesian state estimation.
IEE Proceedings F Radar and Signal Processing , 140(2):107–113, 1993.[12] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.
J. Comput. Phys. , 73:325–348, December 1987.[13] J. E. Handschin and D. Q. Mayne. Monte Carlo techniques to estimate the conditional expectation in multi-stage nonlinearfiltering.
International Journal of Control , 9(5):547–559, 1969.[14] M. Huang, P. E. Caines, and R. P. Malhame. Large-population cost-coupled LQG problems with nonuniform agents:Individual-mass behavior and decentralized ε -Nash equilibria. IEEE Trans. Automat. Control , 52(9):1560–1571, 2007.[15] G. Kallianpur.
Stochastic filtering theory . Springer-Verlag, New York, 1980.[16] H. Kunita.
Stochastic Flows and Stochastic Differential Equations . Cambridge University Press, Cambridge, 1990.[17] H. J. Kushner. On the differential equations satisfied by conditional probability densities of Markov processes.
SIAM J.on Control , 2:106–119, 1964.[18] A. Leonard. Vortex method for flow simulation.
Journal of Computational Physics , 37:289–335, 1980.[19] S. K. Mitter and N. J. Newton. A variational approach to nonlinear estimation.
SIAM Journal on Control and Optimization ,42(5):1813–1833, 2003.[20] B. Øksendal.
Stochastic differential equations (6th ed.): an introduction with applications . Springer-Verlag, Inc., NewYork, NY, USA, 2005.[21] A. K. Tilton, E. T. Hsiao-Wecksler, and P. G. Mehta. Filtering with rhythms: Application to estimation of gait cycle.
InProc. of American Control Conference , pages 3433–3438, June 2012.[22] J. Xiong. Particle approximations to the filtering problem in continuous time. In D. Crisan and B. Rozovskii, editors,
TheOxford Handbook of Nonlinear Filtering . Oxford University Press, 2011.[23] T. Yang, G. Huang, and P. G. Mehta. Joint probabilistic data association-feedback particle filter for multiple target trackingapplications.
In Proc. of American Control Conference , pages 820–826, June 2012.[24] T. Yang, P. G. Mehta, and S. P. Meyn. Feedback particle filter with mean-field coupling.
In Proc. of IEEE Conference onDecision and Control , pages 7909–7916, December 2011.[25] T. Yang, P. G. Mehta, and S. P. Meyn. A mean-field control-oriented approach to particle filtering.
In Proc. of AmericanControl Conference , pages 2037–2043, June 2011.[26] H. Yin, P. G. Mehta, S. P. Meyn, and U. V. Shanbhag. Synchronization of coupled oscillator is a game.
IEEE Trans.Automatic Control , 57(4):920–935, 2012., 57(4):920–935, 2012.