Maximum Correntropy Ensemble Kalman Filter

Throughout the article, we use boldface lower-case letters for vectors and boldface upper-case letters for matrices. The transpose and mathematical expectation are denoted by $\{\cdot\}^{\top}$ and $\mathbb{E}\left[\cdot\right]$, respectively. The Gaussian distribution with mean $\mu$ and covariance $\bm{\Sigma}$ is denoted by $\mathcal{N}(\mu,\bm{\Sigma})$. $\mathbb{R}^{n}$ and $\mathbb{R}^{m\times n}$ denote, respectively, the $n$-dimensional Euclidean space and the set of all $n\times m$ real matrices. Especially, $\mathbf{I}_{n}$ denotes the n-dimensional identity matrix. weighted $l^{2}$ norm $\|\mathbf{x}\|_{\mathbf{A}}^{2}:=\mathbf{x}^{\top}\mathbf{A}^{-1}\mathbf{x}$ with $\mathbf{x}\in\mathbb{R}^{n}$ and non-singular matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$.

In this article, we consider the nonlinear state-space model given by the following state and observation equations:

where state $\mathbf{x}_{k}\in\mathbb{R}^{n}$, and observation $\mathbf{y}_{k}\in\mathbb{R}^{m}$. $\mathbf{f}_{k}$ and $\mathbf{h}_{k}$ are nonlinear functions called state equation and observation equation, respectively. State noise $\mathbf{w}_{k}$ and observation noise $\mathbf{v}_{k}$ are zero means with nominal covariance $\mathbf{Q}_{k}\in\mathbb{R}^{n\times n}$ and $\mathbf{R}_{k}\in\mathbb{R}^{m\times m}$, respectively. Let $\mathbf{y}_{1:k}$ denote the $\sigma$-algebra generated by noisy observations $\{\mathbf{y}_{1},\dots,\mathbf{y}_{k}\}$. The filtering problem refers to estimating the posterior distribution $p(\mathbf{x}_{k}\mid\mathbf{y}_{1:k})$, which is called filtering distribution.

In the follow-up, we focus on the scenario in which the observation noises are not Gaussian, i.e., they are induced by unknown large outliers. Hence, the real distribution of observation noise is unknown to us in fact. For example, we consider $(1-\epsilon)\ \mathcal{N}(0,\mathbf{R}_{k})+\epsilon\ \mathcal{S}(0,\mathbf{S}_{k})$ , where $0<\epsilon\ll 1$ is the unknown probability, and $\mathcal{S}(0,\mathbf{S}_{k})$ is an arbitrary unknown distribution with large covariance $\mathbf{S}_{k}$. Here $\mathbf{R}_{k}$ is known to us, so it is called the nominal covariance matrix. Additionally, in what follows, we shall develop our novel filter to deal with such noises above by utilizing the MCC as the cost function involving the nominal observation covariance $\mathbf{R}_{k}$.

The EnKF sequentially approximates the filtering distributions $p(\mathbf{x}_{k}\mid\mathbf{y}_{1:k})$ using $N$ equally weighted ensembles $\{\mathbf{x}_{k}^{(1)},\dots,\mathbf{x}_{k}^{(N)}\}$. At prediction steps, each ensemble $\mathbf{x}_{k}^{(i)}$ is propagated using the state equation, while at update steps a Kalman-type update is performed for each ensemble:

and

where $\mathbf{w}_{k}^{(i)}\sim\mathcal{N}(0,\mathbf{Q}_{k})$, $\mathbf{v}_{k}^{(i)}\sim\mathcal{N}(0,\mathbf{R}_{k})$ and the Kalman gain

is defined using the empirical prediction covariance $\hat{\mathbf{C}}_{k}$ of prediction ensembles $\{\mathbf{x}_{k\mid k-1}^{(1)},\dots,\mathbf{x}_{k\mid k-1}^{(N)}\}$, namely

with

and

Here we present how to derive MC-EnKF based on the MCC estimate (3). The prediction step of MC-EnKF is the same as those of EnKF, i.e., (6). Therefore, we shall focus on deriving its update step. Let $\hat{\mathbf{x}}_{k\mid k-1}$ and $\mathbf{P}_{k\mid k-1}$ denote the prediction mean $\mathbb{E}\left[\mathbf{x}_{k}\mid\mathbf{y}_{1:k-1}\right]$ and prediction covariance $\mathbb{E}\left[\left(\mathbf{x}_{k}-\hat{\mathbf{x}}_{k\mid k-1}\right)\left(\mathbf{x}_{k}-\hat{\mathbf{x}}_{k\mid k-1}\right)^{\top}\mid\mathbf{y}_{1:k-1}\right]$, respectively. As shown in [22], for the underlying linear system of (5), i.e., $\mathbf{f}_{k}(\mathbf{x}_{k})=\mathbf{F}_{k}\mathbf{x}_{k}$ and $\mathbf{h}_{k}(\mathbf{x}_{k})=\mathbf{H}_{k}\mathbf{x}_{k}$, the update step of KF can be derived by using least square cost function,

Motivated by (13), we shall consider a new cost function based on MCC. Besides, in view of Remark II.1, we shall replace $\hat{\mathbf{x}}_{k\mid k-1}$ and $\mathbf{P}_{k\mid k-1}$ in (13) by the empirical prediction mean $\hat{\mathbf{m}}_{k}$ in (10) and the empirical prediction covariance $\hat{\mathbf{C}}_{k}$ in (9). Therefore, our modified cost function for ensembles $\mathbf{x}_{k\mid k}^{(i)}$ for $i=1,2,\dots,N$ is given by

Then based on $\mathcal{L}_{2}(\mathbf{x}_{k})$, the update step for MC-EnKF can be obtained solving the following optimization problem for $i=1,2,\dots,N$:

In what follows, the derivation contains some approximations. For the sake of obtaining our algorithm, we heuristically treat them as exact equalities. Let us denote

Then recall $\mathbf{H}_{k}$ defined in (11), we consider this approximation when taking the derivative,

Now letting $\frac{\partial\mathcal{L}_{2}\left(\mathbf{x}_{k}\right)}{\partial\mathbf{x}_{k}}=0$, we have

Adopting the first-order Taylor series to approximate the nonlinear observation function $\mathbf{h}_{k}$ at $\mathbf{x}_{k\mid k-1}^{(i)}$, i.e.,

Substituting (19) into (18), we have

Then we obtain the following stochastic ensemble update rule like (7) for $\mathbf{x}_{k\mid k}^{(i)}$ with drawing $\mathbf{v}_{k}^{(i)}\sim\mathcal{N}(0,\mathbf{R}_{k})$:

where the new Kalman gain $\tilde{\mathbf{K}}_{k}$ is given by

Here we shall discuss the practical algorithm, since $l_{\mathbf{R}_{k}}$ and $l_{\hat{\mathbf{C}}_{k}}$ all contain the item $\mathbf{x}_{k}$. Here we introduce a novel technical approximation. Recall (16), we use $\hat{\mathbf{m}}_{k}$ to approximate $\mathbf{x}_{k}$ contained in $l_{\mathbf{R}_{k}}$ and $l_{\hat{\mathbf{C}}_{k}}$, respectively, i.e.,

It follows that we summarise the specific algorithm steps of MC-EnKF in Algorithm 1.

In general, the kernel bandwidth (scale parameter) in the cost function balances the convergence rates of the regression model and its robustness [23]. Since our algorithm does not involve solving such regression problems, we focus on the how to tune this scale parameter for robustness (with respect to outliers). Intuitively, if the arrived observation $\mathbf{y}_{k}$ contains large outliers, $\|\mathbf{y}_{k}-\mathbf{h}_{k}(\hat{\mathbf{x}}_{k})\|_{2}$ will be large, where $\|\cdot\|_{2}$ denote the $l^{2}$-norm. In this case, we need a smaller $\sigma$ to make our algorithm more robust. Motivated by this intuition, we propose to set $\sigma_{k}=\frac{1}{\|\mathbf{y}_{k}-\mathbf{h}_{k}(\hat{\mathbf{x}}_{k})\|_{2}}$ adaptively, where $\sigma_{k}$ denote the kernel bandwidth used in the $k$-th step for MC-EnKF. As shown in our simulation experiments, this adaptive strategy can achieve a good estimation performance. Additionally, the MC-EnKF will act more and more like the common EnKF algorithm as $\sigma$ increases. In particular, the following convergence theorem holds .

The first benchmark we used here is a linear system, which is given by

where $\alpha=\frac{\pi}{18}$, $\mathbf{x}_{0}\sim\mathcal{N}(0,\mathbf{I}_{2})$ and $\mathbf{w}_{k}\sim\mathcal{N}(0,0.01\mathbf{I}_{2})$. Let $\mathbf{Q}_{1}=0.01\mathbf{I}_{1}$ be the nominal observation covariance, the observation noises are sampled from the mixture of Gaussian, i.e.,

Table I lists the MSEs and average CPU times of EnKF and MC-EnKF with different kernel bandwidths in this example, where the number of ensembles of EnKF and MC-EnKF are set as $100$. Note that the MC-EnKF outperforms EnKF in most cases and their average CPU times are virtually identical to the common EnKF. This demonstrates that our proposed MC-EnKF has better performance than those of EnKF with nearly no extra cost on CPU time. It means that MC-EnKF is more robust and efficient when we deal with observation noises with large outliers. In fact, our proposed adaptive strategy (MC-EnKF-Ada) outperforms most different $\sigma$ cases and is only slightly worse than the best case. This suggests that our adaptive strategy works well for linear system. When we conduct the simulation for large enough $\sigma$, for example, $\sigma=10000$, the result verifies Theorem III.1 for linear system since we note that the performance of MC-EnKF is almost identical to the EnKF in this case.

In order to better present the comparison results, we choose the case where $\sigma=5$ to illustrate the true state and EnKF estimation and MC-EnKF estimation over time, which is shown in Fig 1. Here, the selected dimension is indicated by the - suffix of the legend in the figures; for example, EnKF-1 denotes the EnKF estimation on the first state dimension. This setting will also be used for the following figure. We can clearly see that this MC-EnKF gives accurate estimates but EnKF performs poorly when the observations contain large outliers.

The second benchmark is a nonlinear system, which is given as follows:

where $\mathbf{x}_{0}\sim\mathcal{N}(0,\mathbf{I}_{2})$ and $\mathbf{w}_{k}\sim\mathcal{N}(0,\mathbf{I}_{2})$. $\kappa_{1}=\kappa_{2}=0.1$ are constants controlling the state dynamics. Let $\mathbf{Q}_{2}=\mathbf{I}_{2}$ be the nominal observation covariance, the observation noises are sampled from the mixture of Gaussian, i.e.,

We list the MSEs and average CPU times of EnKF, ML-EnKF and MC-EnKF with different kernel bandwidths in Table II for this nonlinear example, where the number of ensembles of these filters still be set as $100$. It is obvious that for nonlinear system the EnKF degrades greatly in the presence of non-Gaussian observation noises and the linear approximation error of observation function, hence has the worst estimate performance. ML-EnKF performs better than EnKF, but is still affected by outliers. Better results can be obtained with our proposed adaptive strategy (MC-EnKF-Ada) than with some kernel bandwidths, indicating that it is still useful for nonlinear systems. Note that the MC-EnKF outperforms EnKF in all cases, which demonstrate that our proposed MC-EnKF can effectively eliminate the effect of the non-Gaussian observation noises. Additionally, we notice that they share almost identical CPU times, which demonstrates the efficiency of MC-EnKF. Moreover, with suitable kernel bandwidth ($\sigma=5$), MC-EnKF can achieve very accurate estimates while EnKF is heavily affected by large outliers, which is presented in Fig 2. We also note that when $\sigma$ is large enough ($\sigma=10000$), the performance of MC-EnKF is identical to the EnKF, which verifies the result of Theorem III.1 for nonlinear system.

According to the simulation experiments mentioned above, we notice the following two arguments:

• •

  With suitable kernel bandwidth, MC-EnKF can significantly improve its robustness to observation noises containing large outliers with nearly no CPU time loss.

• •

  With large kernel bandwidth, MC-EnKF performs like the common EnKF, which implies its flexibility to handle Gaussian noise cases.

These two arguments make our viewpoint a very fruitful area for further study. Besides, we note that the kernel bandwidth $\sigma$ plays a significant role in the algorithm and will have an impact on its robustness to non-Gaussian observation (such as outliers). Although our proposed adaptive strategy is effective, it cannot achieve the performance with suitable kernel bandwidth. Therefore, a highly intriguing future work direction may be a more effective adaptive strategy to pick the right kernel bandwidth with theoretical guarantee since it is crucial for our proposed MC-EnKF. In Addition, this article only investigates using the special case (MCC) of generalized MCC [26] as a cost function to derive a robust EnKF with stochastic update steps. Hence how to develop robust variants of EnKF [9] based on the generalized MCC will be our next research focus.