Stochastic filtering via L projection on mixture manifolds with computer algorithms and numerical examples

We examine some differential geometric approaches to finding approximate solutions to the continuous time nonlinear filtering problem. Our primary focus is a projection method using the direct L2 metric onto a family of normal mixtures. We compare this method to earlier projection methods based on the Hellinger distance/Fisher metric and exponential families, and we compare the L2 mixture projection filter with a particle method with the same number of parameters. We study particular systems that may illustrate the advantages of this filter over other algorithms when comparing outputs with the optimal filter. We finally consider a specific software design that is suited for a numerically efficient implementation of this filter and provide numerical examples.


rential equat
on.The observations that have been made contain an additional noise term, so one cannot hope to know the true state of the system.However, one can reasonably ask what is the probability density over the possible states.

When the observations are made in continuous time, the probability density follows a stochastic partial differential equation k own as the Kushner-Stratonovich equation.This can be seen as a generalization of the Fokker-Planck equation that expresses the evolution of the density of a diffusion process.Thus the problem we wish to address boils down to finding approximate solutions to the Kushner-Stratonovich equation.

For a quick introduction to the filtering problem see Davis and Marcus (1981) [15].For a more complete treatment from a mathemati al point of view see Lipster and Shiryayev (1978) [29].See Jazwinski (1970) [21] for a more applied perspective.For recent results see the collection of papers [14].

The main idea we will employ is inspired by the differential geometric approach to statistics developed in [2] and [32].One thinks o the probability distribution as evolving in an infinite dimensional space P which is in turn contained in some Hilbert space H.One can then think of the Kushner-Stratonovich equation as defining a vector field in P: the integral curves of the vector field should correspond to the solutions of the equation.To find approximate solutions to the Kushner-Stratonovich equation one chooses a finite dimensional submanifold M of H and approximates the probability distributions as points in M .At each point of M one can use the Hilbert space structure to project the vector field onto the tangent space of M .One can now attempt to find approximate solutions to the Kushner-Stratonovich equations by integrating this vector field on the manifold M .This mental image is slightly innaccurate.The Kushner-Stratonovich equation is a stochastic PDE rather than a PDE so one should imagine some kind of stochastic vector field rather than a smooth vector field.Thus in this approach we hope to approximate the infinite dimensional stochastic PDE by solving a finite dimensional stochastic ODE on the manifold.Note that our approximation will depend upon two choices: the choice of manifold M and the choice of Hilbert space structure H.In this paper we will consider two possible choices for the Hilbert space structure: the direct L 2 metric on the space of probability distributions; the Hilbert space structure associated with the Hellinger distance and the Fisher Inf rmation metric.Our focus will be on the direct L 2 metric since projection using the Hellinger distance has been considered before.As we shall see, the choice of the "best" Hilbert space structure is determined by the manifold one wishes to consider -for manifolds associated with exponential families of distributions the Hellinger metric leads o the simplest equations, whereas the direct L 2 metric works well with mixture distributions.

We will write down the stochastic ODE determined by this approach when H = L 2 and show how it leads to a numerical scheme for finding approximate solutions to the Kushner-Stratonovich equations in terms of a mixture of normal dis ributions.We will call this scheme the L 2 normal mixture projection filter or simply the L2NM projection filter.

The stochastic ODE for the Hellinger metric was considered in [12], [9] and [10].In particular a precise numerical scheme is given in [12] for finding solutions by projecting onto an exponential family of distributions.We will call this scheme the Hellinge exponential projection filter or simply the HE projection filter.

We will compare the results of a C++ implementation of the L2NM projection filter with a number of other numerical approaches including the HE projection filte and the optimal filter.We can measure the goodness of our filtering approximations thanks to the geometric structure and, in particular, the precise metrics we are using on the spaces of probability measures.

What emerges is that the two projection methods produce excellent resul s for a variety of filtering problems.The results appear similar for both projection methods; which gives more accurate results depends upon the problem.

As we shall see, however, the L2NM projection approach can be implemented more efficiently.In particular one needs to perform numerical integration as part of the HE projection filter algorithm whereas all integrals that occur in the L2NM projection can be evaluated analytically.

We also compare the L2NM filter to a particle filter with the best possible combination of particles with respect to the Lévy metric.Introducing the Lévy metric is needed because particles densities do not compare well with smooth densities when using L 2 induced metrics.We show that, given the same number of parameters, the L2NM may outperform a particles based system.

The paper is structured as follows: In Section 2 we introduce the nonlinear filtering problem and the infinite-dimensional Stochastic PDE (SPDE) that solves it.In Section 3 we introduce the geometric structure we need to project the filtering SPDE onto a finite dimensional manifold of probabi

ty densities.In Section 4 we perform the projection of the filterin
SPDE according to the L2NM framework and also recall the HE based framework.In Section 5 we briefly discuss the numerical implementation, while in Section 6 we discuss in detail the software design for the L2NM filter.In Section 7 we look at numerical results, whereas in Section 8 we compare our outputs with a particle method.Section 9 concludes the paper.


The non-linear filtering pr some system is modelled by a process X called the signal.This signal evolves over time t according to an Itô stochastic differential equation (SDE).We measure the state of the system using some observation Y .The obs rvations are not accurate, there is a noise term.So the observation Y is related to the signal X by a second equation.
dX t = f t (X t ) dt + σ t (X t ) dW t , X 0 , dY t = b t (X t ) dt + dV t , Y 0 = 0 .(1)
In these equations the unobserved state process {X t , t ≥ 0} takes values in R n , the observation {Y t , t ≥ 0} takes values in R d and the noise processes {W t , t ≥ 0} and {V t , t ≥ 0} are two Brownian motions.

The nonlinear filtering problem consists in find nal probability distribution π t of the state X t given the observations up to time t and the prior distribution π 0 for X 0 .

Let us assume that X 0 , and the two Brownian motions are independent.Let us also assume that the covariance matrix for V t is invertible.We can th t defined by:
a t = σ t σ T t
With these preliminaries, and a number o , one can show that π t satisfies the a stochastic PDE called the Kushner-Stratonovich equation.This states that for any compactly supported test function φ defined on R n
π t (φ) = π 0 (φ)+ t 0 π s (L s φ) ds+ d k=1 t 0 [π s (b k s φ)−π s (b k s ) π s (φ)] [dY k s −π s (b k s ) ds] ,(2)
where for all t ≥ 0, the backward diffusion operator L t is defined by
L t = n i=1 f i t ∂ ∂x i + 1 2 n i,j=1 a ij t ∂ 2 ∂x i ∂x j .
Equation ( 2 We assume now that π t can be represented by a density p t with respect to the Lebesgue measure on R n for all time t ≥ 0 and that we can replace the term involving L s φ with a term involving its formal adjoint L * stochastic partial differential equation (SPDE):
dp t = L * t p t dt + d k=1 p t [b k t − E pt {b k t }][dY k t − E pt {b k t }dt]
where E pt {•} denotes the expectation with respect to the probability density p t (equivalently the conditional expectation given the observations up to time t).The forward diffusion operator L * t is defined by:
L * t φ = − n i=1 ∂ ∂x i [f i t φ] + 1 2 n i,j=1 ∂ 2 ∂x is necessary to use Stratonovich SDE's rather than Itô SDE's.This is because one does not in general know how to interpret the second order terms that arise in Itô calculus in terms of manifolds.The interested reader should consult [16].A straightforward calculation yields the following Stratonvich SPDE:
dp t }] • dY k t .
We have indicated that this is the Stratonovich form of the equatio nd the Brownian motion of the SDE.We shall use this convention throughout the rest of the paper.

In order to simplify notation, we introduce the following definitions :
γ 0 t (p) := 1 2 [|b t | 2 − E p {|b t | 2 }] p, γ k t (p) := [b k t − E p {b k t }]p ,(3)for k = 1, • • • , d.
The Str form of the Kushner-Stratonovich equation reads now
dp t = L * t p t dt − γ 0 t (p t ) dt + d k=1 γ k t (p t ) • dY k t .(4)
Thus, subject to the assumption that a density t exists for all time and assuming the necessary decay condition to ensure tha the non-linear filtering problem is equivalent to solving thi SPDE.Numerically approximating the solution of e we review the technical conditions r quired in order for equation 2 to follow from (1).

(A) Local l R > 0, there exists K R > 0 such th t
|f t (x)−f t (x )| ≤ K R |x−x | and a t (x)−a t (x ) ≤ K R |x−x | ,
for all t ≥ 0, and for all x, x ∈ B R , the ball of radius R.

(B) Non-explosion : there exists K > 0 such that
x T f t (x) ≤ K (1 + |x| 2 ) and trace a t (x) ≤ K (1 + |x| ll x ∈ R n .

( ) Polynomial growth : there exist K > 0 and r ≥ 0 such that
|b t (x)| ≤ K (1 + |x| r ) ,
for all t ≥ 0, and for all x ∈ R n .

Under assumptions (A) and (B), ther der the additional assumption (C) the following finite energy condition holds < ∞ ,
for all T ≥ 0.

Since the finite energy condition holds, it follows from Fujisaki, Kallianpur and Kunita [22] that {π t , t ≥ 0} satisfies the Kushner-Stratonovich equation 2.

3 Statistical manifolds


Families of distributions

As discussed in the introduction, the idea of a projection filter is to approximate soluti ns to the Kushner-Stratononvich equation 2 using a finite dimensional family of distributions.

Example 3.1 A normal mixture family contains distributions given by:
p = m i=1 λ i 1 σ i √ 2π exp −(x − µ i ) 2 2σ 2 i with λ i > 0 and λ i = 1. It is a 3m − 1 dimensional family of distributions.
Example 3.2 A polynomial exponential family contains distributions given by:
p = exp( m i=0 a i x i )
where a 0 is chosen to ensure that the integral of p is equal to 1.To ensure the convergence of the integral we must have that m is even and a m < 0. This is an m dimensional family of distributions.Polynomial exponential families are a special case of the more general notion of an exponential family, see for example [2].

A key motivation for considering these families is that one can reproduce many of the qualitative features of distributions that arise in practice using these distributions.For example, consider the qualitative specification: the distribution should be bimodal with peaks near −1 and 1 with the peak at −1 twice as high and twice as wide as the peak near 1.One can easily write down a distribution of this approximates form using a normal mixture.

To find a similar exponential family, one seeks a polynomial with: local maxima at −1 and 1; with the maximum values at these points differing by log(2); with second derivative at 1 equal to twice that at −1.These conditions give linear equations in the polynomial coeff

ients.Using degree 6 polynomials it is simple to
find solutions meeting all these requirements.A specific numerical example of a polynomial meeting these requirements is plotted in Figure 1.The associated exponential distribution is plotted in Figure 2. We see that normal mixtures and exponential families have a broadly similar power to describe the qualitative shape of a distribution using only a small number of parameters.Our hope is that by approximating the probability distributions that occur in the Kushner-Stratonovich equation by elements of one of these families we will be able to derive a low dimensional approximation to the full infinite dimensional stochastic partial differential equation.


Two Hilbert spaces of probability distributions

We have given direct parameterisations of our famil es of probability distributions and thus we have implicitly represented them as finite dimensional manifolds.In this section we will see how families of probability distributions can be thought of as being embedded in a Hilbert space and hence they inherit a manifold structure an metric from this Hilbert space.There are two obvious ways of thinking of embedding a probability density function on R n in a Hilbert space.The first is to simply assume that the probability density function is square integrable and hence lies directly in L 2 (R n ).The second is to use the fact that a probability density function lies in L 1 (R n ) and is non-negative almost everywhere.Hence √ p will lie in L 2 (R n ).

For clarity we will write L 2 D (R n ) when we think of L 2 (R n ) as containing densities directly.The D stands for direct.We write D ⊂ L 2 D (R n ) where D is the set of square integrable probability densities (functions with integral 1 which are positive almost everywhere).

Similarly we will write L 2 H (R n ) when we think of L 2 (R n ) as being a space of square roots of densities.The H stands for Hellinger (for reasons we will explain shortly).We will write H for the subset of L 2 H consisting of square roots of probability densities.

We now have two possible ways of formalizing the notion of a family of tions to be either a smooth submanifold of L 2 D which also lies in D or a smooth submanifold of L 2 H which also lies in H. Either way the families we discussed earlier will give us finite dimensional families in this more formal sense.

The Hilbert space structures of L 2 D and L 2 H allow us to define two notions of distance between probability distributions which we will denote d D and d H .Given two probability distributions p 1 and p 2 we have an injection ι into L 2 so one defines the distance to be the norm of ι(p 1 ) − ι(p 2 ).So given two probability densities p 1 and p 2 on R n we can define:
d H (p 1 , p 2 ) = ( √ p 1 − √ p 2 ) 2 dµ 1 2 d D (p , p 2 ) = (p 1 − p 2 ) 2 dµ 1 2
.

Here dµ is the Lebesgue measure.d H defines the Hellinger distance between the two distributions, which explains are use of H as a subscript.We will write •, • H for the inner product associated with d H and •, • D or simply

•, • for the inner product associated with d D .

In this paper we will consider the projection of the conditional density of the true state of the system given the observations (which is assumed o lie in D or H) onto a submanifold.The notion of projection only makes sense with respect to a particular inner product structure.Thus we can consider projection using d H or projection using d D .Each has advantages and disadvantages.

The most notable advantage of the Hellinger metric is that the d H metric can be defined independently of the Lebesgue measure and its definition can be extended to define the distance between measures without density functions (see Jacod and Shiryaev [20] or Hanzon [18]).In particular the Hellinger distance is indepdendent of the choice of parameterization for R n .This is a very attractive feature in terms of the differential geometry of our set up.

Despite the significant theoretical advantages of the d H metric, the

D metric has an obvious advantage when studyin
mixture families: it comes from an inner product on L 2 D and so commutes with addition on L 2 D .So it should be relatively easy to calculate with the d D metric when adding distributions as happens in mixture families.As we shall see in practice, when one performs concrete calculations, the d H metric works well for exponential families and the d D m lies.While the d H metric leads to the Fisher Information and to an equivalence with Assumed Density Filters when used on exponential families, see [10], the d D metric for simple mixture families is equivalent to a Galerkin method, see for example [8].


The tangent space of a family of distributions

To make our notion of smooth families precise we need to explain what we mean by a smooth map into an infinite dimensional space.

Let U and V be Hilbert spaces and let f : U → V be a continuous map (f need only be defined on some open subset of U ).We sa that f is Frećhet differentiable at x if there exists a bounded linear map A : U → V satisfying:
lim h→x f (h) − f (x) − Ah V h U
If A exists it is unique and we d note it by Df (x).This limit is called the Frećhet derivative of f at x.It is the best linear approximation to f at 0 in the se be an infinitely Frećhet differentiable map.We define an immersion of an open subset of R n into V to be a map such that Df ( int where f is defined.The latter condition ensures that the best linear iven an immersion f defined on a neighbourhood of x, we can think of the vector subspace of V given by the image of Df (x) as representing the tangent space at x.

To make these ideas more concrete, let us suppose that p(θ) is a probability distribution depending smoothly on some parameter θ
= (θ 1 , θ 2 , . . . , θ m ) ∈ U where U is some open subset of R m . The map θ → p(θ) defines a map i : U → D. At a given point θ ∈ U and for a vector h = (h 1 , h 2 , . . . , h m ) ∈ R m we can compute the Fréchet derivative to obtain: Di(θ)h = m i=1 ∂p ∂θ i h i
So we can identify the tangent space at θ with the following subspace of
L 2 D : span{ ∂p ∂θ 1 , ∂p ∂θ 2 , . . . , ∂p ∂θ m }(5)
We can formally define a s an immersion of an open subset of R n into D. Equivalently it is a smoothly parameterized probability distribution p such that the above vectors in L 2 are linearly independent.

We can define a smooth m-dimensional family of probability distributions in L 2 H in the same way.This time let q(θ) be a square root of a probability distribution depending smoothly on θ.The tangent vectors in this case will be the partial derivatives of q with respect to θ.Since one normally prefers to work in terms of probability distributions rather than their square roots we use the chain rule to write the tangent space as: pan{ 1 2 √ p ∂p ∂θ 1 , 1 2 √ p ∂p ∂θ 2 , . . . , 1 2 √ p ∂p ∂θ m }(6)
We have defined a family of distributions in terms of a single immersion f into a Hilbert space V .In other words we have defined a family of distributions in terms of a specific parameterization of the image of f .It is tempting to try and phrase the theory in terms of the image of f .To this end

one defines an embedded subman
fold of V to be a subspace of V which is covered by immersions f i from open subsets of R n where each f i is a homeomorphisms onto its image.With this definition, we can state that the tangent space of an embedded submanifold is independent of the choice of parameterization.

One might be tempted to talk about submanifolds of the space of probability distributions, but one should be careful.The spaces H and D are not open subsets of L tly by subtracting a small delta-like function.


The Fisher information metric

Given two tangent vectors at a point to a family of probability distributions we can form their inner product using •, • H .This defines a so-called Riemannian metric on the family.With respect to a particular parameterization θ we can compute the inner product of the i th and j th basis vectors given in Figure 3: An element of L 2 close to the normal distribution but not in H equation 6.We call this quantity 1 4 g ij .
1 4 g ij (θ) := 1 2 √ p ∂p ∂θ i , 1 2 √ p ∂p ∂θ j H = 1 4 1 p log p ∂θ i ∂ log p ∂θ j pdµ = 1 4 E p ( ∂ log p ∂θ i ∂ log p ∂θ j )
Up to the factor of 1  4 , this last formula is the standard definition for the Fisher information matrix.So our g ij is the Fisher information matrix.We can now interp ion metric and observe that, up to the constant factor linger distance.See [2], [30] and [1] for more in depth study on this dif .3

The Gaussian family of densities can be parameterized using parameters mean µ and variance 0 1/(2v)
The representation of the metric as a matrix depends heavily upon the choice of parameterization for the family.


Example 3.4

The Gaussian family may be considered as a particular exponential family with parameters θ 1 and θ 2 given by:
p(x, θ) = exp(θ 1 x + θ 2 x 2 − ψ(θ))
where ψ(θ) is chosen to normalize 2 − θ 1 2 4θ 2
This is related to the familiar parameterization in terms of µ and v by:
µ = −θ 1 /(2θ 2 ), v = σ 2 = (1/θ 2 − θ 2 1 /θ 2 2 )/2
One can c mpute the Fisher information metric relative to the parameterization θ 1 to obtain:
g(θ) = −1/(2θ 2 ) θ 1 /(2θ 2 2 ) θ 1 /(2θ 2 2 ) 1/(2θ 2 2 ) − θ 2 1 /(2θ3
2 )

The particular importance of the metric structure for this paper is t at one has m linearly independent vectors w i spanning some subspace W of a Hilbert space V .By linearity, one can write the orthogonal projection onto W as:
Π(v) = m i=1 [ m i=1 A ij v, w j ]w i
for some appropriately chosen constants A ij .Since acts as the identity on w i we see that A ij must be the inverse of the ma apply this to the basis given in equation 6. Defining g ij to be the inverse of the matrix g ij we obtain the following formula for projection, using the Helli stributions:
Π H (v) = m i=1 m j=1 4g ij v, 1 2 √ p ∂p ∂θ j H 1 2 √ p ∂p ∂θ i 3.5 The direct L 2 metric
The ideas from the previous section can also be applied to the direct L 2 metric.This gives a different Riemannian metric on the manifold.We will write h = h ij to denote the L 2 metric when written with respect to a particular parameterization.

Example 3.5 In coordinates µ, ν, the L 2 metric on the Gaussian family is:
h(µ, ν) = 1 8ν
√ νπ
1 0 0 3 4ν
We can obtain a formula for projection in L 2 D using the direct L 2 metric using the basis given in equation 5. We write h ij for the matrix inverse of
h ij . Π D (v) = m i=1 m j=1 h ij v, ∂p ∂θ j D ∂p ∂θ i .(7)

The projection filter

Given a family of probability distributions parameterised by θ, we wish to approximat a stochastic vector field in D and then project that vector field onto the tangent space of our family.The projected equations can then be viewed as giving a stochastic differen r our family.A curve t → θ(t) in the parameter space corresponds to a c tion 4 can be written:
d t p(•, θ(t)) = m i=1 ∂p(•, θ(t)) ∂θ i d t θ i (t) = m i=1 v i dθ i
where .

Given the projection formula given in equation 7, we can project the terms on the right hand side onto the ta gent space of the manifold using the direct L 2 metric as follows:
Π θ D [L * p] = m i=1 m j=1 h ij L * p, v j v i = m i=1 m j=1 h ij p, Lv j v i Π θ D [γ k (p)] = m i=1 m j=1 h ij γ k (p), v j v i θ i = m i=1 m j=1 h ij p, Lv j dt − γ 0 (p), v j dt + d k=1 γ k (p), v j • dY k v i
Since the v i form a basis of the tangent space, we can equate the coefficients of v

to obtain:
dθ i = m j=1 h
ij p(θ), Lv j dt − γ 0 (p(θ)), v j dt + d k=1 γ k (p(θ)), v j • dY k .
(8) This is the promised finite dimensional stochastic differential equation for θ corresponding to L 2 projection.

If preferred, one could instead project the Kushner-Stratonovich equation using the Hellinger metric instead.This yields the follo ing stochastic differential equation derived originally in [12]:
dθ i = m j=1 g ij L * p p , v j dt − 1 2 |b| 2 , v j dt + d k=1 b k , v j • dY k(9)
Note that the inner products in this equation are the direct L 2 inner products: we are simply using the L 2 inner product notation as a compact notation for integrals.


Numerical implementation

Equations 8 and 9 both give finite dimensional stochastic differential equations that we hope will approximate well the solution to the full Kushner-Stratonovich equation.We wish to solve these finite di

nsional equation
numerically and thereby obtain a numerical approximation to the non-linear filtering problem.

Because we are solving a low dimensional system of equations we hope to end up with a more efficient scheme than a brute-force finite difference approach.A finite difference approach can also be seen as a reduction of the problem to a finite dimensional system.However, in a finite difference approach the finite dimensional system still has a very large dimension, determined by the number of grid points into which one divides R n .By contrast the finite dimensional manifolds we shall consider will be defined by only a handful of parameters.


Software design

The specific solution algorithm will depend upon numerous choices: whether to use L 2 or Hellinger projection; which family of probability distributions to choose; how to parameterize that family; the representation of the functions f , σ and b; how to perform the integrations which arise from the calculation of expectations and inner products; the numerical method selected to solve the finite dimensional equations.

To test the effectiveness of the projection idea, we have implemented a C++ engine which performs the numerical solution of the finite dimensional equations and allows one to make various selections from the options above.Currently our implementation is restricted to the case of the direct L 2 projection for a 1-dimensional state X and 1-dimensional noise W .However, the engine does allow one to experiment with various manifolds, parameteriziations and functions f , σ and b.

We use object oriented programming techniques in order to allow this flexibility.Our implementation contains two key classes FunctionRing and Manifold.

To perform the computation, one must choose a data structure to represent elements of the function space.
owever, the most effective choice of representation depends upon the family of probability distributions one is considering and the functions f , σ and b.Thus the C++ engine does not manipulate the data structure directly but instead works with the functions via the FunctionRing in erface.A UML (Unified Modelling Language [33]) outline of the FunctionRing interface is given in table 1.The other key abstraction is the Manifold.We give a UML representation of this abstraction in table 2. For readers unfamiliar with UML, we remark that the * symbol can be read "list".For example, the computeTangentVectors function returns a list of functions.


FunctionRing

The Manifold uses some convenient internal representation for a point, the most obvious representation being simply the m-tuple (θ 1 , θ 2 , . . .θ m ).On request the Manifold is able to provide the density associated with any point represented as an element of the FunctionRing.

In addition the Manifold can compute the tangent vectors at any point.The computeTangentVectors method returns a list of elements of the Func-tionRing corresponding to each o the vectors v i = ∂p ∂θ i in turn.If the point is represented as a tuple θ = (θ 1 , θ 2 , . . .θ n ), the method updatePoint simply adds the components of the tuple ∆θ to each of the components of θ.If a different internal representation is used for the point, the method should make the equivalent change to this internal representation.

The finalizePoint method is called by our algorithm at the end of every time step.At this point the Manifold implemen ation can choose to change its parameterization for the state.Thus the finalizePoint allows us (in principle at least) to use a more sophisticated atlas for the manifold than just a single chart.

One should not draw too close a parallel between these computing abstractions and similarly named mathematical abstractions.For example, the space of objects that can be represented by a given FunctionRing do not need to form a differential ring despite the diffe

ntiate method.This is bec
use the differentiate function will not be called infinitely often by the algorithm below, so the functions in the ring do not ne d to be infinitely differentiable.

Similarly the finalizePoint method allows the Manifold implementation more flexibility than simply changing chart.From one time step to the next it could decide to use a completely different family of distributions.The interface even allows the dimension to change from one time step to the next.We do not currently take advantage of this possibility, but adapatively choosing the family of distributions would be an interesting topic for further research.


Outline of the algorithm he C++ engine is initialized with a Manifold object, a copy of the initial Point and Function objects representing f , σ and b.

At each time point the engine asks the manifold to compute the tangent vectors given the current point.Using the multiply and integrate functions of the class FunctionRing, the engine can compute the inner products of any two functions, hence it can compute the etric matrix h ij .Similarly, the engine can ask the manifold for the density function given the current point and can then compute L * p. Proceeding in this way, all the coefficients of dt and •dY in equation 8 can be computed at any given point in time.

Were equation 8 an Itô SDE one could now nu n θ over a given time interval ∆ in terms of ∆ and ∆Y , the change in Y .One would then use the updateState method to equation 8 an Itô SDE we could numerically solve the SDE using the Euler scheme.

However, equation 8 is a Stratonovich SDE so the Euler scheme is no longer valid.Various numerical schemes for solving stochastic differential equations are considered in [13] and [25].One of the simplest is the Stratonovich-Heun method described in [13].Suppose that one wishes to solve the SDE:
dy t = f (y t )dt + g(y t ) • dW t
The Stratonvich-Heun method generates an estimate for the solution y n at the n-th time interval using the formulae:
Y n+1 = y n + f (y n )∆ + g(y n )∆W n y n+1 = y n + 1 2 (f (y n ) + f (Y n+1 ))∆ + 1 2 (g(y n ) + g(Y n+1 ))∆W n
In these formulae ∆ is the size of the time interval and ∆W n is the chang in W .One can think of Y n+1 as being a prediction and the value y n+1 as being a correction.Thus this scheme is a direct translation of the standard Euler-Heun scheme for ordinary differential equations.We can use the Stratonovich-Heun method to numerically solve equatio 8. Given the current value θ n for the state, compute an estimate for ∆θ n by replacing dt with ∆ and dW with ∆W in equation 8. Using the updateState method compute a prediction Θ n+1 .Now compute a second estimate for ∆θ n using equation 8 in the state Θ n+1 .Pass the average of the two estimates to the updateState function to obtain the the new state θ n+1 .

At the end of each time step the method finalizeState is called.This provides the manifold implementation the opportunity to perform checks such as validation of the state, to correct the normalization and, if desired, to change the representation it uses for the state.

One small observation worth making is that the equation 8 contains the term h ij , the inverse of the matrix h ij .However, it is not necessary to actually calculate the matrix inverse in full.It is better numerically to multiply both sides of equation 8 by the matrix h ij and then compute dθ by solving the resulting linear equations directly.This is the approach taken by our algorithm.

As we have already observed, there is a wealth of choices one could make for the numerical scheme used to solve equation 8, we have simply selected the most convenient.The existing Manifold and FunctionRing implementations co

d be used directly by many of these schemes
-in particular those based on Runge-Kutta schemes.In principle one might also consider schemes that require e or higher derivatives such as ∂ 2 p ∂θ i ∂θ j .In this case one would need to extend the manifold abstraction to provide this information.

Similarly one could use the same concepts in order to solve equation 9 where one uses the Hellinger projection.In this case the FunctionRing would need to be extended to allo division.This would in turn complicate the implementation of the integrate function, which is why we have not yet i plemented this approach.


Implementation for normal mixture families

Let R denote the space of functions which can be written as finite linear combinations of terms of the form:
±x n e ax 2 +bx+c
where n is non-negative integer and a, b and c are constants.R is clo ed under addition, multiplication and differentiation, so it forms a differential ring.

We have written an implementation of FunctionRing corresponding to R.Although th Firstly, we store elements of our ring in memory as a collection of tuples (±, a, b, c, n).Although one can write:

±x n e ax 2 +bx+c = qx n e ax 2 +bx for appropriate q, the use or such a term in computer memory should be avoided as it will rapidly lead to significant rounding errors.A small amount of care is required throughout the implementation to avoid such rounding errors.

Secondly let us consider e plicitly how to implement integration for this ring.Let us define u n to be the integral of x n e −x 2 .Using integration by parts one has:
u n := ∞ −∞ x n e −x 2 dx = n − 1 2 ∞ −∞ x n−2 e −x 2 dx = n − 1 mpute u n recursively titution, we can now integrate p(x − µ)e −(x−µ) 2 for any µ.By completing the square we can analytically compute the integral of p(x)e ax 2 +bx+c so long as a < 0. Putting all this together one has an algorithm for analytically integrating the elements of R.
Let N i denote the space of probability distributions that can be written s i k=1 c k e a k x 2 +b k x for some real numbers a k , b k and c k with a k < 0. Given a smooth curve γ(t) in N i we can write:
γ(t) = i k=1 c k (t)e a k (t)x 2 +b k (t)x .
We can then compute:
dγ dt = i k=1 da k dt x 2 + db k dt x c k e a k x 2 +b k x + dc k dt e a k x 2 +b k x ∈ R
We deduce that the tangent vectors of any smooth submanifold of N i must al o lie in R. In particular this means that our implementation of FunctionRing will be sufficient to represent the tangent vectors of any manifold consisting of finite normal m xtures.

Combining these ideas we obtain the main theoretical result of the paper.

Theorem 6.1 Let θ be a parameterization for a family of probability distributions all of which can be written as a mixture of at most i Gaussians.Let f , a = σ 2 and b be functions in the ring R. In this case one can carry out the direct L 2 projection algorithm for the problem given by equation ( 1) using analytic formulae for all the required integrations.

Although the condition that f , a and b lie in R may seem somewha restrictive, when this condition is not met one could use Taylor expansions to find approximate solutions.

Although the choice of parameterization does not affect the choice of FunctionRing, it does affect the numerical behaviour of the algorithm.In particular if one choos es given later in this paper we parameterize normal mixtures of k Gaussians with a parameterization defined on the whole of R n .We desc ξ i (with 1 ≤ i ≤ k − 1), x 1 , y i (with 2 ≤ i ≤ k) and s i (with 1 ≤ i ≤ k).This gives a total of 3k − 1 parameters.So we can write θ = (ξ 1 , . . ., ξ k−1 , x 1 , y 2 , . . ., y k , s 1 , . . ., s k ) Given a point θ define variables as follows:
λ 1 = logit −1 (ξ 1 ) λ i = logit −1 (ξ i )(1 − λ 1 − λ 2 − . .

− λ i−1 ) (2 ≤ i ≤ k − 1) λ k = 1 − λ 1 − λ 2 − . . . − λ k−1 x i =
x i−1 + e y i (2 ≤ i ≤ k) σ i = e s i
where the logit function sends a probability p ∈ [0, 1] to its log odds, ln(p/1− p).We can now write the density associat d with θ as:
p(x) = k i=1 λ i 1 σ i √ 2π exp(− (x − x i ) 2 2σ 2 i )
We do not claim this is the best possible choice of param terization, but it certainly performs better than some more naïve parameteriations with bounded domains of definition.We will call the direct L 2 projection algorithm onto the normal mixture family given with this projection the L2NM projection filter.


Comparison with the Hellinger exponential (HE) projection al orithm

A similar algorithm is described in [9,10] for projection using the Hellinger metric onto an exponential family.We refer to this as the HE projection filter.

It is worth highlighting the key differences between our algorithm and the exponential projection algorithm described in [9].

• In [9] only the special case of the cubic sensor was considered.It was clear that one could in principle adapt the algorithm to cope with other problems, but there remained symbolic manipulation tha would have to be performed by hand.Our algorithm automates this process by using the FunctionRing abstracti ly, the stochastic term in equation ( 9) simplifies to a term with constant coefficients.This means it can be viewed equally well as either an Itô or Stratonovich SDE.The practical consequence of this is that the HE algorithm can use the Euler-Maruyama scheme rather than the

ratonvoich-Heun sc
eme to solve the resulting stochastic ODE's.Moreover in this case the Euler-Maruyama scheme coincides with the generally more precise Milstein scheme.

• In the case of the cubic senso , the HE algorithm requires one to numerically evaluate integrals such as:
∞ −∞ x n exp(θ 1 + θ 2 x + θ 3 x 2 + θ 4 x 4 )dx
where the θ i are real numbers.Performing such integrals numerically considerably slows the algorithm.In effect one ends up using a rather fine discretization scheme to evaluate the integral and this somewhat offsets the hoped for advantage over a finite difference method.


Numerical Results

In this section we compare the results of using the direct L 2 projection filter onto a mixture of normal distributions with other numerical methods.In par icular we compare it with:

1.A finite difference method using a fine grid which we term the exact filter.Various convergence results are known ( [26] and [27]) for this method.In the simulations shown below we use a grid with 1000 points on the x-axis and 5000 time points.In our simulations we could not visually distinguish the resulting graphs when the grid was refined further justifying us in considering this to be extremely close to the exact result.The precise algorithm used is as described in the section on "Partial Differential Equations Methods" in chapter 8 of Bain and Crisan [5].

2. The extended Kalman filter (EK).This is a somewhat heuristic approach to solving the non-linear filtering problem but which works well so long as one assumes the system is almost linear.It is implemented essentially by linearising all the functions in the problem and then using the exact Kalman filter to solve this linear problem -the details are given in [5].The EK filter is widely used in applications and o provides a standard benchmark.However, it is well known that it can give wildly innaccurate results for non-linear problems so it should be unsurprising to see that it performs badly for most of the examples we consider.

3. The HE projection filter.In fact we have implemented a generalization of the al orithm given in [12] that can cope with filtering problems where b is an aribtrary polynomial, σ is constant and f = 0. Thus we have been able to examine the performance of the exponential projection filter over a slightly wider range of problems t an have previously been considered.

To compare these methods, we have simulated solutions of the equations 1 for various choices of f , σ and b.We have also selected a prior probability distribution p 0 for X and then compared the numerical estimates for the probability distribution p at subsequent times given by the different algorithms.In the examples below we have selected a fixed value for the intial state X 0 rather than drawing at random from the prior distribution.This should have no more impact upon the results than does the choice of seed for the random number generator.

Since each of the approximate methods can only represent certain distributions accurately, we have had to use different prior distributions for each algorithm.To compare the two projection filters we have started with a polynomial exponential distribution for the prior and then found a nearby mixture of normal distributions.This nearby distribution was found using a gradient search a

orithm to minimize
the numerically estimated L 2 norm of the difference of the normal and polynomial exponential distributions.As indicated earlier, polynomial exponential distributions and normal mixtures are qualitatively similar so the prior distributions we use are close for each algorithm.

For the extended Kalman filter, one has to approximate the prior distribution with a single Gaussian.We have done this by moment matching.Inevitably this does not always produce satisfactory results.

or the exact filter,
e have used the same prior as for the L 2 projection filter.


The linear filter

The first test case we have examined is the linear filtering problem.In this case the probability density will be a Gaussian at all times -hence if we project onto the two dimensional family consisting of all Gaussian distributions there should be no loss of information.Thus both projection filters should give exact answers for linear problems.This is indeed the case, and gives some confidence in the correctness of the computer implementations of the various algorithms.


The quadratic sensor

The second test case we have examined is the quadratic sensor.This is problem 1 wit f = 0, σ = c 1 and b(x) = c 2 x 2 for some positive constants c 1 and c 2 .In this problem the non-injectivity of b tends to cause the distribution at any time to be bimodal.To see why, observe that the sensor provides no information about the sign of x, once the state of the system has passed through 0 we expect the probability density to become approximately symmetrical a out the origin.Since we expect the probability density to be bimodal for the quadratic sensor it makes sense to approximate the distribution with a linear combination of two Gaussian distributions.

In Figure 4 we show the probability density as computed by three of the algorithms at 10 different time points for a typical quadratic sensor problem.To reduce clutter we have not plotted the results for the exponential filter.The prior exponential distribution used for this simulation was p(x) = exp(0.25− x 2 + x 3 − 0.25x 4 ).The initial state was X 0 = 0 and Y 0 = 0.

As one can see the probability densities computed using the exact filter and the L2NM filter become visually indistinguishable when the state moves away from the origin.The extended Kalman filter is, as one would expect, completely unable to cope with these bimodal distributions.In this case the extended Kalman filter is simply representing the larger of the two modes.

In Figure 5 we have plotted the L 2 residuals for the different algorithms when applied to the quadratic sensor problem.We define th L 2 residual to be the L 2 norm of the difference between the exact filter distribution and the estimated distribution.As can be s en, the L2NM projection filter outperforms the HE projection filter when applied to the quadratic sensor problem.Notice that the L 2 residuals are initially small for both the HE and the L2NM filter.The superior performance of the L2NM projection filter in th

case stems from
he fact that one can more accurately represent the distributions that occur using the normal mixture family than using the polynomial exponential family.

If preferred one could define a similar notion of residual using the Hellinger metric.The results would be qualitatively similar.

One interesting feature of Figure 5 is that the error remains bounded in size when one might expect the error to accumulate over time.This suggests that the arrival of new measurements is gradually correcting for the errors introduced by the approximation.


The cubic sensor

A third test case we have considered is the general cubic sensor.In this proble one has f = 0, σ = c 1 for some constant c 1 and b is some cubic function.

The case when b is a multiple of x 3 is called the cubic sensor and was used as the test case for the exponential projection filter using the Hellinger metric considered in [12].It is of interest because it is the simplest case where b is injective but where it is known that the problem cannot be reduced to a finite dimensional stochastic differential equation [19].It is known from earlier work that the exponential filter gives excellent numerical results for the cubic sensor.

Our new implementations allow us to examine the general cubic sensor.In Figure 6, we have plotted example probability densities over time for the problem with f = 0, σ = 1 and b = x 3 − x.Wi h two turning points for b this problem is very far from linear.As can be seen in Figure 6 the L2NM projection remains close to the exact distribution throughout.A mixture of only two Gaussians is enough to approximate quite a var ety of differently shaped distributions with perhaps surprising accuracy.As expected, the extended Kalman filter gives poor results until the state moves to a region where b is injective.The r sults of the exponential filter have not been plotted in Figure 6 to reduce clutter.It gave similar results to the L2NM filter.

The prior polynomial exponential distribution used for this simulation was p(x) = exp(0.5x 2 − 0.25x 4 ).The initial state was X 0 = 0, which is one of the modes of prior distribution.The inital value for Y 0 was taken to be 0.

One new phenomen n that occurs when considering the cubic sensor is that the algorithm sometimes abruptly fails.This is true for both the L2NM projection filter and the HE projection filter.

To show the behaviour over time more clearly, in Figure 7 we have shown a plot of the mean and standard deviation as estimated by the L2NM projection filter against the actual mean and standard deviation.We have also indicated the true state of the system.The mean for th L2MN filter drops to 0 at approximately time 7.It is at this point that the algorithm has failed.

What has happened is that as the state has moved to a region where the sensor is reasonably close to being linear, the probability distribution has tended to a single normal distribution.Such a distribution lies on the boundary of the family consisting of a mixture of two normal distributions.As we approach the boundary, h ij ceases to be invertible causing the failure of the algorithm.Analogous phenomena occur for t filter.

The result of running numerous simulations suggests that the HE filter is rather less robust than the L2NM projection filter.The typical behaviour is that the exponential filter maintains a very low residual right up until the The L 2 residuals of the L2MN method are rather large between times 6 and 7 but note that the accuracy of the estimates for the mean and standard deviation in Figure 7 remain reasonable throughout this time.To understand this note that for two normal distributions with means a distance x apart, the L 2 distance between the distributions increases as the standard deviations of the distributions drop.Thus the increase in L 2 residuals between times 6 and 7 is to a large extent due to the drop in standard deviation between these times.As a result, one may feel that the L 2 residual doesn't capture precisely what it means for an approximation to be "good".In the next section we will s ow how to measure residuals in a way that corresponds more closely to the intuitive idea of them having visually similar distribution functions.In practice one's definition of a good approximation will depend upon the Although one might argue that the filter is in fact behaving reasonably well between times 6 and 7 it does ultimately fail.There is an obvious fix for failures like this.When the current point is sufficiently close to the boundary of the manifold, simply approximate the distribution with an element of the boundary.In other words, approximate the distribution using a mixture of fewer Gaussians.Since this means moving to a lower dimensional family of distributions, the numerical implementation will be more efficient on the boundary.This will provide a temporary fix the failure of the algorithm, but it raises another problem: as the state moves back

nto a region where the problem is
highly non linear, how can one decide how to leave the boundary and start adding additional the mixture?We hope to address this question in a future paper.
0 0.5 1 1.5 -2 0 2 4 6 t=0 L2MN Exact EK X t=1 X t=2 X t=3 X t=4 0 0.5 1 1.5 -2 0 2 4 6 t=5 X t=6 X t=7 X t=8 X t=9

Comparison with Particle Method Particle methods approximate the probability density p using discrete measures of the form:
i a i (t)δ v i (t)
These measures are generated using a Monte Carlo method.The measure can be thought of as the empiri al distributions associated with randomly located particles at position v i (t) and of stochastic mass a i (t).

Particle methods are currently some of the most effective numerical methods for solving the filtering problem.See [5] and the references therein for details of specific particle methods and convergence results.

The first issue in comparing projection methods with particle methods is that, as a linear combination of Dirac masses, one can only expect a particle method to converge weakly to the exact solution.In particular the L 2 metric and the Hellinger metric are both inappropriate measures of the residual between the exact solution and a particle approxi ellinger distance will always take the value √ 2.

To combat this issue, we will measure residuals using the Lévy metric.If p and q are two probability measures on R and P and Q are the associated cumulative distribution functions then the Lévy metric is defined by:
d L (p, q) = inf{ : P (x − ) − ≤ Q(x) ≤ (x + ) + ∀x}
This can be interpreted geometrically as the size of the largest square with sides parallel to the coordinate axes that can be inserted between the completed graphs of the cumulative distribution functions (the completed graph of the distribution function is simply the graph of the distribution function with vertical line segments added at discontinuities).

The Lévy metric can be seen as a special case of the Lévy-Prokhorov metric.This can be used to measu e the distance between measures on a general metric space.For Polish spaces, the Lévy-Prokhorov metric metrises the weak convergence of probability measures [4].Thus the Lévy metric provides a reasonable measure of the residual of a particle approximation.We will call residuals measured in this way Lévy residuals.

A second issue in comparing projection methods with particle methods is dec ding how many particles to use for the comparison.A natural choice is to compare a projection method onto an m-dimensional manifold with a particle method that approximates the distribution using (m + 1)/2 particles.In other words, equate the dimension of the families of distributions used for the approximation.

A third issue is deciding which particle method to choose for the comparison from the many algorithms that can be found in the literature.We can w rk around this issue by calculating the best possible approximation to the exact distribution that can be made using (m + 1)/2 Dirac masses.This approach will substantially underestimate the Lévy residual of a particle method: being Monte Carlo methods, large numbers of particles would be requi ed in practice.

In Figure 9 we have plotted bounds on the Lévy residuals for the two projection methods for the quadratic sensor.Since mixtures of two normal distributions lie in a 5 dimensional family, we have compared these residuals with the best possible Lévy residual for a mixture of three Dirac masses.

To compute the Lévy residual between two functions we have approximated first approximated the cumulative distribution functions using step functions.We have used the same grid for these steps as we used to compute our "exact" filter.We have then used a brute force approach to compute a b und on size of the largest square that can be placed between these step functions.Thus if we have used a grid with n points to discretize the xaxis, we will need to make n 2 comparisons to estimate the Lévy residual.More efficient algorithms are possible, but t is approach is sufficient for our purposes.

The maximum accuracy of the computation of the Lévy metric is constrained by the grid size used for our "exact" filter.Since the grid size in the x direction for our "exact" filter is 0.01, our estimates for the projection residuals are bounded below by 0.02.

The computation of the minimum residual for a particle filter is a little more complex.Let minEpsilon(F, n denote the minimum Lévy distance between a distribution with cumulative distribution F and a distribution of n particles.Let minN(F, ) denote the minimum number of particles required to approximate F with a residual of less than .If we can compute minN we can use a line search to compute minEspilon.

To compute minN(F, ) for an increasing step function F with F (−∞) = 0 and F (∞) = 1, one needs to find the minimum number of steps in a similar increasing step function G that is never further than away from F in the L ∞ metric.One constructs candidate step functions G by starting with G(−∞) = 0 and then moving along the x-axis adding in additional steps as required to remain within a distance .An optimal G is found by addin in steps as late as possible and, when adding a new step, making it as high as possible.

In this way we can compute minN and minEpsilon for step functions F .We can then compute bounds on these values fo

a given distribution by approxim
ting its cumulative density function with a step function.

As can be seen, the exponential and mixture projection filters have similar accuracy as measured by the Lévy residual and it is impossible to match this accuracy using a model containing only 3 particles.


Conclusions and Future Research

Projection onto a family of normal mixtures using the L 2 metric allows one to ap roximate the solutions of the non-linear filtering problem with surprising accuracy using only a small number of component distributions.In this regard it behaves in a very similar fashion to the projection onto an exponential family using the Hellinger metric that has been considered previously.

The L2NM projection filter has one important advantage over the HE projection filter, for problems wi h polynomial coefficients all required integrals can be calculated analytically.Problems with more general coefficients can be addressed using Taylor series.One expects this to translate into a better performing algorithm -particularly if the approach is extended to higher dimensional problems.

We tested both filters against the optimal filter in simpl but interesting systems, and we provided a metric to compare the performance of each filter with the optimal one.We also tested both filters against a pa ticle method, showing that with the same number of parameters the L2NM filter outperforms the best possible particle method in Levy metric.

We designed a software structure and populated it with models that make the L2NM filter quite appealing from a numerical and computational point of view.

Areas of future research that we hope to address include: the relat