Nonlinear Kalman Filtering with Divergence Minimization

I Introduction

Modeling and analysis of time-varying signals is one of the most important subfields of signal processing. The problem arises in many different forms, such as communications data sent over a channel, video and audio data, and real-time tracking. A wide variety of algorithms have been developed in the statistics and engineering communities to deal with such dynamic systems. One classic algorithm is the Kalman filter , which performs minimum mean square error estimation of the hidden state of a time-varying linear system. Kalman filter is recursive and online, making it suitable for real-time signal processing applications. Another advantage is its optimality for a large class of state-space models.

Kalman filtering has been applied extensively in control, communication, and signal processing settings, such as robot motion control and radar target tracking. With the recent explosions in sequential and streaming data, Kalman filters have also become a promising means for approaching machine learning problems, such as natural language processing , collaborative filtering and topic modeling .

An important issue that often arises, requiring modification to basic Kalman filter framework, is nonlinearity. For example, in radar tracking, distance and bearing measurements require a Cartesian-to-polar transformation , whereas dynamic collaborative filtering model contains a bilinear form in two unknown vectors . The nonlinear problem has been studied extensively in the literature, resulting in well-known filtering algorithms such as the extended Kalman filter (EKF) and unscented Kalman filter (UKF) . On the other hand, Monte Carlo methods have been developed , which are non-parametric and can represent any probability distribution using a discrete set of points, also referred to as particles.

While particle filters can approximate arbitrary densities, it may still be important to find the best parametric distribution according to a particular objective function. This has been a major goal in Bayesian learning, where the exact posterior distribution is usually intractable and approximated by a known, “simpler” distribution. Two established ways to handle this problem are variational inference and expectation-propagation , in which the Kullback-Leibler (KL) divergence between the true posterior and the approximating distribution are minimized. Ideas from approximate inference have also been used in the Kalman filtering framework . However, a thorough analysis of posterior optimization for nonlinear Kalman filters have not yet been made.

In this paper we fill this gap by presenting three algorithms for nonlinear Kalman filtering based on three respective divergence measures for posterior approximation, each based on a parametric form (in our case, a multivariate Gaussian). These approximations are obtained by algorithms for approximation-free divergence minimization. The divergence measures we consider are: 1) the forward KL divergence as used in variational inference; 2) the reverse KL divergence as used in expectation-propagation; and 3) the $\alpha$-divergence, which is a generalized family that contains the former two as special cases. We also show that well-known algorithms such as the EKF and UKF are actually solving approximations to KL divergence minimization problems. This further motivates our study to address these shortcomings.We emphasize that our methods are still approximate in that the true non-Gaussian posterior will be approximated by a Gaussian. It is approximation-free in that the three algorithms directly optimize the three divergences.

The main machinery we use for obtaining these unbiased minimization algorithms is importance sampling. However, the resulting algorithms are all computationally cheaper than particle filtering since 1) no resampling is necessary, and 2) the number of unnecessary samples can be reduced by our proposed adaptive sampling procedure. We show advantages of our algorithms for target tracking and options pricing problems compared with the EKF, UKF and particle filter.

We organize this paper as follows: In Section II we define our filtering framework by reviewing the Kalman filter and discussing its non-linear variants. In particular, we discuss parametric approaches, also called assumed density filters, and nonparametric approaches, also called particle filters. In Section II-B we present three divergence minimization problems based on the forward and reverse KL divergence, and $\alpha$-divergence. For each case we propose an algorithm which minimizes the corresponding objective function. Our algorithms are based on Monte Carlo integration techniques. Section IV contains a number of experiments to show how these divergence measures compare with each other and with standard approaches. Finally we conclude in Section 5.

II Kalman Filtering

The Kalman filter has been developed and motivated as an optimal filter for linear systems. A key property is that this optimality is assured for general state-space models. This has made Kalman filtering widely applicable to a wide range of applications that make linearity assumptions. The Kalman filter can be written compactly at time step $t$ as

The two main tasks of Kalman filtering are prediction and posterior calculation ,

When the initial distribution on $p(x)$ is Gaussian all these calculations are in closed form and are Gaussians, which is an attractive feature of the linear Kalman filter.

II-B Nonlinear framework

In many problems the measurements $y_{t}$ involve nonlinear functions of $x_{t}$. In this case the Kalman filter becomes nonlinear and the closed-form posterior calculation discussed above no longer applies. The nonlinear process is

where the noise process is the same as in Eq. (1), but $h(\cdot)$ is a nonlinear function of $x_{t}$.We focus on measurement nonlinearity in this paper, assuming the same state space model. The techniques described in this paper can be extended to nonlinearity in the state space as well. While formally Bayes’ rule lets us write

the normalizing constant is no longer tractable and the distribution $p(x_{t}|y_{1:t})$ is not known. Although the nonlinearity in $h$ may be required by the problem, a drawback is the loss of fast and exact analytical calculations. In this paper we discuss three related techniques to approximating $p(x_{t}|y_{1:t})$, but first we review two standard approaches to the problem.

II-C Parametric approach: Assumed density filtering

To address the computational problem posed by Eq. (4), assumed density filters (ADF) project the nonlinear update equation to a tractable distribution. Building on the linear Gaussian state-space model, Gaussian assumed density filtering has found wide applicability . The main ingredient here is an assumption of joint Gaussianity of the latent and observed variables. This takes the form,

(We’ve suppressed some time indexes and conditioning terms.) Under this joint Gaussian assumption, by standard computations the conditional distribution $p(x_{t}|y_{t})$ is

In this case, the conditional distribution is also the posterior distribution of interest. Using this approximation, Kalman filtering can be carried out. For reference we provide predictive update equations in Appendix A.

There are several methods for making this approximation. We briefly review the two most common here: the extended Kalman filter (EKF) and the unscented Kalman filter (UKF). The EKF approximates $h$ using the linearization

where $H(x_{0})$ is the Jacobian matrix evaluated at the point $x_{0}$. For example, $x_{0}$ could be the mean of the prior $p(x_{t}|y_{1:t-1})$. By plugging this approximation directly into the likelihood of $y_{t}$, the form of a linear Kalman filter is recovered and a closed form Gaussian posterior can be calculated.

There are many extensions to the UKF framework such as cubature Kalman filtering (CKF) and QMC Kalman filtering , which use different numerical quadratures to carry out the approximation, but still correspond to the joint Gaussian assumption of Eq. (9). With that said, however, not all Gaussian ADFs make a joint Gaussianity assumption. For example, methods based on expectation-propagation use moment matching (e.g., ) to obtain a Gaussian posterior approximation without modifying the joint likelihood distribution. We focus on an EP-like method in Section III-B.

II-D Nonparametric approach: Particle filtering

The positive weights $w_{t}^{s}$ sum to one, and $\delta_{x_{t}^{s}}$ is a point mass at the location $x_{t}^{s}$. The main approach is to use particle filters, a method base on importance sampling. In case of particle filtering using sequential importance resampling (SIR) , updating an empirical approximation of $p(x_{t}|y_{1:t})$ uses a uniform-weighted prior approximation, $p(x_{t})\approx\sum_{s=1}^{N_{s}}\frac{1}{N_{s}}\delta_{x_{t}^{s}}$ to calculate the posterior importance weights

It then constructs the uniform-weighted prior approximation by sampling $N_{s}$ times

While SIR particle filters can adaptively approximate any posterior density, the double sampling has computational cost, making these filters considerably slower compared to the above ADF approaches. Another potential issue is the need to propagate particles between time frames, which can be prohibitively expensive in communication-sensitive distributed applications, such as sensor networks .

III Three divergence minimization approaches

Given two distributions $p(x|y)$ and $q(x)$, the forward KL divergence is defined as

The KL divergence is always nonnegative, becomes smaller the more $q$ and $p$ overlap, and equals zero if and only if $q=p$. These properties of the KL divergence make it a useful tool for measuring how “close” two distributions are. It is not a distance metric however, as $\text{KL}[q\|p]\neq\text{KL}[p\|q]$; we discuss the latter in detail in Section III-B. In Bayesian machine learning, minimizing an objective of this form over $q$ is know as variational inference (VI) . In this case, $p(x|y)$ corresponds to an unknown posterior distribution of the model parameters, and $q$ is its simpler approximation.

For the nonlinear Kalman filtering problem, the posterior is on the latent state vector $x_{t}$ and so is intractable. Therefore, as is often the case, $\text{KL}[q\|p]$ is not calculable. Variational inference instead uses the identity

This often is tractable since the joint distribution $p(y,x)$ is defined by the model. Since the marginal $\ln p(y)$ is constant and $\text{KL}\geq 0$, variational inference instead maximizes $\mathcal{L}$ with respect to parameters of $q(x)$ to equivalently minimize KL.

The terms in the second line are tractable, but in the first line the nonlinearity of $h(x_{t})$ will often result in an integral not having a closed form solution.

In the variational inference literature, common approaches to fixing this issue typically involve making tractable approximations to $h(x_{t})$. For example, one such approximation would be to pick a point $x_{0}$ and make the first-order Taylor approximation $h(x_{t})\approx h(x_{0})+H(x_{0})(x_{t}-x_{0})$. One then replaces $h(x_{t})$ in (III-A) with this approximation and optimizes $q(x_{t})$. In fact, in this case the resulting update of $q(x_{t})$ is identical to the EKF. This observation implies a correspondence between variational inference and commonly used approximations to the non-linear Kalman filters such as the EKF. We make this formal in the following theorem.

Theorem 1: Let the joint Gaussian ADF correspond to the class of filters which make the joint distribution assumption in (9). Then, all filters in this class optimize an approximate form of the variational lower bound in (III-A).

We present a complete proof in Appendix B. Theorem 1 is general in that it contains the most successfully-applied ADFs such as the EKF and UKF, among others. For the special case of EKF, the nature of this approximation is more specific.

Corollary 2: The EKF corresponds to optimizing the objective (III-A) using a first order Taylor approximation of $h$.

Please see Appendix C for a proof. Consequently, the existing algorithms modify $\mathcal{L}$ and the optimization of this approximation to $\mathcal{L}$ over the parameters of $q(x_{t})$ is no longer guaranteed to minimize $\text{KL}[q\|p]$. Instead, in this paper we are motivated to fill in this gap and find ways to directly optimize objectives such as (III-A), and thus minimize divergence measures between $q$ and the intractable posterior $p(x_{t}|y_{t})$. We next devise a method for $\text{KL}[q\|p]$.

Recently Paisley, et al. proposed a stochastic method for sampling unbiased gradients of $\mathcal{L}$, allowing for approximation-free minimization of the forward KL divergence using stochastic gradient descent. We derive this technique for the nonlinear Kalman filter, which will allow for approximate posterior inference having smaller KL divergence than the EKF and UKF. Using simpler notation, we seek to maximize an objective of the form,

We let $\mu_{t}$ be the current value of the mean of $q(x_{t})$ at a given iteration of time $t$. If we define $\widetilde{y}_{t}=y_{t}-h(\mu_{t})+H(\mu_{t})\mu_{t}$, then equivalently we can write

The expectation is now in closed form. While a better approximation may have greater variance reduction for a fixed number of MC-samples, we emphasize this would not make the algorithm more “correct.” Where the EKF simply replaces $f(x)$ with $g(x)$, our stochastic gradient approach then corrects the error of this approximation.

III-B Approach 2: Reverse KL divergence minimization

As mentioned in Section III-A, KL divergence is not a distance measure since it is not symmetric. The complement of the forward KL divergence defined in (14) is the reverse KL divergence:

We can see that (29) offers an alternative measure of how similar two probability distributions are; therefore we can use it to approximate an intractable posterior distribution.

Note that for either objective function, (14) or (29), the optimal solution will be $q(x)=p(x|y)$. However, since the approximating distribution is typically different from the exact posterior distribution, the two optimization problems will give different solutions in practice. In particular, reverse KL divergence has shown to be a better fit for unimodal approximations, while forward KL works better in multimodal case . Consequently, we can expect that optimizing the reverse KL will be a better choice for the nonlinear Kalman filtering problem (this is supported by our experiments). In Section III-A, finding a fixed point of the forward KL problem required an iterative scheme for maximizing the variational objective function. The fixed point of the reverse KL has a more interpretable form, as we will show.

To this end, we first note that an exponential family distribution has the form

where $\eta$ is the natural parameter and $s(x)$ is the sufficient statistic. Therefore inference in exponential families correspond to determining $\eta$. Substituting this parametrized form in (29) and setting the derivative with respect to the natural parameter equal to zero, one can show that

This moment matching is well-known in statistics, machine learning, and elsewhere . In machine learning it appears prominently in expectation propagation .

A common choice for the approximating exponential family distribution is again Gaussian because it is the maximum entropy distribution for the given first and second order moments . Since a Gaussian is completely specified by its mean and covariance, when the approximating distribution $q(x)$ is selected to be Gaussian, the optimal solution is simply found by matching its mean and covariance to that of $p(x|y)$.

III-C Approach 3: α𝛼\alpha-divergence minimization

In Sections III-A and III-B we showed how nonlinear Kalman filtering can be performed by minimizing the forward and reverse KL divergence. A further generalization is possible by considering the $\alpha$-divergence, which contains both KL divergences as a special case. Following we define the $\alpha$-divergence to be

where the parameter $\alpha$ can take any value in $(-\infty,\infty)$. Some special cases are

where $\text{Hel}[p\|q]$ is the Hellinger distance. We see that when $\alpha=1/2$ we get a valid distance metric. Similar as before, we now seek a $q$-distribution which approximates $p(x|y)$, where approximation quality is now measured by the $\alpha$-divergence.

Again assuming that the approximating distribution is in the exponential family, $q(x)=h(x)\exp\{\eta^{\top}s(x)-\log A(\eta)\}$. The gradient of the $\alpha$-divergence shows that

Note that we defined a new probability distribution $\widetilde{p}(x)=p(x)^{\alpha}q(x)^{1-\alpha}/Z_{\widetilde{p}}$ where the denominator term is the cumulant function. This leads to a generalized moment matching condition,

This problem is more complicated than the reverse KL because the left hand side also depends on the $q$-distribution. The $\alpha$-divergence generalizes a number of known divergence metrics. In context of EP, it is possible to obtain a generalization which is called Power EP . More recently, used a similar black-box optimization, where they showed that by varying the value of $\alpha$ the algorithm varies between variational inference and expectation propagation. It turns out that, for many practical problems, using a fractional value of $\alpha$ can give better performance than the limiting cases $\alpha=0$ or $\alpha=1$. This motivates our following $\alpha$-divergence minimization scheme.

A similar importance sampling methodology can be used for this optimization as for the reverse KL divergence. Using similar notation, we can write

where $x^{s}\sim_{iid}q(x)$. Again we define

We see that the procedure in (31) is a special case of this when we set $\alpha=1$. However, there is a significant difference in that the moment matching of (30) can be done in one iteration since it only depends on $p$. In (36) the $q$ distribution appears on both sides of the equality. This is similar to of EP and Power-EP algorithms, where multiple iterations can be run to update $q$. Upon convergence we know that the solution is a fixed point of (32), but convergence of the procedure is not guaranteed and multiple iterations might degrade the performance. In our experiments we will only iterate once to avoid possible diverging and also to keep the cost of the algorithm the same as that of MKF in the previous section. We call this algorithm $\alpha$-divergence Kalman filter ($\alpha$KF ) and summarize it in Algorithm 3. We note that the only difference between $\alpha$KF and MKF is in step 4.

We can get a better understanding of employing $\alpha$-divergence by analyzing the weight coefficients. In particular, lets assume that we choose our proposal distribution as the prior, i.e. $\pi(x)=p(x)$. Then, the MKF weights become $w^{s}\propto p(y|x^{s})$ in the Kalman filter. The $\alpha$KF weights, on the other hand are $w^{s}\propto p(y|x^{s})^{\alpha}$; therefore, the likelihood term is scaled by alpha and as $\alpha\rightarrow 0$ all the particles generated will have equal contribution. For very low values of $\alpha$ this will discard all the information, which is clearly unwanted, but for intermediate values this can alleviate the effects of sharply fluctuating likelihood factors. As we will show in our experiments, when the measurement noise is strong, choosing an intermediary $\alpha$ value provides robustness.

III-D Adaptive Sampling

The main parameter in the implementation of sampled filters such as particle filters and the three filters proposed here is the number of particles $S$ that will be used. Hence it is desirable to have a method of estimating the minimum number of samples necessary for a given degree of accuracy. Then, for each round of filtering we can use this computed sample size to reduce the computation as much as possible, but still be able to increase the sample size when necessary. In Figure 1 we illustrate the problem of tracking a moving target. At time $t$ this target makes an abrupt maneuver where we need more particles for accurate tracking, but we can reduce the size afterwards.

where $\chi_{d}^{2}(p)$ is the quantile function of chi-squared distribution with d degrees of freedom (which equals the state-space dimension here), and $p$ is the probability value (for 95% confidence intervals this is set to $p=0.95$). $\chi_{d}^{2}$ indicates the chi-squared distribution. The region described by (38) is a hyper-ellipsoid, so the maximum possible radius will correspond to the major axis, which is given by

Note that this is a conservative estimate, as the hypersphere with radius $r_{\text{max}}$ will typically be much larger than the hyper-ellipsoid. An illustration of the bounding circle for 2D multivariate normal distribution is given in Figure 1.

Now assume that using a small sample set $S_{\text{base}}$ we wish to estimate the minimum number of samples $S_{\text{min}}$ required to achieve a certain $r_{\text{max}}$. We have the relation and result that

IV Numerical Results

We experiment with all three proposed nonlinear Kalman filter in algorithms, as well as the EKF, UKF and particle filter, on radar and sensor tracking problems, as well as an options pricing problem.

The first problem we consider is target tracking. This problem arises in various settings, but here we consider two established cases: radar and sensor networks. The radar tracking problem has been a primary application area for nonlinear Kalman filtering. The target is typically far away from the radar, for example an airplane. Wireless sensor networks are another emerging area where nonlinear filtering is useful. Driven by the advances in wireless networking, computation and micro-electro-mechanical systems (MEMS), small inexpensive sensors can be deployed in a variety of environments for many applications .

For both problems the state-space will have the form

Here, $F_{t}$ and $Q_{t}$ model the dynamics of target motion and are usually time-varying. On the other hand, $h(\cdot)$ specifies the equipment that performs the measurements, and the environment and equipment based inaccuracies are represented by $R_{t}$. In the radar setting, when the target is far away and the angle measurement noise is strong enough, the problem can become highly nonlinear. For sensor networks, the nonlinearity is caused by the small number of active sensors (due to energy constraints) with large measurement noise (due to the attenuation in received signal) . While the value of $R_{t}$ can be determined to some extent through device calibration, it is more challenging to do this for $Q_{t}$ .

The radar measures the distance and bearing of the target via the nonlinear function $h(\cdot)$ of the target location,

i.e. the Cartesian-to-polar transformation . For the sensor networks, we will consider a scenario which uses range-only measurements from multiple sensors. This yields the model in (41) where $h(\cdot)$ is the measurement function such that the $i$-th dimension (i.e. measurement of sensor $s_{i}$) is given by

and the length of $h(x_{t})$ will be the number of activated sensors at time $t$.

We consider two types of problems: tracking with uncertain parameters and tracking with known parameters. For the case of uncertain parameters, we set the radar and sensor simulation settings as follows. First, we note that for both simulations we assume a constant measurement rate, and so set $\Delta t=1$. For radar we sweep the process noise values in (42) as $\sigma_{CV}\in\{10^{-3},~{}2\times 10^{-3},~{}\ldots~{},~{}10^{-2}\}$. We generate 20 data sets for each value of $\sigma_{CV}$, yielding a total of 200 experiments. For the measurement noise we use a diagonal $R$ with entries $\sigma_{r}^{2}=10^{-1}$ and $\sigma_{\theta}^{2}=10^{-2}$ which dictates the noise of distance and bearing measurements respectively. The initial state is selected as $x_{0}=[1000,~{}10,~{}1000,~{}10]^{\top}$; this distance from origin and angle noise variance results in a severely nonlinear model, making filtering quite challenging. For sensor network simulations, we use the same constant-velocity model of (42) with $\sigma_{CV}=10^{-2}$. We deploy 200 sensors and at each time there are exactly 3 distinct ones responsible for range measurements. The measurement covariance matrix is $R=\sigma_{R}^{2}I$ where we set $\sigma_{R}=20$. We select the initial state as $x_{0}=[1000,~{}1,~{}1000,~{}1]^{\top}$. With this, once again, we obtain a highly nonlinear system, albeit less severe than the radar case. We also consider the case where the generating parameters are known to the filter. In this case, we assign the performance of the filter as a function of process and measurement noise covariances. For this one, we sweep $\sigma_{\text{CV}}\in\{0.001,0.005,0.01,0.05,0.1\}$ and $\sigma_{r}\in\{10,15,20,25,30\}$. We report the results for the sensor network case.

We implemented EKF, UKF, sampling-importance-resampling particle filter (PF), and our proposed SKF, MKF, and $\alpha$KF for $\alpha=0.5$. For SKF we use $500$ particles/iteration, whereas we consider $10^{4}$ particles for PF, MKF and $\alpha$KF . When there is parameter uncertainty, the exact value of $Q$ is not known to the filter, therefore we consider a scaled isotropic covariance of form $\sigma_{Q}^{2}I$.

In Table I we show mean square error (MSE) for radar tracking as a function of the selected scale value ($\sigma_{Q}$). Here, the base error corresponds to the estimations based on measurements only, and its order-of-magnitude difference from filter MSE values show the severity of nonlinearity. Now, comparing MSE values, first we see that MKF and $\alpha$KF overperforms EKF and UKF for all settings of $\sigma_{Q}$ which shows that the Gaussian density obtained from these filters is indeed more accurate. SKF also gets better results, particularly for $\sigma_{Q}=10^{-1}$ but it is less robust to the changes in scale value. This is due to the iterative gradient scheme employed by SKF, which could give worse results depending on parameter changes or covariance initializations. Since MKF/$\alpha$KF are based on importance sampling, they do not exhibit the same sensitivity. As for PF, this algorithm also produces competitive results when $\sigma_{Q}=10^{-2}$; however its performance significantly deteriorates (even more than that of SKF) as $\sigma_{Q}$ increases, which shows the nonparametric inference of particle filtering is more sensitive to parameter uncertainty. We also mark the best overall MSE with boldfaces, which is given by $\alpha$KF for $\sigma_{Q}=10^{-1}$. Furthermore, $\alpha$KF has the highest robustness to parameter changes, therefore it is a better candidate to choose when parameters are not known and measurements are very noisy, since the $\alpha$ coefficient has the capability of mitigating excess measurement noise, as discussed in Section III-C.

Table II presents MSE results for sensor networks. Unlike the radar problem, all particle-based filters are better than EKF/UKF for all values of $\sigma_{Q}$. This reduced sensitivity is due to the reduced nonlinearity in the problem. The performance of SKF, MKF, and PF are similar to each other, MKF being the favorable choice for most of the cases. On the other hand $\alpha$KF is the best performer in all cases, and as $\sigma_{Q}$ increases, the margin increases. The best overall MSE is again achieved by this filter for $\sigma_{Q}=10^{-1}$, where using $\alpha$KF provides a clear benefit.

In Figure 2 we show qualitative tracking results from sensor networks. The top, middle, and bottom rows correspond to SKF, MKF, and $\alpha$KF respectively. For each two we pick four different paths (shared across different rows) and for each plot we plot the true trajectory along with EKF, UKF, and one of our filters, depending on the row. First we see that our simulation settings encompass a wide variety of paths which exhibit multimodality such as, for example, a combination of constant velocity and constant turn models . By visual inspection we can see that our algorithms provide more accurate tracking compared to EKF/UKF in all cases. Furthermore, moving down the rows we can see that the accuracy of our filtering algorithms also increase and the $\alpha$KF estimated paths are more robust to measurement errors, as clearly demonstrated in second column.

So far, for $\alpha$KF we only considered the case when $\alpha=0.5$, which used the symmetric Hellinger distance metric, as given in (III-C). Now we focus on varying the value of $\alpha$ and analyzing its effects. For this we use the sensor experiments with $\sigma_{Q}=10^{-1}$, which corresponds to the mid column of Table II. The mean squared error as a function of $\alpha$ is plotted in Figure 3. We see that low-mid ranges of $\alpha$ (i.e $0.3-0.5$) give the best MSE results. This improvement is obtained since lower values of $\alpha$ help mitigate the effects of strong measurement noise. There is, however, a tradeoff here since choosing a too small value for this parameter will discard all the measurement information and give poor results. This is seen for lower values of $\alpha$, where decreasing the parameter degrades performance.

As discussed in Section III-D we can use adaptive sampling to choose the minimum possible sample size to achieve a certain confidence region radius, $r_{\text{max}}$. We implemented adaptive sampling for $\alpha$KF using an initial batch size of $S_{\text{base}}=500$. We picked four different values of $r_{\text{max}}$ from $\{0.5,1,1.5,2\}$. Figure 4 displays the results for this experiment. In the left panel we compare the MSE results as a function of $r_{\text{max}}$ for $\alpha$KF and PF for the sensor tracking problem with $\sigma_{\text{CV}}=0.1$. Note that, for PF, adaptive sampling is not a choice as all particles should be propagated, resampled, and updated at every time step. So for PF we simply set the sample size as the average $S_{\text{min}}$ for the $\alpha$KF for each case. We can see that, the MSE performances differ very little across different cases, showing even for larger target values of $r_{\text{max}}$ both methods can still produce accurate estimates of the true state. We also see that $\alpha$KF overperforms PF in all cases. On the other hand, the right panel shows the number of samples required to achieve a certain confidence radius. From this figure we can see the $\mathcal{O}(1/r^{2})$ decaying rate of $S_{\text{min}}$ as implied by (40). Given the high accuracies in the left panel, we see that several hundred samples can be sufficient to obtain high-quality estimates, which makes $\alpha$KF competitive for real time applications. Another point is, as $r_{max}$ increases, the variance of the sample size also decreases, which means the runtime per round will have small discrepancy, as opposed to using a smaller $r_{\text{max}}$.

We now turn to the case where the process noise parameter is known. In Figure 4 we show the filter MSEs as a function of $\sigma_{CV}$ and $\sigma_{R}$. For the measurement noise, as $\sigma_{R}$ increases the overall MSEs also increase, while for process noise this trend is not present. For both cases we see that the particle filter gives the best result overall. This is expected, since when the parameters are known perfectly, particle filter can approximate the posterior with more accuracy, as it is nonparametric. With that said, $\alpha$KF is also competitive in this setting. In fact, for several cases such as $\sigma_{\text{CV}}=0.001$ and $\sigma_{r}=25$ performance of $\alpha$KF and PF are equal, and for the remaining cases the particle filter does not improve much compared to $\alpha$KF , while both filters can perform much better than SKF and MKF. This means $\alpha$KF can be preferred over PF, since it does not require resampling. As a second observation, note that SKF/MKF perform much better than EKF/UKF, and $\alpha$KF perform even better compared to the rest. This means, by minimizing different forms of divergence one can indeed get significantly better Gaussian approximations of the posterior, which supports our theoretical analysis in Section III.

IV-B Options Pricing

We also consider a problem in options pricing. In finance, an option is a derivative security which gives the holder a right to buy/sell (call/put option) the underlying asset at a certain price on or before a specific date. The underlying asset can be, for example, a stock. The price and date are called the strike price and expiry date respectively. The value of the option, called premium, depends on a number of factors. Let $C$ and $P$ denote the call and put prices. We use $\sigma$ and $r$ to denote volatility and risk-free interest rate respectively; the values of these variables are not directly observed, hence they need to be estimated. Let $S$ denote the price of underlying asset and $X$ denote the strike price. Finally, let $t_{m}$ denote the time to maturity; this is the time difference between the purchase and expiry dates which is written as a fraction of a year. For example, an option which expires in two months will have $t_{m}=1/6$.

Accurate pricing of options is an important problem in mathematical finance. For a European style option, the price as a function of all these parameters can be modeled using the well-known Black-Scholes equation

Following the approach of , let $x={[\sigma~{}r]}^{\top}$ be the state and $y={[C~{}P]}^{\top}$ be the measurement. We get the following state space representation

where the nonlinear mapping $h(\cdot)$ is given by (44). In this case we model the process and measurement noises with time-invariant covariance matrices $Q$ and $R$. We consider two tasks: 1) predicting the one-step ahead prices, and 2) estimating the values of hidden state variables. This problem is also considered in to assess the performance of particle filtering algorithms.

Here we use the Black-Scholes model as the ground truth. In order to synthesize the data, we use historical values of VIX (CBOEINDEX:VIX), which measures the volatility of S&P 500 companies. From this list we pick Microsoft (NASDAQ:MSFT), Apple (NASDAQ:AAPL), and IBM (NYSE:IBM) as underlying assets and use their historical prices. The interest rate comes from a state-space model with a process noise of zero mean and variance $10^{-4}$. We set $\sigma_{Q}=\sigma_{R}=10^{-2}$. In Table IV-B we show the next-day prediction performance of all algorithms. We can see that the prediction performance imporves as we move towards MKF. This, again shows the difference between Gaussian approximations of the methods we employ. For MKF and PF we used $10^{3}$ particles, and their results were similar so we only report MKF here; however we also note that MKF can achieve this performance without using resampling, and it can leverage adaptive sampling to reduce sample size, which makes it preferable over PF. On the other hand, for SKF we need to use a large number of particles per iterations (around $1,000$). Even though this gives better results then EKF and UKF it is much slower than MKF/PF, and its performance can vary significantly between iterations, which makes it less competitive in this case. On the other hand, since the measurement noise is small in this case, choosing $\alpha<1$ for $\alpha$KF does not provide improvement over MKF in this case, which is consistent with our previous intuition. Therefore $\alpha=1$ is the best choice in this case.

Figure 5 shows the volatility estimation for three filters: Usually EKF tends to over/under-shoot a lot and UKF is significantly better in that respect; however MKF improves even further as it gives the most robust estimates. The plot of SKF output is similar to MKF. Also, similar to the target tracking experiments, we see that MKF has better performance than SKF, which once again agrees with the observation that expectation-propagation typically outperforms variational inference for unimodal posterior.

For the proofs in the following appendices we need the predict-update equations of the joint Gaussian ADFs. Note that this corresponds to the model in (9). The equations are summarized as

We emphasize that these hold for any joint Gaussian ADF. When EKF is employed, $H_{t}$ is the Jacobian at prior mean, and $S_{t}$ and $K_{t}$ are calculated accordingly.

The matrix inversion lemma asserts $(A+UCV)^{-1}=A^{-1}-A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}$; applying this to (52) we obtain

Substituting (54) into (53) and expanding we get

Once again, using the matrix inversion lemma we get