0:00:13thanks thanks very much to very nice to be here and
0:00:15thank you for organising
0:00:17special
0:00:18session um
0:00:20between talking about particle filters
0:00:23uh uh uh to sell not your filtering problem so uh
0:00:29problem is to estimate X
0:00:30which is D dimensional
0:00:32and were given a dynamical system which can be a linear
0:00:36in X
0:00:37as well as measurement sequence of measurements
0:00:39which can be nonlinear in in and
0:00:42the usual approach and particle filtering is to compute
0:00:46the entire are
0:00:47conditional density of X conditioned on all measurement
0:00:50which is a very very and vicious thing
0:00:53uh to to uh and hence the
0:00:57basic problem is computational complexity
0:01:00so here's a plot that we made about eight or nine years ago for the um bootstrap particle filter
0:01:07uh what i call the classical particle filter it um
0:01:11illustrates the problem in high dimensions
0:01:14uh
0:01:16uh a different dimensional state vectors a one dimensional problems to easy and high dimensional problems are hard it's the
0:01:22message of the chart the
0:01:24vertical axes is the error estimation error roll to the optimal so this is optimal performance
0:01:31and the horizontal axis is the number of particles
0:01:34uh a million particles uh
0:01:36ten dimensional state vector
0:01:39that ten orders of magnitude worse than optimal
0:01:42that's the problem now
0:01:45uh it's also been the focus of the research and particle filters for the last fifteen
0:01:51or uh
0:01:53so here's the
0:01:55the theory are we talking about today um
0:01:58fix this problem at least for a certain class of examples so here is the same
0:02:03problem uh for a thirty dimensional state vectors stable plant
0:02:08uh you get essentially optimal performance
0:02:11even for just the few that's was a hundred
0:02:13particle so if you contrast
0:02:16um
0:02:19uh ten dimensional state vector with a million particles ten orders of magnitude worse than optimal
0:02:25to the theory you'll be talking about today
0:02:28or uh about a four order of magnitude and prove meant
0:02:31at least for this example
0:02:33uh that example was a stable plant here is an extremely unstable plants meeting all the eigenvalues are in the
0:02:40left have
0:02:41in the right have
0:02:43plane and uh five dimensional state vector essentially optimal
0:02:48for just a hundred particles
0:02:49if the plant is uh a higher dimensional say thirty dimensional
0:02:54you need more particles to approach optimality
0:02:57uh
0:02:59but one shouldn't get too excited about that that's sort of result these are just linear system
0:03:05um it's an assurance that at least for linear system the new three L be talking about
0:03:11is interesting here is a highly nonlinear example
0:03:15hundred particles six dimensions
0:03:17the i'm talking about here
0:03:20each the extended kalman filter by roughly
0:03:23uh and order of magnitude so there's a hint
0:03:27that there something interesting useful new in that theory um talking about all show you more
0:03:32not an your examples later
0:03:35uh
0:03:36second question is the computational complexity the theory i'm talking about here which be call exact flow
0:03:44as compared to say the bootstrap particle filter and other particle filters um
0:03:51for are for different dimensional state vectors going from five to thirty
0:03:56uh for a hundred thousand particles
0:03:59the new theory is roughly four orders of magnitude faster that's running in matlab on a single
0:04:06P C uh and the question is i is that that
0:04:09quite surprising
0:04:11because we don't resample we have no proposal density and we never ever resample particles
0:04:16and the
0:04:17the bottleneck in a normal particle filter is the resampling step by eliminating
0:04:23that we've
0:04:24uh improve performance by
0:04:27about
0:04:27for orders of magnitude
0:04:29speed up so we have
0:04:31about three or four orders of magnitude fewer particles required in addition to
0:04:37three or four orders of magnitude faster execution of the code
0:04:42um
0:04:43so these put together
0:04:45is
0:04:46six to eight orders of magnitude faster computation for the same
0:04:51accuracy
0:04:52in addition to that
0:04:53um the theory all be telling you about this afternoon is uh
0:04:59embarrassingly parallel liza meaning that
0:05:03each individual particle has its own individual computation not interacting
0:05:09in a strong way with other particles uh
0:05:12this allows one to code these algorithms and say gpus thousand gpus you might give the speed up
0:05:19factor
0:05:20of a hundred this is unlike other particle filters where there's a very strong
0:05:24bottleneck because of the resample
0:05:27and again because we do not
0:05:29resample we've eliminated that
0:05:31well
0:05:32so um how does this work
0:05:35uh what's the theory
0:05:37um well actually the first question is what's the problem the problem is sort of surprising
0:05:43um any nonlinear filter would consist of these two parts a prediction of the probably density
0:05:50and the application of bayes rule
0:05:52so this is a pretty complicated operation it amounts to solving a partial differential equation
0:05:58where is this operation is uh not as simple as you can imagine is just multiplying two functions of
0:06:03point
0:06:04and the very interesting surprising
0:06:07but true
0:06:08situation is that that
0:06:10this is where the problem is
0:06:12uh for the most part
0:06:14well known
0:06:16oh it's been well known for many years and has goes by the name of particle the generous C and
0:06:21this is just a partition the double straight
0:06:24what it means uh
0:06:26if the prior density of the state is well represented at these points school particles
0:06:32but you wish to multiply the red curve times the blue curve it turns out that uh
0:06:37quite often it's unlikely to have particles in the right position
0:06:41to represent that problem of the red and the blue curve
0:06:44so somehow or other one would like to
0:06:47position the product
0:06:49we're get a new set of particles or moves the particles
0:06:53so was to be able to get got a good representation of the product
0:06:56and that's what we'll do
0:06:59a do it starting with an idea yeah that is uh uh not our our idea uh not a new
0:07:05idea uh it
0:07:06first occurs
0:07:08in a paper about ten years ago by simon got soul will be giving a talk later
0:07:14here uh
0:07:15it's the to find of flow of the probably city
0:07:18from the prior G
0:07:20to the post E the product of G times a
0:07:23or this single value parameter lambda a lambda is like time
0:07:27lan to go from zero to one
0:07:30from the prior to the posteriori
0:07:32and you can see when the is equal to one
0:07:35uh we
0:07:35P will be equal to the product of G times
0:07:39what you really like to do however is somehow or other computer flow of particle
0:07:44and that's what i'll tell you how to do flow the means
0:07:47move the particles
0:07:49using an ordinary differential equation
0:07:52uh from lambda at
0:07:53zero to lambda at one
0:07:56so the name of the game is to compute the flow
0:07:58a
0:07:59we can do that in a large number of examples and uh the way to do it is to uh
0:08:05expressed that uh flow as a function
0:08:09and tight turns out the functions related to the solution of this partial differential equation
0:08:15partial differential equation is uh
0:08:18stream really simple if you had to have a partial differential equation it's about the nice one
0:08:24to get it's linear
0:08:26it's first water and it has constant color fish
0:08:31um
0:08:32on the other hand
0:08:34we're particular about the solution that is we want the to be a stable dynamical system and secondly the left
0:08:40hand side is a known function
0:08:42but it's only known a random points in particular at the particle
0:08:47uh but we can solve but nevertheless uh using many different methods
0:08:51it's a um it's furthermore i should point out it's a
0:08:55say
0:08:57it's a
0:08:58it's a linear partial differential equation that's highly under determined
0:09:03at at it's a single scalar valued partial differential equation
0:09:06but it's in the known
0:09:08so uh
0:09:09is generically only the term
0:09:12uh there are many many ways to solve it and i'll show you a few
0:09:16in this brief
0:09:18period of time that i have
0:09:19the talk about um
0:09:22subject the first one is the most concrete the simplest to understand and the fast
0:09:27for a special case
0:09:28the special cases
0:09:30if the prior and the likelihood are both gaussian
0:09:33then we can solve the partial differential equation this function that's uh affine index
0:09:39the matrix a a
0:09:40and the vector V or explicitly computable
0:09:43and the very good news is that that dynamical flow
0:09:47is automatically stable
0:09:49well the eigenvalues are in the left half
0:09:52right
0:09:53a second simple explicit
0:09:55solution
0:09:57that's also quite fast is given here
0:10:00suppose you were given the divergence of F then you can pick the minimum norm solution F
0:10:05as the generalized inverse as gradient of the log of the uh
0:10:10the density times the uh log likelihood blessed virgin
0:10:15and uh although that's quite
0:10:17possibly
0:10:18come suitable um
0:10:20it turns out that the um
0:10:25yeah generalized inverse is nothing of this gradient is nothing more than the transpose divided by the square of the
0:10:31norm and so with
0:10:32like only amounts to this
0:10:34so the flow oh
0:10:36particles goes in the direction of the gradient
0:10:39of a lot of the density
0:10:41and uh if that
0:10:43gradient along the density is zero
0:10:46the flow stuff
0:10:47which is in agreement with ones
0:10:50intuition
0:10:51uh yeah the third solution
0:10:54which is fairly easy to grasp
0:10:56uh intuitively
0:10:57and quickly is that uh we can solve that P D we can turn the pde into plus sounds equation
0:11:03famous and Z
0:11:05by asking that the vector feel be the gradient of the scalar potential
0:11:10and then um to greens function for plus sounds equation in even in D dimensions as is well it's that
0:11:17differentiated to get the gradient
0:11:19we can pick either either of these two solution
0:11:22and it a
0:11:23turns into an algorithm
0:11:25with the algorithm is the following
0:11:27the direction of
0:11:28the flow of the uh particles is given as the sum over particles
0:11:33uh equal to the uh basically the electron density in this uh
0:11:38charged a particle feel
0:11:41times uh something that looks an awful lot like uh
0:11:44a well of D was equal to three this to be an inverse square
0:11:48law so it looks like cool ones a cool hoodlums law famous from physics
0:11:53uh in particular lecture min attics and fluid flow and also
0:11:57gravity
0:11:59so here we have particles which are mathematical constructs to solve a nonlinear filtering problem
0:12:05that should slow according to a famous physical law called cool homes
0:12:11law uh in addition to being interesting from a physical viewpoint it's very fast
0:12:17other ways to solve
0:12:19at uh
0:12:20concrete example
0:12:23for a highly nonlinear problem
0:12:25is
0:12:26to solve the problem where you were given the
0:12:29measurements of the cosine of an angle and then asked to estimate the and of course that's highly in figure
0:12:35a single measurement
0:12:36um
0:12:37so one has to depend on the motion i'll of
0:12:39the angle
0:12:40over time
0:12:42and the knowledge of that dynamics
0:12:44and um
0:12:45the dense the conditional density is highly multimodal uh bimodal at the beginning
0:12:51and so a kalman filter or extended kalman filter would die verge the error would grow with time
0:12:57um but if you wait about twenty seconds the particle filter that we design on version
0:13:04and B the extended kalman filter
0:13:06and uh the reason for this is of course extended kalman filters do not cannot
0:13:12represent multi-modal dance
0:13:14uh i here as a movie
0:13:16of the particle flow uh
0:13:19based on the very first measurement
0:13:21a part there's about a thousand particles it's dot as the particle
0:13:25it's in this uh as two dimensional state space
0:13:28and it there's a uniform distribution of uncertainty in angle and angle rate
0:13:34and as a result of making a single measurement
0:13:37what the differential equations that we derive
0:13:40automatically do
0:13:42is create that
0:13:43conditional them
0:13:45so we go from a uniform density to a bimodal
0:13:49then
0:13:51wondering what uh
0:13:53five but five minutes
0:13:56um
0:14:01most general solution for this uh um P D E can be written out this way
0:14:06terms the generalized inverse C sharp sharper of this when you're difference operator
0:14:11which is a matrix post the divergence operator
0:14:14a plus a projection and to the D minus one dimensions orthogonal to that
0:14:19uh
0:14:20direction and Y is anything
0:14:23so this is a very explicit
0:14:25simple representation of the
0:14:27um
0:14:29uh ambiguity of the arbitrariness in the solution
0:14:32a
0:14:33the solution of this equation is highly arbitrary um
0:14:37and
0:14:38a question is can you pick a Y
0:14:41that is a good
0:14:42uh solution
0:14:43take it good solution the answer is yes you can pick Y
0:14:47stabilise this flow to make that
0:14:49that's stable dynamical system
0:14:51and you can do it in a very very simple way can just picked white why the minus
0:14:56ten X
0:14:57as an exam
0:14:59uh quickly to nonlinear examples
0:15:02twelve dimensional state vector
0:15:04you measure three components the a quadratic measurements
0:15:08uh the error of the particle filter decreases as the number of particles increase
0:15:13or
0:15:14and you finally be the extended kalman filter performance by about two
0:15:19orders of magnitude similar for the cubic
0:15:22linearity
0:15:24one summarise by saying um
0:15:28like just repeating what i said
0:15:30uh are ready
0:15:32emphasising if you fact
0:15:34um
0:15:35speed is of the essence here
0:15:38um
0:15:39and what i've channel you is an is a new algorithm based on a new theory that's orders of magnitude
0:15:45faster than
0:15:46certainly the bootstrap filter
0:15:48and other standard particle filter
0:15:52at least for certain non linear examples that i've shown you it's orders of magnitude more accurate than the extended
0:15:59kalman filter
0:16:01uh it focuses on solving
0:16:03the key
0:16:04a problem with standard particle filters that is
0:16:08the problem particle to C
0:16:10and it does it
0:16:12uh by moving the particles to the right
0:16:15part of state space
0:16:16to compute
0:16:17the product of two functions corresponding to bases
0:16:21rule
0:16:23we never resample particles because we don't have to and we'd rather not because that's a time consuming operation
0:16:30since we never resample we never need a proposal density
0:16:34uh we don't use any
0:16:36a markov chain monte carlo that method
0:16:39uh we never read represent the probably density but rather we represent the log
0:16:44of the a normalized probably density
0:16:47i've asserted that this is
0:16:48highly parallel lies able
0:16:50and uh
0:16:51you might believe need but if you don't you can
0:16:54try for yourself when your favourite seven
0:16:57you use
0:16:58and a final point is um
0:17:00a a is this magic were fraud or is there some reason
0:17:05that one can grasp intuitively Y
0:17:09we should be getting better performance
0:17:11uh and i think
0:17:12i hint at
0:17:14a reason is
0:17:15where lot solving the most general filtering problem
0:17:18we are not solving wall of the
0:17:20i only full during problems that particle filters could solve
0:17:24we solving the special
0:17:26class of
0:17:27problems where the probably densities are nowhere vanishing
0:17:31and continuously french
0:17:33the state
0:17:34which is a more restricted class
0:17:36so that we can write differential equation
0:17:40and exploit the good idea is uh a mathematicians of
0:17:45cooked top over the last two hundred years
0:17:48to a van
0:17:50a at least that's my curve
0:17:52maybe be think something else after uh
0:17:54a few years of reflection
0:17:56questions uh please
0:18:07simon
0:18:16uh nowhere vanishing
0:18:19yeah
0:18:21a course in terms of computational complexity
0:18:23the structure of the problem um um
0:18:27can greatly influence
0:18:28as you well know the um resulting computational complexity that from a theoretical
0:18:33you point a it is nowhere vanishing continuously differentiable
0:18:37then sit
0:18:42peter
0:19:02the particles are asked to follow the conditional density
0:19:06uh throughout the predict time prediction step
0:19:09as in any particle filter and the particles are asked to follow the conditional density
0:19:15uh at the big
0:19:17but four and after bayes rule
0:19:19so if the algorithm is doing what i've asked to do the particles do in fact
0:19:24satisfy the conditional then
0:19:27yeah every single point
0:19:28in the L
0:19:29and at each
0:19:30step
0:19:32there's no proof that that happen
0:19:34i have no idea what these sufficient conditions would be on such or
0:19:38proof obviously you need enough particles
0:19:41you need an of um
0:19:44regularity of the densities and the dynamics
0:19:47and the likelihood
0:19:48uh we're no close to having a theory
0:19:51that would be useful in guaranteeing such a result
0:19:54it's proof by map
0:19:57on a honest on the subset of
0:19:59examples that we chose
0:20:19for for all examples that we've were we we have worked out uh six different classes of
0:20:24exam
0:20:26lots of linear examples
0:20:27uh
0:20:28quadratic
0:20:29Q
0:20:30uh rate examples both with and without range measurements
0:20:34and the um
0:20:36measurements of phase
0:20:38the cosine
0:20:39so when all those examples that we
0:20:41run our
0:20:42we get
0:20:43better
0:20:44much better
0:20:45accuracy than the extended
0:20:47comes to
0:20:52right
0:20:55for the linear K yes for the linear your cases um
0:20:58and for the quadratic and Q
0:21:01for a given number
0:21:03particles we be the
0:21:05uh
0:21:07well
0:21:08a selected
0:21:09set of standard particle
0:21:12you probably know that there's lots and lots of particle filter
0:21:14and will probably hear about new ones
0:21:17sept F known that i
0:21:18that one of that so
0:21:20an exhaustive check be
0:21:22the
0:21:23uh
0:21:24uh suppose you not believe me i'll give you the code and you can check it on your very for
0:21:28uh set of
0:21:29particle
0:21:32a point of comparison
0:21:33i think that's a way to make progress
0:21:38or questions
0:21:43yeah