Nonlinear Kalman Filtering with Divergence Minimization

TL;DR

This paper introduces novel Monte Carlo-based algorithms for nonlinear Kalman filtering that directly optimize divergence measures like KL and alpha-divergence, improving accuracy over traditional filters.

Contribution

It proposes unbiased algorithms for divergence optimization in nonlinear Kalman filtering without approximations, enhancing filtering performance in complex applications.

Findings

Improved accuracy over UKF and EKF in tracking and pricing tasks.

Competitive performance with particle filtering.

Demonstrated effectiveness on radar, sensor tracking, and options pricing.

Abstract

We consider the nonlinear Kalman filtering problem using Kullback-Leibler (KL) and $\alpha$-divergence measures as optimization criteria. Unlike linear Kalman filters, nonlinear Kalman filters do not have closed form Gaussian posteriors because of a lack of conjugacy due to the nonlinearity in the likelihood. In this paper we propose novel algorithms to optimize the forward and reverse forms of the KL divergence, as well as the alpha-divergence which contains these two as limiting cases. Unlike previous approaches, our algorithms do not make approximations to the divergences being optimized, but use Monte Carlo integration techniques to derive unbiased algorithms for direct optimization. We assess performance on radar and sensor tracking, and options pricing problems, showing general improvement over the UKF and EKF, as well as competitive performance with particle filtering.

Figures

Tables

Table 1. TABLE I: Radar tracking problem: Mean Square Error (MSE) of various filtering schemes as a function of process noise parameter σ Q subscript 𝜎 𝑄 \sigma_{Q} . The boldfaces show the best performers for small/large particle sizes.

	${\sigma }_Q$
	${10}^{-2}$	$5\times {10}^{-2}$	${10}^{-1}$	$5\times {10}^{-1}$	$1$
SKF	41.4100	34.6611	29.9952	42.1360	38.0507
MKF	31.3088	27.6861	29.0376	35.2422	39.2536
$\alpha$KF	30.8783	27.9475	27.4130	31.0271	34.9420
PF	28.5429	32.3768	35.1842	44.3704	48.9767
EKF	33.8611	35.8086	37.7808	42.6595	45.9788
UKF	31.7528	31.8616	33.7625	41.1282	45.4806
BASE	223.5281	223.5281	223.5281	223.5281	223.5281

Table 2. TABLE II: Sensor network tracking problem: Mean Square Error (MSE) of various filtering schemes as a function of process noise parameter σ Q subscript 𝜎 𝑄 \sigma_{Q} . The boldfaces show the best performers for small/large particle sizes.

	${\sigma }_Q$
	${10}^{-2}$	$5\times {10}^{-2}$	${10}^{-1}$	$5\times {10}^{-1}$	$1$
SKF	10.4674	9.5812	9.5038	10.1664	10.5996
MKF	10.5572	9.2879	9.1684	9.8175	10.3307
$\alpha$KF	9.9441	8.0913	8.0623	9.1002	9.7055
PF	9.5661	9.3464	9.4726	10.0422	10.3834
EKF	14.0034	13.9357	14.5161	15.5277	16.1438
UKF	11.5303	10.3639	10.2068	10.8830	11.5845

Table 3. TABLE III: Mean Absolute Error(MAE) values of various filtering schemes for three different call/put option pairs; calculated for σ Q = 10 − 2 subscript 𝜎 𝑄 superscript 10 2 \sigma_{Q}=10^{-2} . For Option 3, EKF loses track so MAE is not reported.

		EKF	UKF	SKF	MKF
Option 1 MAE	Call	0.1352	0.0788	0.0658	0.0654
	Put	0.1528	0.0789	0.0642	0.0654
Option 2 MAE	Call	0.0425	0.0354	0.0312	0.0319
	Put	0.0478	0.0355	0.0368	0.0331
Option 3 MAE	Call	-	0.2155	0.1573	0.1586
	Put	-	0.2158	0.1574	0.1586

Equations

$$x_t=F_t{x}_{t-1}+w_t,y_t=H_t{x}_t+v_t,$$

$$p(x_t|y_{1:t-1})$$

$$p(x_t|y_{1:t})$$

$$x_t=F_t{x}_{t-1}+w_t,y_t=h(x_t)+v_t,$$

$$p(x_t|y_{1:t})=\frac{p(y_t|x_t)p(x_t|y_{1:t-1})}{\int p(y_t|x_t)p(x_t|y_{1:t-1})d{x}_t},$$

$$\displaystyle p(x_{t},y_{t})\sim N\left(\left[\begin{array}[]{c}\mu_{x}\\ \mu_{y}\end{array}\right],\left[\begin{array}[]{cc}\mathrm{\Sigma}_{xx}&\mathrm{\Sigma}_{xy}\\ \mathrm{\Sigma}_{yx}&\mathrm{\Sigma}_{yy}\end{array}\right]\right).$$

$$p(x_t|y_t)$$

$${\mu }_{x|y}$$

$${\mathrm{\Sigma }}_{x|y}$$

$$h(x_t)\approx h(x_0)+H(x_0)(x_t-x_0)$$

$$x^s$$

$$p(x_t|y_{1:t})\approx \sum _{s=1}^{N_s}{w}_t^s{\delta }_{x_t^s}.$$

$$w_t^s\propto \frac{1}{N_s}p(y_t|x_t^s).$$

$$x_{t+1}^s\stackrel{iid}{\sim }\sum _{s=1}^{N_s}{w}_t^{s}N(F_t{x}_t^s,Q_t),p(x_{t+1})\approx \sum _{s=1}^{N_s}\frac{1}{N_s}{\delta }_{x_{t+1}^s}$$

$$\text{KL}[q\Vert p]=\int q(x)\ln \frac{q(x)}{p(x|y)}dx.$$

$$\ln p(y)=\mathcal{L}(q,p(y,x))+\text{KL}[q(x)\Vert p(x|y)],$$

$$\mathcal{L}(q,p(y,x))$$

$$\mathcal{L}$$

$$f(x_t)$$

$$\displaystyle\nabla\mathcal{L}=\mathbb{E}_{q}[f(x_{t})\nabla\ln q(x_{t})]+\nabla\mathbb{E}_{q}\Big{[}\ln\frac{p(x_{t})}{q(x_{t})}\Big{]},$$

$${\mathbb{E}}_q[f(x_t)\nabla \ln q(x_t)]\approx \frac{1}{S}\sum _{s=1}^{S}f(x_t^s)\nabla \ln q(x_t^s),x_t^s\stackrel{iid}{\sim }q(x_t).$$

$$\nabla \mathcal{L}$$

$${\hat{\mu }}_t^{(i+1)}$$

$${\hat{\mathrm{\Sigma }}}_t^{(i+1)}$$

$$g(x_t)$$

$$\tilde{h}({\mu }_t,x_t)$$

$$g(x_t)=-\frac{1}{2}({\tilde{y}}_t-H({\mu }_t)x_t{)}^{\top }{R}^{-1}({\tilde{y}}_t-H({\mu }_t)x_t).$$

$${\nabla }_{{\hat{\mu }}_t}\mathcal{L}$$

$$+{\mathrm{\Sigma }}_t^{-1}({\mu }_t-{\hat{\mu }}_t)+H(x_0{)}^{\top }{R}^{-1}({\tilde{y}}_t-H(x_0){\hat{\mu }}_t),$$

$${\nabla }_{{\hat{\mathrm{\Sigma }}}_t}\mathcal{L}$$

$$\displaystyle\quad\frac{1}{2}\big{[}\widehat{\mathrm{\Sigma}}_{t}^{-1}(x^{s}-\widehat{\mu}_{t})(x^{s}-\widehat{\mu}_{t})^{\top}\widehat{\mathrm{\Sigma}}_{t}^{-1}-\widehat{\mathrm{\Sigma}}_{t}^{-1}\big{]}$$

$$+\frac{1}{2}({\hat{\mathrm{\Sigma }}}_t^{-1}-{\mathrm{\Sigma }}_t^{-1})-\frac{1}{2}H(x_0{)}^{\top }{R}^{-1}H(x_0).$$

$${\hat{\mu }}_t^{(i+1)}$$

$${\hat{\mathrm{\Sigma }}}_t^{(i+1)}$$

$$\text{KL}[p\Vert q]={\int }_{\mathcal{X}}p(x)\log \frac{p(x)}{q(x)}dx.$$

Full text

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 [1], 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 [2], collaborative filtering [3] and topic modeling [4].

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 [5], whereas dynamic collaborative filtering model contains a bilinear form in two unknown vectors [6]. The nonlinear problem has been studied extensively in the literature, resulting in well-known filtering algorithms such as the extended Kalman filter (EKF) [7] and unscented Kalman filter (UKF) [8]. On the other hand, Monte Carlo methods have been developed [9], 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 [10] and expectation-propagation [11], 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 [12, 13, 14]. 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.111We 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

II-A Basic Linear Framework

The Kalman filter [1] 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

[TABLE]

where $w_{t}$ and $v_{t}$ are independent zero-mean Gaussian random vectors with covariances $Q_{t}$ and $R_{t}$ respectively.222Kalman’s formulation in [1] is optimal for more general noise models, but Gaussian noise is the most common choice, which we also use in this paper. The latent variable $x_{t}\in\mathbb{R}^{d}$ is the unobserved state of the system. The vector $y_{t}\in\mathbb{R}^{p}$ constitutes the measurements made by the system.

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

[TABLE]

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

[TABLE]

where the noise process is the same as in Eq. (1), but $h(\cdot)$ is a nonlinear function of $x_{t}$.333We 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

[TABLE]

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 [15, 7, 8, 16, 17]. The main ingredient here is an assumption of joint Gaussianity of the latent and observed variables. This takes the form,

[TABLE]

(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

[TABLE]

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

[TABLE]

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.

As discussed in [8], the first-order approximation made by the EKF is often poor and performance can suffer as a result. Instead, they propose to estimate the quantities in (9) with an unscented transform—a numerical quadrature method. The result is the UKF, which has similar computational cost as the EKF and higher accuracy. Based on the calculated Gaussian prior $p(x_{t}|y_{1:t-1})=N(x_{t}|\mu_{x},\mathrm{\Sigma}_{xx})$, the UKF selects a discrete set of sigma points at which to approximate $\mu_{y}$, $\mathrm{\Sigma}_{yy}$ and $\mathrm{\Sigma}_{xy}$. Let $d_{x}=\mathrm{dim}(x_{t})$ and $N_{s}=2\times d_{x}+1$. These sigma points $x^{1},\dots,x^{N_{s}}$ are

[TABLE]

The vector $[\sqrt{(d_{x}+\lambda)\mathrm{\Sigma}}]_{s}$ corresponds to the $s$th column of the Cholesky decomposition of the matrix $\sqrt{(d_{x}+\lambda)\mathrm{\Sigma}}$. Positive weights $w^{s}$ are also defined for each $x^{s}$. The constant $\lambda$ controls these sigma point locations, as well as the weights (along with additional fixed parameters). These $N_{s}$ locations are used to empirically approximate all means and covariances in Eq. (9). Once $y_{t}$ is measured, the approximation of $p(x_{t}|y_{t})$ can then be calculated using Eq. (10).

There are many extensions to the UKF framework such as cubature Kalman filtering (CKF) [18] and QMC Kalman filtering [17], 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 [11] use moment matching (e.g., [19]) 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

We have seen that the main theme of ADF is approximating the posterior with a pre-specified joint probability density; when this joint density is Gaussian then $p(x_{t}|y_{t})\approx N(\mu_{t},\mathrm{\Sigma}_{t})$. On the other hand, nonparametric versions use sampling for posterior approximation without making any density assumptions on the form of this posterior,

[TABLE]

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) [9], 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

[TABLE]

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

[TABLE]

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 [14].

III Three divergence minimization approaches

In this section we discuss the three proposed divergence minimization approaches to the nonlinear Kalman filtering problem. These include the two directions of the Kullback-Leibler (KL) divergence as well as the related $\alpha$-divergence that contains both KL divergence measures as limiting cases. In all cases, our goal is to approximate the intractable posterior distribution $p(x_{t}|y_{t})$ with a multivariate Gaussian distribution $q(x_{t})=N(\mu_{t},\mathrm{\Sigma}_{t})$, using these three divergences as potential quality measures. We again note that our contribution is to provide three unbiased methods for directly optimizing these divergences without introducing additional approximations. We therefore anticipate an improvement over the standard EKF and UKF approximations. In each of the following three subsections, we present the divergence objective, review its tractability issues, and then present our approach to resolve this issue.

III-A Approach 1: Forward KL divergence minimization

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

[TABLE]

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) [22]. 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 [23, 22] instead uses the identity

[TABLE]

where

[TABLE]

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.

While nonlinear Kalman filters have a simply-defined joint likelihood $p(y_{t},x_{t}|y_{1:t-1})$ at time $t$, a significant problem still arises in calculating $\mathcal{L}$ due to the nonlinear function $h$. That is, if we define $q(x_{t})=N(\mu_{t},\mathrm{\Sigma}_{t})$, then for the Gaussian generative process of Eq. (3) we optimize $\mu_{t}$ and $\mathrm{\Sigma}_{t}$ over the function

[TABLE]

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. [24] 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,

[TABLE]

over the parameters of $q(x_{t})=N(x_{t}|\mu_{t},\mathrm{\Sigma}_{t})$, and thereby minimize $\text{KL}[q\|p]$. This can be done by gradient ascent. However, since $\mathbb{E}_{q}[f(x_{t})]$ does not have a closed form solution, $\nabla\mathcal{L}$ can not be evaluated analytically. The proposed solution in [24] is to instead step in the direction of an unbiased stochastic gradient. To this end, the observation is made that

[TABLE]

where the identity $\nabla q(x_{t})=q(x_{t})\nabla\ln q(x_{t})$ is used. While the second gradient can be calculated analytically with respect to either $\mu_{t}$ or $\mathrm{\Sigma}_{t}$, the first gradient can be sampled using Monte Carlo integration,

[TABLE]

A second observation is made by [24] that the variance of these samples may be so large that $S$ needs to be set to too large a number to make this approximation computationally feasible. For this reason employing variance reduction methods is crucial. Paisley, et al. [24] propose introducing a control variate $g(x_{t})$ that is highly correlated with $f(x_{t})$, but has an analytic expectation $\mathbb{E}_{q}[g(x_{t})]$. The gradient of $\mathcal{L}$ with a control variate is equal to

[TABLE]

Though this leaves the gradient unchanged, MC sampling of the first term has much smaller variance when $|corr(f,g)|$ is large (calculated using $q(x_{t})$). The parameter $\lambda\in\mathbb{R}$ is set to minimize the variance.444As shown in [24], when $\lambda\equiv cov(f,g)/var(g)$ (approximated by sampling) the variance reduction ratio is $var(f-\lambda g)/var(f)=1-corr(f,g)^{2}$. Intuitively, this can be seen by noting that if $f(x_{t}^{s})\approx\lambda g(x_{t}^{s})$ at the sampled values $x_{t}^{s}$, then $|f(x_{t}^{s})-\lambda g(x_{t}^{s})|\ll|f(x_{t}^{s})|$. In this case, the analytic gradient $\lambda\nabla\mathbb{E}_{q}[g(x_{t})]$ gives an initial approximation of $\mathbb{E}_{q}[f(x_{t})\nabla\ln q(x_{t})]$, which is then corrected to be made unbiased by the MC-sampled $\mathbb{E}_{q}[(f(x_{t})-\lambda g(x_{t}))\nabla\ln q(x_{t})]$. Since $g(x_{t})$ is a good approximation of $f(x_{t})$ in the region of high probability defined by $q(x_{t})$, the analytic approximation captures most information, but is refined by the MC-sampled gradient to make the method approximation-free.

The requirements on $g(x_{t})$ to be a good control variate for $f(x_{t})$ are that: 1) it is an approximation of $f(x_{t})$, and 2) the expectation $\mathbb{E}_{q}[g(x_{t})]$ is solvable. There are many possible control variates for the function $(y-h(x))^{\top}R^{-1}(y-h(x))$. However, building on the EKF framework we propose setting

[TABLE]

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

[TABLE]

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.

Next we derive the unbiased gradients. In this case, these gradients are $\nabla_{\widehat{\mu}_{t}}\mathcal{L}$ and $\nabla_{\widehat{\mathrm{\Sigma}}_{t}}\mathcal{L}$. We note that, if we were using the EKF framework, by replacing $f(x_{t})$ with $g(x_{t})$ in Eq. (19), the roots of these gradients could be solved and the EKF solutions for $\widehat{\mu}_{t}$ and $\widehat{\mathrm{\Sigma}}_{t}$ would result. However, since we have the additional stochastic gradient term, we must perform gradient ascent. The final expressions for the unbiased gradients using samples $x_{t}^{s}\sim_{iid}q(x_{t})$ are:

[TABLE]

Comparing the gradients (25)-(26) we see that stochastic search acts as a correction step to the EKF updates resulting from the last lines alone. On the other hand, unlike EKF, we cannot simply solve for $\widehat{\mu}_{t}$ and $\widehat{\mathrm{\Sigma}}_{t}$ by settings the gradients equal to zero, so we use gradient ascent. Without proper scaling we can easily have a numerically unstable algorithm and the covariance matrix can lose its positive definiteness. To fix this we pre-condition the gradients with a symmetric positive definite matrix $C$ and perform the following updates

[TABLE]

We note the difference between index $t$ and $i$. The first is the time frame we are currently processing, while the second is the iteration number at time $t$ since we are using a gradient optimization method. For the conditioning matrix we choose $C^{(i)}=[\mathrm{\Sigma}_{t}^{(i)}]^{-1}$. Using this setting, we get approximate natural gradients [25] for $\widehat{\mu}$ and $\widehat{\mathrm{\Sigma}}$. When the step size satisfies the Robbins-Monro conditions, $\sum_{i=1}^{\infty}\rho^{(i)}=\infty$ and $\sum_{i=1}^{\infty}[\rho^{(i)}]^{2}<\infty$, the gradients in (27)-(28) converge to a fixed point of the exact variational lower bound. In practice we can, for example, choose $\rho^{(i)}=(w+i)^{-\eta}$ with $\eta\in(0.5,1]$ and $w\geq 0$. In simulations we observed that, when natural gradients are employed a generic schedule for step sizes can be used, and no further hand-tuning is necessary. We refer to this algorithm as stochastic search Kalman filtering (SKF) and summarize it in Algorithm 1 for a single time step.

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:

[TABLE]

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 [23]. 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

[TABLE]

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

[TABLE]

which follows from the exponential family identity $\nabla_{\eta}\log A(\eta)=\mathbb{E}_{q}[s(x)]$. Therefore the fixed points of the objective are given by

[TABLE]

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

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 [22]. 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)$.

Therefore, in the context of exponential families the task of finding the optimal distribution for the reverse KL reduces to the task of matching moments. However, there is still a difficulty in the need to compute the moments of an unknown posterior distribution. Fortunately, Monte Carlo methods prove useful here as well. Let $I(f)=\mathbb{E}_{q}(f(x))$ be the expectation we wish to compute. For example, choosing $f(x)=x$ and $f(x)=xx^{\top}-\mathbb{E}[x]\mathbb{E}[x]^{\top}$ gives the mean and covariance respectively. This expectation can be approximated as

[TABLE]

We will define

[TABLE]

and so $\mathbb{E}_{q}[f(x)]\approx\frac{1}{W}\sum_{s=1}^{S}f(x^{s})w^{s}$. This is related to importance sampling, with the added normalizer $W$. As we can see from (31) this procedure is biased as it is a ratio of two approximations, yet it converges to the true expectation $\mathbb{E}_{q}[f(x)]$ almost surely. Therefore, we have an asymptotically unbiased divergence minimization procedure. We call this the moment matching Kalman filter (MKF) and summarize it in Algorithm 2 for Gaussian distributions, as Gaussian approximations are our focus in this paper. We observe that a major difference between the MKF and SKF of the previous section is that the MKF only needs to sample once to obtain the moment estimates for a time step. Therefore, the MKF is not an iterative algorithm and is much faster. Also, the MKF is slightly faster than particle filtering as it eliminates the need for resampling.

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 [27] we define the $\alpha$-divergence to be

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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 [29]. More recently, [27] 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

[TABLE]

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

[TABLE]

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.

Since the approximate Gaussian posterior distribution has a parametric form, we are able to use an adaptive sampling method for the MKF and $\alpha$KF .555For the SKF the per-iteration sample size is much smaller, so there is less benefit in using this technique. To determine the appropriate number of samples, we measure the uncertainty of our mean approximation for $q(x_{t})=\mathcal{N}(\mu_{t},\mathrm{\Sigma}_{t})$, where $\mu_{t}=\sum_{s=1}^{S}w^{s}x^{s}/W$. For importance sampling, the variance of this estimator is approximately

[TABLE]

Here, if $S$ is large enough the estimator can be approximated as normal by the central limit theorem [20]. We use this to compute the radius of a 95% confidence region. Without loss of generality assume that the estimator is zero mean, which is justified by the asymptotic unbiasedness of the unnormalized importance sampling procedure. We denote the estimator by $\widehat{X}\sim N(0,\mathbb{V}(\mu_{t}))$. Then we have

[TABLE]

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

[TABLE]

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

[TABLE]

As expected, the smaller radius we desire, the larger sample size we need. We note in passing that $\mathbb{V}(\mu_{t})$ is our confidence in estimating the mean of the true posterior, and not the ground truth. The accuracy of estimating the latter is dictated by the measurement noise, and cannot be made arbitrarily small by increasing the sample size.

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.

IV-A Target Tracking

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 [30, 31].

For both problems the state-space will have the form

[TABLE]

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) [32]. 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}$ [33].

Our experiments are based on synthetic data using a constant velocity model in $\mathbb{R}^{2}$ which corresponds to the state vector vector $x_{t}=[x_{1},~{}\dot{x_{1}},~{}x_{2},~{}\dot{x_{2}}]^{\top}$; the second and fourth entries correspond to the velocity of the target in each dimension. Following [34], we set the parameters for the state variable equation to

[TABLE]

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

[TABLE]

i.e. the Cartesian-to-polar transformation [5]. 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

[TABLE]

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 [33]. 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 [35]

[TABLE]

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

[TABLE]

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 [37] 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.