An Efﬁcient Particle Filtering Algorithm based on Ensemble Kalman Filter Proposal Density

—This paper describes a novel Particle Filter al- gorithm that is both computationally efﬁcient and accurate. The proposed Particle Filter addresses the issue of degeneracy by using the Ensemble Kalman Filter as the proposal density function. The suggested methodology combines the strength of the Ensemble Kalman Filter’s Linear Gaussian update with the Monte Carlo technique. The technique considerably decreases processing effort while also resolving the problem or degeneracy inherent in traditional Sequential Importance Algorithms. The problem of Bayesian state estimation in a non-linear/non-Gaussian framework is considered. The proposed algorithms is studied in comparison with the other available proposal densities based on Gaussian approximation. It is demonstrated that the suggested method does not require explicit computation of error covariance for each sample, signiﬁcantly reducing computational effort. The proposed algorithm’s suitability and effectiveness have been thoroughly validated in the simulation with benchmark non- linear estimation problem for various scenarios.


I. INTRODUCTION
B AYESIAN estimation is a method to combine the evidence or likelihood with prior knowledge to get an update or posterior probability using Bayes' theorem. Once the posterior probability is available, the "measurement of the central tendency" is used in most cases to estimate the state (variable or parameter). Bayesian state estimation is an important technique in Bayesian statistics and has found applications in areas ranging from science, engineering, medicine, robotics and so on.
In almost all the complex applications of the real world, posterior probability density is never available in closed form, therefore approximate techniques are the only alternatives. Some algorithms are approximations of the process model while some attempt to approximate the state in probability space. In complex applications like meteorology or machine learning, the state and measurement dimensions are enormous, and state evolution model are highly non-linear. Thus, nonlinearity, non-Gaussianity and high dimensionality has increasingly become an inescapable elements to take into account in the state estimation to precisely model the dynamics of the physical systems. The state estimation for such complex scenarios entirely depends upon the accuracy with which the model can be approximated. The solution is straightforward for a simplified case where the system evolves linearly and the posterior density is exactly represented by mean and variance (for a Gaussian posterior). An exact analytical solution can be found for a sequentially evolving linear Gaussian dynamical system, which is known as the Kalman Filter.
Over the last few decades many different schemes were developed to solve the aforementioned issues. One of the most popular algorithm is Extended Kalman Filter (EKF) which is based on the Taylor series expansion of state evolution function. Other solutions are Gaussian Sum approximation and Grid based filters. However these schemes were still not sufficient to solve complex nonlinear problems. Recent development of high-speed digital computing has opened a new arena of the class of filters known as particle filters or sequential Monte Carlo methods using approximation in probability space [1], [2], [3], [4]. These algorithms have the Bayes' theorem at its heart, which has proven to be the most natural and logical way to combine prior knowledge to make a correct estimate of the posterior. Particle Filtering algorithms represent posterior density by a set of samples. They are drawn from the posterior density function and are represented in combination with the associated weights. However, it may not be always straightforward to draw from the posterior density, hence a proposal density is used to draw the samples. The sample weights are updated in a recursive manner according the most recently obtained observation data.
Despite being a powerful tool, Particle Filtering still has some serious issues. When proposal density is not properly selected, after a few iterations, all but one sample carries negligible weights. This happens due to lack of overlap between the likelihood and proposal density from which the samples are drawn. The phenomenon in particle filter estimation is also known as "degeneracy". Conventional particle filtering algorithm, also known as "Bootstrap" filter, "Condensation", uses transition prior as the proposal density [1]. This technique is severely effected by degeneracy, since the samples are drawn without the considerations of observation, and there is a very high probability that the many of the sample may fall on the low likelihood are zone. Our goal in this article is to focus on the issue of degeneracy in Particle Filtering and compare different type of proposal densities.
This study provides, first, an overview of the various methods of choosing the importance density functions to draw the samples, Second, we suggest the Ensemble Kalman Filter (EnKF) as a suitable proposal density function. The Ensemble Kalman Filter (EnKF) first appeared in Evensen (1994b) [5] as a new data assimilation technique over the computationally demanding Extended Kalman Filter (EKF). It also appeared in many other articles [6], [7], [8], [9]. The EnKF is an approximate computational technique which represents the probability distribution of system state with the help of a set of samples (ensemble) drawn from the posterior distribution. The samples drawn are restricted by the assumption of Gaussian posterior distribution parametrized by "sample mean" and "sample variance". The EnKF has become quite popular in Geophysics and is very convenient to implement for real-time applications. In this paper, we present the proposed algorithm's formulation as well as a comparison of performance with various proposal densities. In terms of sample size and critical characteristics such as convergence, estimation error, and computation time, we draw some useful conclusions for EnKF-based PF.
The rest of the paper is organised as follows. Section II discusses Particle Filtering, as well as different implementation strategies and non-linear state-space estimation difficulties. Section III covers over Ensemble Kalman Filter formulation. Section IV provided the description of the proposed EnKF based Particle Filter methodology. In Section V, we compare the performance of the proposed methodology to that of current approaches in terms of error, convergence, and computing time, and in Section VI we conclude.

II. PARTICLE FILTER
Particle Filters are a class of recursive filter that approximates the posterior density based on Monte Carlo integration [10]. The main idea behind Particle Filtering is to recursively estimate the approximate posterior density in a non-linear environment by using a set of samples and their associated importance weights. Once a new measurement is available, each iteration entails updating the samples and their importance weights. Consider a discrete-time non-linear state evolution model, and measurements as where F : R n → R n is the sate transition function. G : R n → R p is the state observation function. x t , η t ∈ R n×1 are n-dimensional system state and process disturbance vectors respectively at a discrete time instant t. Similarly y t , t ∈ R p×1 are p-dimensional measurement and noise sequence respectively. It is assumed that η t and t are zero mean white noise sequences, mutually and serially independent to each other with known PDFs.
For simplicity, we assume that the state evolution is a 1st order Markovian process. Assuming we have the posterior density p(x t−1 |Y t−1 ) at the time instant t − 1 where Y t−1 = {y j : j = 1, 2, . . . , t − 1} represent a set of measurement up to time instant t − 1. We are also generally aware of the model for evaluating likelihood density function p(y t |x t ). The prior density p(x t |Y t−1 ) at time t can be evaluated using the following integral where the probabilistic model p(x t |x t−1 ) for the state transition equation using Eq. 1 is described by where δ(·) is the Dirac-delta function. Once the latest measurement y t is available, the posterior is evaluated by the Bayesian update where denominator is normalising constant Following recursion can be built for the sequential case Further, it is assumed that observation y t at any time is conditionally independent of its past values. The posterior density becomes This leads to a simplified recursive expression using only the last observation The particle filter approximates the posterior density function using a set of finite samples {x (i) t ; i = 1, 2 . . . , N } with their discrete probability masses ω (i) t called importance weights, such that The weights ω (i) t are chosen as per the importance density principal. Since it is not always straightforward to draw samples from p(x t |y t ), the sample are drawn from the importance density, or proposal density π(x t |y t ) from which they can be easily drawn. The required conditions for the importance density π(x t |y t ) is that it must contain the support of target density p(x t |y t ). The weights can be defined as [11] Further factorising of eqn. (11) along with the posterior density recursion established in eqn. (7) lead to a recursive expression for the importance weights [11].
Thus, if we can generate the samples from proposal density π x t−1 , y t and evaluate the weights using (12), we can recursively estimate the approximate posterior density given by (10). The derivation for eqn. (12) is provided in The next most important thing is to choose the right proposal density. It is an established fact that if there is lack of support points, variance of the weights is bound to increases and a phenomenon known as degeneracy occurs. Consequently, all the computation effort is spent to propagate only one sample that decreases the overall estimation accuracy as well is a waste of computational efforts. There are infinite choice of proposal density. An optimal proposal density will keep the variance of sample weights to minimum, and hence all the weights will be almost equal. It is shown that the optimal proposal density is p(x t |x t−1 , y t ) [12], however, it is not always easy to draw samples from this as stated earlier.
An important step to tackle degeneracy, is re-sampling of the samples according to their importance weights. As the number of significant samples drop considerably, re-sampling is performed to increase the occurrence of samples with more weights. The basic idea is to avoid samples that does not contribute to the estimation accuracy and to utilize computation power more effectively. The pseudo-code for the basic Sequential Importance Re-sampling (SIR) algorithm is described in Algorithm 1.

A. Bootstrap Particle Filter
An early solution was given by Gorden et. al. (1993) [1], which was to draw the samples from transitional prior p(x t |x t−1 ). The technique is known as "bootstrap" or "Condensation" algorithm. With transitional prior as the proposal, the weight recursion simply becomes The above technique makes computation quite simpler, and we only need to evaluate the likelihood of samples after they are propagate in time through the stochastic process model, to generate new weights.
Although the technique is quite simple to implement, its success heavily depends upon the amount of overlap between the prior and the likelihood. If the prior samples doesn't fall close to the likelihood area, most of the sample weights will soon have almost no importance, and much of computation will be wasted, causing severe degeneracy.

B. Observation Driven Proposal Densities
To reduce the degeneracy; or in other words, to increase the amount of overlap between the samples and likelihood, it is required to generate the samples from a proposal density that incorporates the most recent measurement. The conventional bootstrap Particle Filter moves blindly without using the latest measurement available to it and hence is prone to severe degeneracy, and sensitive to outliers.
Obviously, the most important point in the design of particle filters is to incorporate the most recent measurement into the proposal. A variant of the SIR algorithm was introduced by Pitt and Shephard (1994) [13] known as Auxiliary Particle Filter (APF) that attempts to draw the samples from a proposal density p(x t |y t ). The APF algorithm is derived with the help of an auxiliary variable -hence its name. The importance density considered in APF is π(x t , j|y t ) from whih the sample pair (x (i) t , j (i) ) are drawn, where the auxiliary variable j is used as an index. The technique can be considered as resampling from one step previous time instant t − 1, taking the latest measurement y t into consideration. The samples are propagated to evaluate the mean µ k is some value that is able to characterize the prior density, although other indicators are also possible. The likelihood p(y t |µ (i) t ) is evaluated to calculate the weights that are used for re-sampling. The resampled particles are then processed using the generic SIR algorithm, and the weights are updated. The weights produced by this method should be much more evenly distributed, and the estimated state should be close to the true state. However, APF has the drawback of performing well only when the process noise is low. Because a single point cannot characterise the posterior density p(x t |x (i) t−1 ) in the presence of substantial process noise, the APF will yield poor results and be prone to degeneracy due to the poor approximation.
Another class of likelihood-based particle filters draw samples directly from proposal densities based on Kalman filterbased estimators, which are specifically intended for estimating the state in a non-linear environment. Two major variants of Kalman Filter that deal with the non-linear process and measurement models are Extended Kalman Filter (EKF) and Unscented Kalman Filter. There are similar particle filter hybridisation schemes has been in the literature for various applications [14], [15]. Particle Filters variants using these approaches are known as EKF Particle Filter (EPF) and UKF Particle Filter (UPF) [16] respectively. Since these filters recursively update the state using likelihood and provide an estimate of the state, it is quite simple to incorporate them in the SIR algorithm and draw samples. Although they are Gaussian approximation to the posterior p(x t |x t−1 , y t ), they tend to successfully reduce variance among the importance weights, and keep the number of effective samples large enough.
The EKF is most widely used estimator for non-linear systems and is the result of a linearisation of the state evolution model using the Taylor series expansion. However the major drawback of EKF is that it can handle only mild non-linearity (approximately linear within a time update) and hence not reliable. Also the algorithm requires computation of the Jacobian which is can be quite computationally intensive for large dimensional systems.
So far, the UKF has proven to be the best version of Kalman Filter for the non-linear process model. It is a deterministic sampling method whereby sampling points are selected strategically. The sampling points when they propagate they capture the first and second moment up to the 3rd order for a Taylor serial expansion of the non-linear function. It is based on the simple concept that it is better to approximate the posterior and use exact function rather than approximating the function and using the exact posterior. Hence incorporation of UKF in SIR algorithm improves the performance greatly.
One of the most important concerns and the main motivation behind our study is the computational effort required in a nonlinear state estimation. It can be seen that if we choose to draw samples by propagating through EKF or UKF based estimator, we have to run the estimation algorithm for each sample separately. Although the estimation accuracy will improve significantly but so the computation time. In this paper we propose to integrate the Ensemble Kalman Filter (EnKF) into the Particle filter due to their similarity and also due to numerous advantages of EnKF. The EnkF similar to the particle Filter algorithm works on a set of sample called as ensemble, and is based on the Monte-Carlo integration. The subsequent sections are dedicated to the design of an EnKF based Particle Filter algorithm that greatly reduces the computation without compromising the accuracy.

III. ENSEMBLE KALMAN FILTER
The EnKF is very similar to Particle Filter in a way that a continuous variable (states) defined by posterior density function p(x t |y t ) is represented using discrete random samples except the associated importance weights. The basic distinction between EnKF and Particle Filter is that with EnKF, the samples are updated by moving them toward the likelihood area once the current measurement is received. The Particle Filter, on the other hand, recalculates weights depending on their likelihood density making EnKF's tracking ability better than basic Particle Filtering algorithm since the samples always overlap with likelihood. There are some recent development in Ensemble Kalman Filtering such as Augmented EnKF [17] to overcome the limitations of linear updating. In EnKF, the error statistics are predicted by Monte Carlo integration to solve the Fokker-Plank equation. Consider a discrete-time non-linear state evolution model, and measurements as where x t , η t ∈ R n×1 are n-dimensional system state and process disturbance respectively at a discrete time instant t. y t , t ∈ R p×1 are p-dimensional measurement variable and noise signal respectively. It is assumed that η t and t are mutually and serially independent sequence of data. The goal is to obtain the posterior distribution p(x t |y t ) at time t given the measurement data up to point t. Once the posterior distribution is available, an updated state estimate x t−1 can be easily obtained.
Assume that the samples where µ t−1 and Σ t−1 are sample mean and sample variance at time t. The samples from the forecast distribution can be obtained by simply applying the state evolution as The sample mean of the forecast state is obtained as Similarly, the forecast measurement and the corresponding sample mean can be obtained as Let ε The Kalman gain is evaluated by using the following covariance matrices The updated samples can be obtained as The samples thus obtained represent the posterior distribution of the state at time t. Notice that the linear update equation (25) is the same as the one used for analytical Kalman Filtering, with the exception that the predicted observation is perturbed by the observation noise. The analysed samples can be said to be distributed as N x i t : µ t|t , Σ t|t , where Σ t|t is the empirical variance computed from analysed samples.

IV. ENSEMBLE KALMAN FILTER PARTICLE FILTER
Using the descriptions given in previous sections, we now present a hybrid SIR algorithm that uses EnKF as the proposal density. Because the basic EnKF works with samples with equal weights, it is relatively more straight forward to incorporated into the SIR algorithm if we choose to re-sample at each iteration. The re-sampling after the weight update stage will result in all weights being equal. However, we present here an alternative solution that employs a slightly modified version of EnKF that works with weighted samples and eliminates the need for weight equalisation at each iteration.
We will assume that the evolution of state vector is described as per the model given in Eq. 1 and 2. As described in Section-II, let the posterior density p(x t |y t ) at time instant t be again represented by samples x The predicted ensemble measurement y (i) t perturbed with Gaussian noise (i) t ∼ N (0, R) can be obtained by passing the forecasted ensemble members through the measurement equation The forecasted sample mean of the weighted states and measurement can be obtained as Let ε (i) t and e (i) t represent the state and measurement estimation errors at time t Following covariance matrices for weighted samples can be evaluated (with Bessel's correction) If the weights ω (i) t−1 are chosen to be normalised such that The covariance matrices can be redefined as The Kalman gain can be evaluated using the above matrices The updated state (analysed ensemble members) are obtained by the following linear combination equation The resultant samples can be said to be distributed as per the proposal distribution where The weight update equation in Particle Filter can be now written as Thus, the proposal density of EnKFPF is obtained using a common empirical mean and variance, as opposed to EPF or UPF, where the mean µ  t|t . The proposal density for EPF and UPF is given by Another crucial component of particle filtering is the convergence property, which assures that the estimate error approaches the precise analytical answer as the number of samples approaches infinity. The convergence in Particle Filtering has been discussed in general by Crisan and Doucet (2002) [18], [19], [20], and [21]. The sufficient condition for empirical convergence for a Particle Filter in L p -sense requires to assume that the process and observation models are both bounded, as are the importance weights (unnormalized). The following conditions must be met for any Borel test function φ(x).
for some constant c k . The suggested EnKFPF method, like the EPF and UPF algorithms, pulls samples that adequately overlap the likelihood p(y t |x i t ) and hence weights are not expected to grow unbound. The EnKFPF algorithm thus fulfils the sufficient condition defined above.
The pseudo-code for the basic Sequential Importance Resampling (SIR) algorithm is described in Algorithm 2.

V. RESULTS AND DISCUSSIONS
In this section, we present numerical results illustrating performance of the proposed EnKFPF algorithm for various scenarios, and compare with EKF and UKF proposals density based Particle Filtering algorithms. We consider the following non-linear time series model that is often used in the literature to compare numerical filter techniques [22], [1].
where η t ∼ N (0, Σ η ) and t ∼ N (0, Σ ) represent zero mean Gaussian process and measurement noise. The above state evolution model quickly bifurcates the prior p(x t |x t−1 ) transforming it into a bi-model density function if the posterior density p(x t−1 |y t−1 ) is assumed to be zero mean Gaussian. It's also interesting to note that when the measurement yields negative values, the likelihood density's shape shifts from bimodel to zero mean uni-model. For such a bi-modal, nonlinear scenario, any algorithms based on the assumption of Gaussianity involving linear update may generally result in a less accurate estimate and will quickly diverge. The first set of simulations were run with Gaussian assumptions for both the process and the observation model for various sample sizes and uncertainties. For simplicity, the simulations were run without the Metro-police Hastings step. The process and measurement error are assumed to be zero mean white Gaussian with σ η = 1 and σ = 1. The results for all three algorithms are demonstrated in Fig. 1, for sample sizes = 1000. For this particular case, UPF appears to be providing better estimate than EnKFPF. For further quantification of the convergence in MSE error, simulations with various sample sizes were carried out. The results are shown in Fig. 2.  Fig. 2 illustrates the efficacy of proposed EnKFPF algorithm along with EPF and UPF in terms of error convergence as a function of sample size for process and observation noise standard deviations σ η and σ equal to 1. The larger sample sizes is showing to improve the results as expected and asymptomatically converges to certain value as the sample size is infinite (N → ∞). In case of EnKFPF, all the ensemble member are shifted towards the observation using the empirical covariances. As a results the EnKFPF is expected to generate proposal densities q(x t |x t , y t ) that sufficiently overlaps with the p(y t |x t )p(x t |x t−1 ) resulting into faster convergence.
It is reasonable to assume that the MSEs of all the algorithms converge to the analytical method as N approaches to infinity (N → ∞). It can be observed that the UPF provide better results for almost all case as compare with EnKFPF. This is due to the fact that UPF propagates each sample through UKF estimator and hence samples are drawn from individually evaluated density parameters. When large number of samples are used, however, it is clear that EnKF is more efficient in achieving practically identical accuracy. When compared to EPF, both EnKFPF and UPF perform much better. The inaccurate approximation of the extremely non-linear model is to account for the EPF's poor estimation accuracy. The computation time for all the three algorithm and MSE values for 1000 samples with σ η = 1 and σ = 1 is shown in terms of a bar graph in Fig. 3. The crucial element to note here is the EnKFPF's substantially shorter computation time, which corresponds to the assertions stated in the preceding sections. As state earlier, the fact that EnKFPF does not run the computationally intensive estimation algorithms for each sample and instead uses empirical variances results in a significant reduction in computation time. It should be noticed that the EPF has a computation time that is nearly identical to the EnKFPF. However, the computation time in EPF is significantly influenced by the model's structure and dimensions, and hence can vary greatly depending on the case.

EPF
UPF EnKFPF  Another set of trials was conducted by changing the process uncertainty from Gaussian to non-Gaussian density while maintaining Gaussian measurement noise. We choose the process noise uncertainty to be Gamma distributed with parameters α = 3 and θ = 1/2. The mean and variance are αθ = 3/2 and αθ 2 = 3/4, respectively, which are used to generate estimated states in the EPF and UPF methods. The EnKFPF, on the other hand, makes use of the empirical variances calculated from the samples. With all three algorithms, the estimated states are presented in Figure-4 for a sample size of 1000. The analysis reveal that UPF gives superior accuracy for small measurement noise in case where process uncertainty is described by non-Gaussian non-zero mean density function case. Due to broader overlap with likelihood, the EnKFPF outperforms UPF in terms of accuracy and computation time for large measurement noise values (e.g., σ = 10). For a clear understanding, we provided the estimation results for several Gaussian and non-Gaussian scenrios Table-I, which is in terms of MSE error and computation time averaged over 100 trials. Based on the inaccuracy and computation time in the table, EnKFPF looks to be an efficient alternative that provides both accuracy and a significant reduction in computing cost. EnKFPF's accuracy improves as the level of uncertainty increases, which is another desirable aspect.   In this article, a version of Ensemble Kalman Filter-based proposal density has been suggested in order to realise an efficient Particle Filter that does not require large computational effort. The proposed algorithm draws the samples straight away from the ensemble mean and variance evaluated in EnKF's analysis step and uses them to update the weights. The proposed approach thus effectively combines the benefits of both analytical and Monte Carlo approaches to improve accuracy. We ran several scenarios on the benchmark nonlinear problem to demonstrate that the proposed methodology is on par with the UPF in terms of accuracy and convergence while requiring far less computational effort. The deterministic shift in the samples towards measurement makes it suitable for many real world on-line application where tracking ability is required along-with lower computation effort.
Finally, we note that the efficacy of our proposed algorithm is shown for specific problem using scaler-non-linear case only with Gaussian and bi-modal uncertainties and can be extended to multi-modal cases. We believe that the proposed method may be further studies in details with respect to the large dimensional problems in general.