An Efficient Meshfreee Implicit Filter for Nonlinear Filtering Problems

In this paper, we propose a meshfree approximation method for the implicit filter developed in [2], which is a novel numerical algorithm for nonlinear filtering problems. The implicit filter approximates conditional distributions in the optimal filter over a deterministic state space grid and is developed from samples of the current state obtained by solving the state equation implicitly. The purpose of the meshfree approximation is to improve the efficiency of the implicit filter in moderately high-dimensional problems. The construction of the algorithm includes generation of random state space points and a meshfree interpolation method. Numerical experiments show the effectiveness and efficiency of our algorithm.

where Φ(·) is a test function, and Z k denotes the space of all Y 1:k -measurable and square integrable random variables.
Tremendous efforts have been made to solve nonlinear filtering problems in the last few decades. Two of the well-known Bayesian filters are extended Kalman filters (EKFs) [3,10,16,17,19], and particle filters [5,7,8,13]. The key ingredient of the EKFs is the linearization of both f and g in (1.1), so that the standard Kalman filter can be applied directly. Thus, if the nonlinearity of the state and the observation systems is not severe, then the EKFs can provide efficient and reasonable inferences about the state, otherwise, the performance of the EKFs can be very poor. For particle filters, the central theme is to approximate the desired posterior PDF of the state by the empirical distribution of a set of adaptively selected random samples (referred to as "particles"). The particle filter method is essentially a sequential Monte Carlo approach, which requires no assumption on the linearity of the underlying system. As such, with sufficiently large number of particles, it is capable of providing an accurate approximation of the posterior PDF for a highly nonlinear filtering problems. However, there are some fundamental issues concerning the efficiency and robustness of particle filters [9]. For example, since the empirical PDF is constructed based on particles with equal weights after resampling, the particle filter still needs a lot of samples in order to accurately approximate the target distribution.
To overcome such a disadvantage, the authors proposed a new nonlinear filter named "implicit filter" [2] . This approach adopts the framework of Bayesian filtering, which has two stages at each time step, i.e., prediction and update. At the prediction stage, we estimate the prior PDF p(X k |Y 1:k−1 ) of the future state X k given the current available observation information Y 1:k−1 ; at the update stage, we update the prior PDF by assimilating the newly received data Y k to obtain the estimate of the posterior PDF p(X k |Y 1:k ). The implicit filter is distinguished from the particle filters by the use of interpolatory approximations to the prior and posterior PDFs. Specifically, in the particle filter, p(X k |Y 1:k−1 ) is approximated by explicitly propagating the samples of the current state X k−1 |Y 1:k−1 through the nonlinear state equation X k = f k (X k−1 , w k−1 ), and constructing the empirical PDF of X k |Y 1:k−1 . In the implicit filter, the interpolation of p(X k |Y 1:k−1 ) requires its function values at a set of grid points of the future state X k . Under the condition that X k = x ∈ R d , we solve implicitly the state equation x = f k (X k−1 , w k−1 ) given a set of Monte Carlo samples of w k−1 , so that the value of p(X k = x|Y 1:k−1 ), at the grid point of x, can be estimated by averaging the function values of p(X k−1 |Y 1:k−1 ) at all the solutions of the state equation. As an implicit scheme, the implicit filter has a stabilizing effect which provides more accurate numerical approximations to the solution of the nonlinear filtering problem than the particle filter method [2].
The main challenge of the implicit filter method is that the conditional PDF of the nonlinear filtering solution is estimated at grid points. As such the method suffers the so called "the curse of dimensionality" when the dimension of the state variable is high. In addition, the efficiency of the method may be significantly reduced when the domain of the PDF is unbounded. In this paper, we propose to construct a meshfree implicit filter algorithm to alleviate the aforementioned challenges. Motivated by the particle filter method, we first generate a set of random particles and propagate these particles through the system model and use these particles to replace the grid points in the state space. After that we generate other necessary points through the Shepard's method which constructs the interpolant by the weighted average of the values on state points [12]. In order to prevent particle degeneracy in the generation of random state points, we introduce a resample step in the particle propagation. In addition we choose state points according to the system state, which make them adaptively located in the high probability region of the PDF of state. In this way, we solve the nonlinear filtering problem in a relatively small region in the state space at each time step and approximate the solution on a set of meshfree state points distributed adaptively to the desired PDF of the state. Furthermore, since we approximate the PDF as a function on each state point, instead of using state points themselves to describe the empirical distribution, the implicit filter algorithm requires much fewer points than the particle filter method to depict the PDF of the state.
The rest of this paper is organized as follows. In §2, we introduce the mathematical framework of the Bayesian optimal filter. In §3, we construct meshfree implicit algorithm. In §4, we demonstrate the efficiency and accuracy of our algorithm through numerical experiments. Finally, §5 contains conclusions and directions for the future research.
2. Bayesian optimal filter. For m, n ∈ N + , let X m:n and Y m:n denote the σ fields generated by {X m , X m+1 , . . . , X n } and {Y m , Y m+1 , . . . , Y n }, respectively. For k = N + , we use x k to represent a realization of the random variable X k , and define for notational simplicity. It is easy to see that the dynamical model in (1.2) is Markovian in the sense that We also know that the measurements are conditionally independent given x k , i.e., The Bayesian optimal filter constructs the conditional distribution p(x k |Y 1:k ) recursively in two stages: prediction stage and update stage.
For k = 1, 2, · · · , assume that p(x k−1 |Y 1:k−1 ) is given. In the prediction stage p(x k |Y 1:k−1 ) is evaluated through the Chapman-Kolmogorov formula: (2.1) In the update stage, the prior PDF obtained in (2.1) is used to obtain the posterior PDF p(x k |Y 1:k ) via the Bayes' formula: 3. The meshfree implicit filter. In this section, we construct the meshfree implicit filter algorithm. The algorithm is based the implicit filter algorithm on grid points [2]. The implicit filter algorithm introduced in [2] is developed from the general framework of the Bayesian optimal filter discussed above, in which the primary computational challenge is the numerical approximation of the term p(x k |x k−1 ) in (2.1).
3.1. The prediction stage. For k = 1, 2, · · · , the goal of this stage is to approximate the prior distribution p(x k |Y 1:k−1 ) of the state X k , given the posterior distribution p(x k−1 |Y 1:k−1 ) of the state X k−1 . Due to the the fact that the prior PDF p(x k |Y 1:k−1 ) derived in identity (2.1) can be rewritten as represents the expectation with respect to the white noise w k−1 , and the PDF p(x k |x k−1 , w k−1 ) is 2) can be viewed as a generalization of the Dirac delta function in the space R d × R d × R r , where the mass is located according to the state equation Note that the estimation of (3.1) requires an approximation to the expectation of the white noise w k−1 , and define an approximation to p(x k |x k−1 , w k−1 ) as To construct an interpolation of p(x k |Y 1:k−1 ), the next step is to approximate is the approximation error. Then, by further fixing w k−1 = w j k−1 , the location of the mass can be obtained by implicitly solving the state equation which is the reason we named the approach the implicit filter. Now substituting x i k into (3.4), and using the same sample set {w j k−1 } M j=1 as above to approximate the integral on the right hand side of (3.4), we obtain , we obtain the following iterative numerical scheme for constructing an approximation, denoted by (x i k |Y 1:k−1 ), of the prior PDF p(x i k |Y 1:k−1 ), i.e., In our previous work [2], the subsets H k , for k = 0, 1, . . ., were defined by a full tensor product mesh, denoted by It is simple to implement but has several significant disadvantages. First, at each time step, one needs to approximate the prior PDF p(x k |Y 1:k−1 ) at a total of N (1) × · · · × N (d) grid points which grows exponentially as the dimension d increases. This is also known as "the curse of dimensionality". On the other hand, since the construction of M is not informed by the target PDF, the domain [a 1 , b 1 ] × · · · × [a d , b d ] needs to be defined sufficiently large, so as to capture the statistically significant region of the PDF. This may lead to a great waste of computation effort in the low probability region of p(x k |Y 1:k−1 ).
To alleviate such disadvantages, we propose to develop a distribution-informed meshfree interpolation approach to efficiently approximate the prior PDF. The central idea of the generation of random points for the state variable is to build a set of points, denoted by H k , according to the state distribution. To begin with, we generate of N random samples from the initial PDF p 0 of the initial state: If the initial PDF p 0 is close to the true state distribution, it's obvious that our random state points are more concentrated near the target state. For k = 1, 2, · · · , K, we propagate points and approximate the conditional PDF p(x k |Y 1:k−1 ) on H k with the scheme given by (3.8). In this way, the random points in H k move according to the state model. As opposed to particle filter methods, which use the number of particles to represent empirical distributions and require a large number of particles to follow the state distribution, in the implicit filter method we provide an approximation of the value of the PDF at each state point. Therefore, much fewer points are needed to describe the state PDF and the random state points are not necessary to accurately follow the state distribution.
3.2. The update stage. By incorporating the new data Y k , we update the prior PDF p(x k |Y 1:k−1 ) at each grid point x i k , using the Bayesian formula, to obtain where (x i k |Y 1:k−1 ) is given in (3.8), C k is the normalization factor, and R i k|k : is the approximation error. By neglecting the error term R i k|k in (3.10), we obtain the following iterative numerical scheme for the update stage on H k , i.e., where (x i k |Y 1:k ) is desired the approximation of the posterior PDF p(x i k |Y 1:k ). Next, we use interpolation methods to construct the approximation (x k |Y 1:k ) of is the set of basis functions. Since the state points in H k are generated randomly in the meshfree framework, standard polynomial interpolation [2] is unstable due to the uncontrollable Lebesgue constant. Instead, we propose to use the Shepard's method [12], which is an efficient meshfree interpolation technique, to construct the interpolant (x k |Y 1:k ). The basic idea of the Shepard's method is to use the weighted average of { (x i k |Y 1:k )} x i k ∈H k in the interpolating approximation. Specifically, for a given point x k ∈ R d , we re-order the points in H k by the distances to x k to get a sequence {x where · is the Euclidean norm in R d . Then, for a pre-chosen integer L we use the where the weight h l (x k ) is defined by is the error of the Shepard's interpolation. We assume that p(x k |Y 1:k ) has bounded first order derivative. For each pair p(x k |Y 1:k ) and p(x k . It is reasonable to assume that in high probability region of the derivative p (z|Y 1:k ) is large. It's worth pointing out that the random state points generated in this algorithm are concentrated in the high probability region. Thus, if x k lies in the high probability region, the distance is small, which balances the error brought by the large derivative. On the other hand, if x k lies in the low probability region, although the distance x k − x (l) k is relatively large, the approximation error (3.14) is still small due to the small value of the derivative p (z|Y 1:k ).

3.3.
Resampling. Similar to the particle filter method, the above random state points generation suffers from the degeneracy problem for long term simulations, especially for high-dimensional problems. After several time steps, the probability density tends to concentrate on a few points which dramatically reduces the number of effective sample points in H k .
In this work, we propose an occasional resampling procedure to address these problems and rejuvenate the random points cloud. At the time step k − 1, the resampling procedure takes place after we obtain (x k−1 |Y 1:k−1 ), in order to remove the degenerated points in H k−1 using the information provided by (x k−1 |Y 1:k−1 ). Specifically, the first step is to develop a degeneracy metric to determine the necessity of doing resampling. To this end, we define the following degenerated subset S k−1 ⊂ H k−1 , where ε > 0 is a user-defined threshold. We also define J (S k−1 ) := {i = 1, . . . , N |x i k−1 ∈ S k−1 } to be the index set of S k−1 . Then, the degeneracy of H k−1 can be measured by the ratio #(S k−1 )/#(H k−1 ) ∈ [0, 1], where #(·) denotes the number of points in a set. If the ratio is smaller than a threshold τ ∈ [0, 1], then we will skip the resampling step and propagate H k−1 to get H k ; otherwise, the set H k−1 is considered degenerated, and the resampling procedure is needed.
In resampling, instead of propagating H k−1 to H k , we aim at constructing an and propagate H k− 1 2 through the state model (1.2) to obtain H k . According to the definition of S k−1 in (3.15), we consider the state points in H k−1 \S k−1 are in the statistically significant region of (x k−1 |Y 1:k−1 ), so that we first put those points in H k− 1 2 , i.e., For the state points in S k−1 , we replace them by generating new samples from (x k−1 |Y 1:k−1 ) using the importance sampling [6], i.e., As a result, the resampling procedure helps us remove the state points with low probabilities, and makes the state point set H k concentrated in the high probability region of the posterior PDF (x k−1 |Y 1:k−1 ) at each time step.
3.4. Summary of the algorithm. Finally, we summarize the entire meshfree implicit filter algorithm introduced in §3.1- §3.3 in Algorithm 1 below.

Numerical experiments.
In this section, we present two numerical examples to examine the performance of our meshfree implicit filter method. In Example 1, we use a two dimensional nonlinear filtering problem to show the distributions of the random points H k . In Example 2, we solve a three dimensional bearing-only tracking problem, which is a six dimensional nonlinear filtering problem. For this higher dimensional problem, we compare the accuracy and efficiency of our meshfree implicit filter method with the extended Kalman filter and the particle filter. Example 1. In this example, we consider the two dimensional noise perturbed tumoral growth model [21] where W t is a two dimensional standard Brownian motion and σ = (0.01, 0.01) T . The state process X t = (X 1 t , X 2 t ) T is a two dimensional vector, F (X t ) := (f 1 (X t ), f 2 (X t )) T is defined as .
Here, f 1 models the Gompertzian growth rate of the tumor and f 2 gives the degree of vascularization of the tumor which is also called " angiogenic capacity". To approximate the state variables, we discretize the dynamic system (4.1) in time and obtain a discrete state model Here, ω k is a two dimensional zero mean Gaussian white noise process with covariance Q = I∆, where I is the 2 × 2 identity matrix and ∆ is the time partition stepsize. The measurement of the state model is given by where v k is a two dimensional zero mean Gaussian white noise process with covariance Λ = I∆, I is a 2 × 2 identity matrix and R = (0.1, 0.1) T . In the numerical experiment, we use uniform time partition with stepsize ∆ = 0.2 and simulate the state process for K = 40 with initial state X 0 = (0.8, 0.3) T and parameters α 1 = 1, α 2 = 0.2, α 3 = 0.2. At time step k = 0, we initialize the prior PDF p 0 by N (X 0 , Σ), whereX 0 = (0.78, 0.32) T and In Figure 4.1, we plot 1500 random samples generated from the initial PDF p 0 , which are our initial random points H 0 . To present the accuracy of the algorithm, we show the simulation of the tumoral growth states in Figure 4.3. The black curves are the true X 1 and X 2 coordinate values of the tumoral growth states, respectively. The blue curves show the simulated states obtained by using the meshfree implicit filter method.
Example 2. In this example, we study a six dimensional target tracking problem. In Figure 4.4, the target, denoted by the red line, moves in the three dimensional space and two platforms on the ground, denoted by pentagons, take angular observations of the moving target.
The state process X k = (X 1 k , X 2 k , X 3 k , X 4 k , X 5 k , X 6 k ) T is described by the following dynamic model where (X 1 , X 2 , X 3 ) describes the position of the moving target which is controlled by parameters (X 4 , X 5 , X 6 ). The system noise ω k = (ω 1 k , ω 2 k , ω 3 k , ω 4 k , ω 5 k , ω 6 k ) T is a zero mean Gaussian white noise process with covariance Q . = I∆, I is the 6 × 6 identity matrix and ∆ is a given time period, σ = (0.1, 0.1, 0.1, 0.01, 0.01, 0.01) T is a constant vector and f is given by where v k is a 4 dimensional zero mean Gaussian white noise process with covariance Λ = I∆, I is a 4 × 4 identity matrix, R = (0.6, 0.6, 0.6, 0.6) T , (a 1 , b 1 ) and (a 2 , b 2 ) are locations of two observers. We choose ∆ = 0.3, α = 3, v 1 = v 2 = v 3 = 0.05 . Also, we assume that platforms are located at (a 1 , b 1 ) = (16, 6), (a 2 , b 2 ) = (8, 15) and the initial sate is given by a Gaussian N (X 0 , Σ) where X 0 = (2, 2, 1, 0.4, 0.4, 0) T and The target will be observed over the time period 0 ≤ t ≤ 15. In the numerical experiments, we compare the performance of our meshfree implicit filter with the extended Kalman filter and the particle filter. In particular, we compare the estimated mean values of the states process along each dimension in Figure 4.5. In the particle filter method, we choose 15, 000 particles. In the meshfree implicit filter method, we choose the number of state points to be N = 4, 000 and the number of random samples in the implicit filter Monte Carlo simulation to be M = 6. The black curves in Figure 4.5 show the real states process along each direction, the green curves give the estimated means obtained by the extended Kalman filter method, the red curves give the estimated means obtained by the particle filter method, and the blue curves give the estimated means obtained by the meshfree implicit filter. We also plot the L 2 error err k corresponding to all three methods in figure 4.6. As we can see from figure 4.5 and 4.6, the implicit filter and the particle filter are much more accurate than the extended Kalman filter and the implicit filter is the most accurate approximation in this experiment.
To further compare the efficiency between the meshfree implicit filter and the particle filter, we repeat the above experiment over 50 realizations and show the average CPU time and the corresponding global root mean square error err G defined by where err k (j) is the L 2 error of the j-realization at time step k. In table 4.1, we can see that with 15, 000 particles, the CPU time of the particle filter method is comparable to that of the implicit filter with 4, 000 random state points, but the global RMSE of the particle filter is more than doubled the RMSE of the implicit filter. With 25, 000 particles, the particle filter method achieves an accuracy comparable to the implicit filter, but at a significantly higher cost.

Conclusions.
In this work, we proposed an efficient meshfree implicit filter algorithm by evaluating the conditional PDF on meshfree points in the state space. These meshfree points are chosen adaptively according to the system state evolution. We also apply Shepard's method as the meshfree interpolation method to compute interplants with random state points. In order to address the degeneracy of the random points, we use importance sampling method to construct a resample step. Numerical examples demonstrate the effectiveness and efficiency of our algorithm. In   the future, we plan to perform a rigorous numerical analysis for the meshfree implicit filter algorithm.