# EMORF/S: EM-Based Outlier-Robust Filtering and Smoothing With Correlated
Measurement Noise

Aamir Hussain Chughtai, Muhammad Tahir, , and Momin Uppal The authors are with
Department of Electrical Engineering, Lahore University of Management
Sciences, DHA Lahore Cantt., 54792, Lahore Pakistan. (email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this article, we consider the problem of outlier-robust state estimation
where the measurement noise can be correlated. Outliers in data arise due to
many reasons like sensor malfunctioning, environmental behaviors,
communication glitches, etc. Moreover, noise correlation emerges in several
real-world applications e.g. sensor networks, radar data, GPS-based systems,
etc. We consider these effects in system modeling which is subsequently used
for inference. We employ the Expectation-Maximization (EM) framework to derive
both outlier-resilient filtering and smoothing methods, suitable for online
and offline estimation respectively. The standard Gaussian filtering and the
Gaussian Rauch–Tung–Striebel (RTS) smoothing results are leveraged to devise
the estimators. In addition, Bayesian Cramer-Rao Bounds (BCRBs) for a filter
and a smoother which can perfectly detect and reject outliers are presented.
These serve as useful theoretical benchmarks to gauge the error performance of
different estimators. Lastly, different numerical experiments, for an
illustrative target tracking application, are carried out that indicate
performance gains compared to similarly engineered state-of-the-art outlier-
rejecting state estimators. The advantages are in terms of simpler
implementation, enhanced estimation quality, and competitive computational
performance.

###### Index Terms:

State-Space Models, Approximate Bayesian Inference, Nonlinear Filtering and
Smoothing, Outliers, Kalman Filters, Variational Inference, Expectation-
Maximization, Robust Estimation, Statistical Learning, Stochastic Dynamical
Systems.

##  I Introduction

State estimation is a key fundamental task in analyzing different dynamical
systems with subsequent decision-making and control actions arising in a
variety of fields including cybernetics, robotics, power systems, sensor
fusion, positioning, and target tracking [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] etc.
The states describing the system dynamics can evolve intricately. Moreover,
these are not directly observable only manifesting themselves in the form of
external measurements. Mathematically, it means that for state estimation in
general, the inference is performed considering stochastic nonlinear equations
making it a nontrivial task.

Filtering is the common term used for online state estimation where inference
is carried out at each arriving sample. Kalman filter with its linear and
nonlinear versions [11, 12, 13] are considered the primary choices for
filtering given their ease of implementation and estimation performance. Other
options for nonlinear filtering are also available including methods based on
Monte-Carlo (MC) approximations e.g. Particle Filters (PFs) [14], ensemble
Kalman filter (EnKF) [15] etc.

Smoothing, on the other hand, refers to offline state estimation where the
primary concern is not to work on a per-sample basis. We are rather interested
in state inference considering the entire batch of measurements. Different
options for smoothing exist including the famous Rauch–Tung–Striebel (RTS),
two-filter smoothers [16, 17] etc.

The standard state estimators are devised with the assumption that the
dynamical system under consideration is perfectly modeled. The estimators
assume the availability of system and observation mathematical models
including the process and measurement noise statistics. However, any modeling
mismatch can result in deteriorated performance even possibly crippling the
functionality of the regular estimators completely.

In this work, we are interested in coping with the modeling discrepancy and
the associated estimation degradation that results from the occurrence of
outliers in the measurements. Data outliers can arise due to several factors
including data communication problems, environmental variations, and effects,
data preprocessing front-end malfunctioning, inherent sensor defects, and
degradation, etc [9]. We keep our consideration generic by taking into account
the possibility of correlated measurement noise with a fully enumerated
nominal noise covariance matrix. This is in contrast to the existing
approaches where noise in each data dimension is assumed to be independent
targeting a specific class of applications [18]. However, this leaves out
several important application scenarios where noise correlation exists which
should be taken into account. For example, due to double differencing of the
original measurements in Real Time Kinematic (RTK) systems, noise correlation
appears [19]. Likewise, a significant negative correlation exists between the
range and range rate measurement noise in radar data [20]. Similarly, due to
the use of a common reference sensor to extract the time difference of arrival
(TDOA), correlated range measurement noise arises [21]. Besides, in different
sensor networks, correlated observation noise also emerges [22, 23, 24].

(a) Static loss functions in traditional methods

(b) Adaptive loss functions in learning-based methods

Figure 1: Typical loss functions for outlier-robust state estimators. In
traditional approaches, the loss function is static. In learning-based
methods, the loss function adapts e.g. between a quadratic function and a
constant to weight the data during inference.

The problem of neutralizing data outliers during state estimation has been
approached with various proposals. The traditional way of dealing with
outliers is based on assuming fixed statistics for measurement noise or the
residuals between predicted and actual measurements. For example, different
methods resort to describing observation noise using heavy-tailed
distributions like the Student-t and Laplace densities [25, 26]. Similarly,
the theory of robust statistics suggests the use of prior models for residuals
to downweigh the effect of outliers during inference [27, 28, 29]. Moreover,
some techniques are based on rejecting the data sample by comparing the
normalized measurement residuals with some predefined thresholds [30, 31].

Literature survey indicates that performances of the conventional approaches
are sensitive to the tuning of the design parameters which affect the static
residual error loss functions during estimation [32]. Therefore, tuning-free
learning-based techniques have been justified in prior works that make the
error loss function adaptive [32, 33, 34, 18, 19, 35]. These approaches
consider appropriate distributions for the measurement noise and subsequently
learn the parameters describing the distributions and consequently the loss
functions during state estimation. Fig. 1 depicts the comparison of typical
static loss functions in traditional approaches and dynamic loss functions in
learning-based methods considering a uni-dimensional model for visualization
(see Section III-D [32] for more details). Resultingly, learning-based robust
state estimators offer more advantages by reducing user input, being more
general, and suiting better for one-shot scenarios.

Several learning-based methods for robust state estimation have been reported
in the relevant literature. As exact inference is not viable for developing
these approaches, approximate inference techniques like PFs and variational
Bayesian (VB) methods can be used in their design. Since PFs can be
computationally prohibitive, VB-based techniques are the appealing alternative
considering these can leverage the existing standard filtering and smoothing
results. Our focus in this work remains on the learning-based outlier
mitigation approaches designed using VB.

In previous works, we observe that various outlier-robust state estimators,
devised using VB, treat the entire measurement vector collectively during
estimation owing to under-parameterized modeling [32, 33, 34]. Instead of
treating each dimension individually, the complete vector is either considered
or downweighted by varying the noise covariance matrix by a scalar
multiplicative factor. This leaves room for improvement considering useful
information is unnecessarily lost during inference. In this regard, we offer a
vectorial parameterization to treat each dimension individually in [18].
Therein we also suggest a way to make the estimators in [32, 33, 34]
selective. However, these proposals are based on the assumption of independent
noise for each measurement dimension. Another learning-based outlier-resilient
filter has been presented in [35]. But the authors only consider linear
systems and test the method with diagonal measurement noise covariance
matrices. With the possibility of correlated measurements, the Variational
Bayes Kalman Filter (VBKF) has been devised [19] by extending the work in
[33]. However, we observe that VBKF assumes a complex hierarchical model. As a
result, along with updating the state densities, it involves updating the
nuisance parametric distributions and their hierarchical distributions during
the VB updates. This includes evaluation of the digamma function to find the
expectation of logarithmic expressions [33, 19]. Therefore, implementing VBKF
can get complicated e.g. within an embedded computing device where access to
such functions is not inherently available and additional libraries are
required. Moreover, extending VBKF to outlier-robust smoothing also gets
cumbersome. This calls for simpler state estimation approaches, for systems
with correlated noise, that can weather the effect of outliers.

With this background, considering the possibility of correlated measurements
in nonlinear dynamical systems, we make the following contributions in this
work.

  * •

Using a suitable model and VB (more specifically Expectation-Maximization
(EM)) we devise an outlier-robust filter availing the standard Gaussian
filtering results. The results are further utilized in deriving an outlier-
robust smoother based on the standard Gaussian RTS smoothing. Since our
proposed method is inspired by our prior work which considers independent
measurement noise [18], we also present insightful connections.

  * •

We derive Bayesian Cramer-Rao Bounds (BCRBs) for a filter and a smoother which
can perfectly detect and reject outliers. This provides a useful benchmark to
assess the estimation ability of different outlier-mitigating estimators.

  * •

We evaluate the performance of the devised estimators as compared to the other
similarly devised outlier-discarding methods. Different scenarios of a
relevant TDOA-based target tracking application are considered in numerical
experiments indicating the merits of the proposed methods.

The rest of the article is organized as follows. Section II provides the
modeling details. In Section III, we present the derivation of the proposed
filter. Thereafter, the derivation of the proposed smoother is given in
Section IV. In Section V, BCRBs for a filter and a smoother with perfect
outlier detecting and rejecting capabilities are provided. Subsequently, the
performance evaluation results have been discussed in Section VI. The paper
ends with a conclusive commentary in Section VII.

Notation: As a general notation in this work, $\mathbf{r}^{\top}$ is the
transpose of the vector $\mathbf{r}$, ${{r}^{i}}$ denotes the $i$th element of
a vector $\mathbf{r}$; ${\mathbf{r}^{i-}}$ is the vector $\mathbf{r}$ with its
$i$th element removed; the subscript $k$ is used for time index;
$\mathbf{r}_{k{}}$ is the vector $\mathbf{r}$ at time instant $k$;
$\mathbf{r}_{k{-}}$ is the group of vectors $\mathbf{r}$ considering the
entire time horizon except the time instant $k$; ${{R}^{i,j}}$ is the element
of the matrix $\mathbf{R}$ present at the $i$th row and $j$th column;
$\mathbf{R}^{-1}$ is the inverse of $\mathbf{R}$; $|\mathbf{R}|$ is the
determinant of $\mathbf{R}$; $\bm{\mathfrak{R}}$ is the swapped form of
$\mathbf{R}$ where the swapping operation is defined in a particular context;
$\delta(.)$ represents the delta function;
$\langle.\rangle_{q(\bm{\psi}_{k})}$ denotes the expectation of the argument
with respect to the distribution $q(\bm{\psi}_{k})$; $\mathrm{tr}(.)$ is the
trace operator; $a\mod b$ denotes the remainder of $a/b$; the superscripts $-$
and $+$ are used for the predicted and updated filtering parameters
respectively; the superscript $s$ is used for the parameters of the marginal
smoothing densities. Other symbols are defined in their first usage context.

##  II State-space modeling

###  II-A Standard modeling

Consider a standard nonlinear discrete-time state-space model (SSM) to
represent the dynamics of a physical system given as

$\displaystyle\mathbf{x}_{k}$
$\displaystyle=\mathbf{f}(\mathbf{x}_{k-1})+\mathbf{q}_{k-1}$ (1)
$\displaystyle\mathbf{y}_{k}$
$\displaystyle=\mathbf{h}(\mathbf{x}_{k})+\mathbf{r}_{k}$ (2)

where $\mathbf{x}_{k}\in\mathbb{R}^{n}$ and $\textbf{y}_{k}\in\mathbb{R}^{m}$
denote the state and measurement vectors respectively; the nonlinear functions
$\textbf{f}(.):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and
$\textbf{h}(.):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ represent the process
dynamics and observation transformations respectively;
$\textbf{q}_{k}\in\mathbb{R}^{n}$ and $\textbf{r}_{k}\in\mathbb{R}^{m}$
account for the additive nominal process and measurement noise respectively.
$\textbf{q}_{k}$ and $\textbf{r}_{k}$ are assumed to be statistically
independent, White, and normally distributed with zero mean and known
covariance matrices $\textbf{Q}_{k}$ and $\textbf{R}_{k}$ respectively. We
consider that $\textbf{R}_{k}$ can be a fully enumerated matrix capturing the
correlations between the measurement noise entries.

###  II-B Modeling outliers for inference

The model in (1)-(2) assumes that the measurements are only affected by
nominal measurement noise $\textbf{r}_{k}$. However, the observations in every
dimension can be corrupted with outliers leading to the disruption of standard
state estimators as the measurement data cannot be described by the regular
model. Therefore, data outliers need to be appropriately modeled within the
generative SSM with two basic objectives. Firstly, the model should
sufficiently capture the effect of outlier contamination in the data.
Secondly, the model should remain amenable to inference.

To model the outliers in SSM, we consider an indicator vector
$\bm{\mathcal{I}}_{k}\in\mathbb{R}^{m}$ having Bernoulli elements where its
$i$th element ${{\mathcal{I}}}^{i}_{k}$ can assume two possible values:
$\epsilon$ (close to zero) and 1. ${{\mathcal{I}}}^{i}_{k}=\epsilon$ denotes
the presence of an outlier, whereas ${{\mathcal{I}}}^{i}_{k}=1$ indicates no
outlier in the $i$th dimension at time $k$. Since outliers can occur
independently at any instant, we assume that the elements of
$\bm{{\mathcal{I}}}_{k}$ are statistically independent of their past.
Additionally, we assume that the entries of $\bm{{\mathcal{I}}}_{k}$ to be
independent of each other since generally no knowledge of correlations between
outliers is available which are not easy to model anyway. Moreover, this
choice is motivated by the goal of inferential tractability. We also consider
$\bm{\mathcal{I}}_{k}$ and $\mathbf{x}_{k}$ to be statistically independent
since the outlier occurrence does not depend on the state value. The assumed
distribution of $\bm{\mathcal{I}}_{k}$ is given as

$p(\bm{\mathcal{I}}_{k})=\prod_{i=1}^{m}p({{\mathcal{I}}}^{i}_{k})=\prod_{i=1}^{m}(1-{\theta^{i}_{k}})\delta({{{\mathcal{I}}}^{i}_{k}}-\epsilon)+{\theta^{i}_{k}}\delta({{{\mathcal{I}}}^{i}_{k}}-1)$
(3)

where $\theta^{i}_{k}$ denotes the prior probability of no outlier in the
$i$th observation at time $k$. Further, the conditional measurement likelihood
given the current state $\mathbf{x}_{k}$ and the indicator
$\bm{\mathcal{I}}_{k}$, is proposed to be normally distributed as

Figure 2: Probabilistic graphical model for the proposed method

$\displaystyle
p(\mathbf{y}_{k}|\mathbf{x}_{k},\bm{\mathcal{I}}_{k})={\mathcal{N}}\big{(}\mathbf{y}_{k}|\mathbf{h}(\mathbf{x}_{k}),\mathbf{R}_{k}({\bm{\mathcal{{I}}}_{k}})\big{)}$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{m}|\mathbf{R}_{k}({\bm{\mathcal{{I}}}_{k}})|}}\mathrm{exp}\Big{\\{}{-}\mfrac{1}{2}{(\mathbf{y}_{k}-\mathbf{h}(\mathbf{x}_{k}))}^{\top}\mathbf{R}_{k}^{-1}(\bm{\mathcal{I}}_{k}$
$\displaystyle\ \ \ \ (\mathbf{y}_{k}-\mathbf{h}(\mathbf{x}_{k}))\Big{\\}}$
(4)

where $\mathbf{R}_{k}(\bm{\mathcal{I}}_{k})$

$\displaystyle\begin{bmatrix}{R}^{1,1}_{k}/{{{\mathcal{I}}}^{1}_{k}}&\dots&{R}^{1,m}_{k}\delta({{{\mathcal{I}}}^{1}_{k}}{-}1)\delta({{{\mathcal{I}}}^{m}_{k}}{-}1)\\\
\vdots&\ddots&\vdots\\\
{R}^{m,1}_{k}\delta({{{\mathcal{I}}}^{m}_{k}}{-}1)\delta({{{\mathcal{I}}}^{1}_{k}}{-}1)&\cdots&{R}^{m,m}_{k}/{{{\mathcal{I}}}^{m}_{k}}\end{bmatrix}$
(5)

$\mathbf{R}_{k}({\bm{\mathcal{{I}}}_{k}})$ is the modified covariance matrix
of the measurements considering the effect of outliers. The effect of
$\mathbf{R}_{k}({\bm{\mathcal{{I}}}_{k}})$ on the data generation process can
be understood by considering the possible values of ${\mathcal{I}}^{i}_{k}$.
In particular, ${\mathcal{I}}^{i}_{k}=\epsilon$ leads to a very large $i$th
diagonal entry of $\mathbf{R}_{k}({\bm{\mathcal{{I}}}_{k}})$, while placing
zeros at the remaining $i$th row and column of the matrix. Resultingly, when
an outlier occurs in the $i$th dimension its effect on state estimation is
minimized. Moreover, the $i$th dimension no longer has any correlation with
any other entry, ceasing to have any effect on any other dimension during
inference. This is in contrast to ${\mathcal{I}}^{i}_{k}=1$ which ensures the
diagonal element and the off-diagonal correlation entries with other non-
affected dimensions are preserved. Lastly, note that the conditional
likelihood is independent of the batch of all the historical observations
$\mathbf{y}_{1:{k-1}}$. Fig. 2 shows how the standard probabilistic graphical
model (PGM) is modified into the proposed PGM for devising outlier-robust
state estimators. The suggested PGM meets the modeling aims of describing the
nominal and corrupted data sufficiently while remaining docile for statistical
inference.

##  III Proposed Robust Filtering

In filtering, we are interested in the posterior distribution of
$\mathbf{x}_{k}$ conditioned on all the observations $\mathbf{y}_{1:{k}}$ that
have been observed till time $k$. For this objective, we can employ the Bayes
rule recursively. Given the proposed observation model, the analytical
expression of the joint posterior distribution of $\mathbf{x}_{k}$ and
$\bm{\mathcal{I}}_{k}$ conditioned on the set of all the observations
$\mathbf{y}_{1:{k}}$ is given as

$p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}})=\frac{p(\mathbf{y}_{k}|\bm{\mathcal{I}}_{k},\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{y}_{1:{k-1}})p(\bm{\mathcal{I}}_{k})}{p(\mathbf{y}_{k}|\mathbf{y}_{1:{k-1}})}$
(6)

Theoretically, the joint posterior can further be marginalized to obtain the
required posterior distribution $p(\mathbf{x}_{k}|\mathbf{y}_{1:{k}})$.
Assuming $p(\mathbf{x}_{k}|\mathbf{y}_{1:{k-1}})$ as a Gaussian distribution,
we need to run $2^{m}$ Kalman filters corresponding to each combination of
$\bm{\mathcal{I}}_{k}$ to obtain the posterior. This results in computational
complexity of around $\mathcal{O}(2^{m}m^{3})$ where $m^{3}$ appears due to
matrix inversions (ignoring sparsity). Moreover, we need to resort to a
debatable approximation of the resulting Gaussian mixture distribution for
$p(\mathbf{x}_{k}|\mathbf{y}_{1:{k}})$ as single Gaussian density for
recursive tractability. Therefore, this approach clearly becomes impractical
and unsuitable.

To get around the problem, we can possibly employ the standard VB method where
the product of VB marginals is conveniently used to approximate the joint
posterior. We assume the following factorization of the posterior

$p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}})\approx
q^{f}(\mathbf{x}_{k})\prod_{i}q^{f}({\mathcal{I}}^{i}_{k})$ (7)

The VB approximation aims to minimize the Kullback-Leibler (KL) divergence
between the product approximation and the true posterior and leads to the
following marginals [36]

$\displaystyle q^{f}(\mathbf{x}_{k})$
$\displaystyle\propto\exp\big{(}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}}))\big{\rangle}_{q^{f}({{\bm{\mathcal{I}}}_{k}})}\big{)}$
(8) $\displaystyle q^{f}({\mathcal{I}}^{i}_{k})$
$\displaystyle\propto\exp\big{(}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:k}))\big{\rangle}_{q^{f}(\mathbf{x}_{k})q^{f}(\bm{\mathcal{I}}^{i-}_{k})}\
\big{)}\ \forall\ i$ (9)

Using (8)-(9) alternately the VB marginals can be updated iteratively until
convergence. The procedure provides a useful way to approximate the true
marginals of the joint posterior by approximating these as
$p(\mathbf{x}_{k}|\mathbf{y}_{1:{k}})\approx{q^{f}(\mathbf{x}_{k})}$ and
$p(\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}})\approx\prod_{i}q^{f}({\mathcal{I}}^{i}_{k})$.

For our model, (8) becomes computationally unfriendly if we use the standard
VB approach. In fact, the same complexity order of $\mathcal{O}(m^{3}2^{m})$
appears as with the basic marginalization approach making this approach
intractable too. We elaborate more on it in the upcoming subsection.

### Expectation-Maximization as a particular case of variational Bayes

To deal with the complexity issue, instead of considering distributions we can
resort to point estimates for ${{\mathcal{I}}}^{i}_{k}$. In particular,
consider
$q^{f}({\mathcal{I}}^{i}_{k})=\delta({{\mathcal{I}}}^{i}_{k}-{\hat{{\mathcal{I}}}}^{i}_{k})$
where ${\hat{{\mathcal{I}}}}^{i}_{k}$ denotes the point approximation of
${{{\mathcal{I}}}}^{i}_{k}$. Consequently, the variational distributions can
be updated in an alternating manner in the Expectation (E) and Maximization
(M) steps in the EM algorithm given as [37]

#### E-Step

$\displaystyle{q^{f}(\bm{\mathbf{x}_{k}})}$
$\displaystyle=p(\mathbf{x}_{k}|\mathbf{y}_{1:{k}},\hat{\bm{\mathcal{I}}}_{k})\propto
p(\mathbf{x}_{k},\hat{\bm{\mathcal{I}}}_{k}|\mathbf{y}_{1:{k}})$ (10)

#### M-Step

${\hat{\mathcal{I}}}^{i}_{k}=\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},{\mathcal{I}}^{i}_{k},\hat{\bm{\mathcal{I}}}^{i-}_{k}|\mathbf{y}_{1:{k}})\big{\rangle}_{q^{f}(\mathbf{x}_{k})}$
(11)

where all ${\hat{\mathcal{I}}}^{i}_{k}$ in the M-Step are successively updated
using the latest estimates.

###  III-A Prediction

For filtering, we first obtain the predictive distribution
$p(\mathbf{x}_{k}|\mathbf{y}_{1:{k-1}})$ using the posterior distribution at
the previous instant $p(\mathbf{x}_{k-1}|\mathbf{y}_{1:{k-1}})$ approximated
as Gaussian
$q^{f}(\mathbf{x}_{k-1})\approx{\mathcal{N}}\left(\mathbf{x}_{k-1}|\mathbf{m}^{+}_{k-1},\mathbf{P}^{+}_{k-1}\right)$.
Using Gaussian (Kalman) filtering results we make the following approximation
[17]

$p(\mathbf{x}_{k}|\mathbf{y}_{1:{k-1}})\approx{\mathcal{N}}\left(\mathbf{x}_{k}|\mathbf{m}^{-}_{k},\mathbf{P}^{-}_{k}\right)$
(12)

where

$\displaystyle\mathbf{m}^{-}_{k}=$
$\displaystyle\int\mathbf{f}(\mathbf{x}_{k-1})\
{\mathcal{N}}\left(\mathbf{x}_{k-1}|\mathbf{m}^{+}_{k-1},\mathbf{P}^{+}_{k-1}\right)d\mathbf{x}_{k-1}$
(13) $\displaystyle\mathbf{P}^{-}_{k}=$
$\displaystyle\int\Big{\\{}(\mathbf{f}(\mathbf{x}_{k-1})-\mathbf{m}^{-}_{k})(\mathbf{f}(\mathbf{x}_{k-1})-\mathbf{m}^{-}_{k})^{\top}$
$\displaystyle{\mathcal{N}}({\mathbf{x}_{k-1}}|\mathbf{m}^{+}_{k-1},\mathbf{P}^{+}_{k-1})\Big{\\}}d\mathbf{x}_{k-1}+\mathbf{Q}_{k-1}$
(14)

###  III-B Update

Using the expressions of the prior distributions from (3) and (12) along with
the conditional measurement likelihood in (4), we can express the joint
posterior distribution as

$\displaystyle\
p\left(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:k}\right)\propto$
$\displaystyle\mfrac{\mathcal{N}(\mathbf{x}_{k}|\mathbf{m}_{k}^{-},\mathbf{P}_{k}^{-})}{\sqrt{(2\pi)^{m}\left|{\mathbf{R}}_{k}(\bm{\mathcal{I}}_{k})\right|}}\exp\Big{\\{}{-}\mfrac{1}{2}\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)^{\top}$
$\displaystyle{\mathbf{R}}_{k}^{-1}(\bm{\mathcal{I}}_{k})\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)\Big{\\}}\Big{\\{}\prod_{i}\left(1-\theta_{k}^{i}\right)$
$\displaystyle\delta\left(\mathcal{I}_{k}^{i}-\epsilon\right)+\theta_{k}^{i}\delta\left({\mathcal{I}}_{k}^{i}-1\right)\Big{\\}}$
(15)

####  III-B1 Derivation of $q^{f}(\mathbf{x}_{k})$

With the E-Step in (10) we can write

$\displaystyle\ q^{f}\left(\mathbf{x}_{k}\right)\propto$
$\displaystyle\exp\Big{\\{}{-}\mfrac{1}{2}\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)^{\top}{\mathbf{R}}_{k}^{-1}(\bm{\hat{\mathcal{I}}}_{k})\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)$
$\displaystyle{-}\mfrac{1}{2}\left(\mathbf{x}_{k}-\mathbf{m}_{k}^{-}\right)^{\top}({\mathbf{P}_{k}^{-}})^{-1}\left(\mathbf{x}_{k}-\mathbf{m}_{k}^{-}\right)\Big{\\}}$
(16)

where ${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$ assumes a
particular form resulting from the inversion of
${\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$ as described in Appendix A.

Note that we avoid evaluating
$\big{\langle}{\mathbf{R}}_{k}^{-1}(\bm{{\mathcal{I}}}_{k})\big{\rangle}_{q^{f}({{\bm{\mathcal{I}}}_{k}})}$
that would be required in the standard VB approach. This means that we are
able to evade the complexity level of around $\mathcal{O}(m^{3}2^{m})$ since
matrix inversion for each of the $2^{m}$ combinations are required to evaluate
the expectation. However, thanks to EM, we are now working with
${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$ which can be evaluated
with the maximum complexity of $\mathcal{O}(m^{3})$ (considering a fully
populated matrix).

To proceed further, we use the results of the general Gaussian filter, to
approximate $q^{f}(\mathbf{x}_{k})$ with a Gaussian distribution,
${\mathcal{N}}\left(\mathbf{x}_{k}|\mathbf{m}^{+}_{k},\mathbf{P}^{+}_{k}\right)$,
with parameters updated as

$\displaystyle\mathbf{m}^{+}_{k}$
$\displaystyle=\mathbf{m}^{-}_{k}+\mathbf{K}_{k}(\mathbf{y}_{k}-\bm{\mu}_{k})$
(17) $\displaystyle\mathbf{P}^{+}_{k}$
$\displaystyle=\mathbf{P}^{-}_{k}-\mathbf{C}_{k}\mathbf{K}^{\top}_{k}$ (18)

where

$\displaystyle\mathbf{K}_{k}=$ $\displaystyle\
\mathbf{C}_{k}(\mathbf{U}_{k}+{\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}}))^{-1}=\\{\mathbf{C}_{k}({\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$
$\displaystyle{-}{\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})(\mathbf{I}+\mathbf{U}_{k}{\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}}))^{-1}\mathbf{U}_{k}{\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}}))\\}$
$\displaystyle\bm{\mu}_{k}=$ $\displaystyle\int\mathbf{h}(\mathbf{x}_{k})\
{\mathcal{N}}\left(\mathbf{x}_{k}|\mathbf{m}^{-}_{k},\mathbf{P}^{-}_{k}\right)d\mathbf{x}_{k}$
$\displaystyle\mathbf{U}_{k}=$
$\displaystyle\int(\mathbf{h}(\mathbf{x}_{k})-\bm{\mu}_{k})(\mathbf{h}(\mathbf{x}_{k})-\bm{\mu}_{k})^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\mathbf{m}^{-}_{k},\mathbf{P}^{-}_{k})d\mathbf{x}_{k}$
$\displaystyle\mathbf{C}_{k}=$
$\displaystyle\int(\mathbf{x}_{k}-\mathbf{m}^{-}_{k})(\mathbf{h}(\mathbf{x}_{k})-\bm{\mu}_{k})^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\mathbf{m}^{-}_{k},\mathbf{P}^{-}_{k})d\mathbf{x}_{k}$

####  III-B2 Derivation of $\hat{\mathcal{I}}^{i}_{k}$

With the M-Step in (11) we can write

$\displaystyle{\hat{\mathcal{I}}}^{i}_{k}=$ $\displaystyle\
\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},{\mathcal{I}}^{i}_{k},\hat{\bm{\mathcal{I}}}^{i-}_{k}|\mathbf{y}_{1:{k}})\big{\rangle}_{q^{f}(\mathbf{x}_{k})}$
(19)

Using the Bayes rule we can proceed as

$\displaystyle{\hat{\mathcal{I}}}^{i}_{k}=$ $\displaystyle\
\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{y}_{{k}}|\mathbf{x}_{k},{\mathcal{I}}^{i}_{k},\hat{\bm{\mathcal{I}}}^{i-}_{k},\text{\sout{$\mathbf{y}_{1:{k-1}}$}}))\big{\rangle}_{q^{f}(\mathbf{x}_{k})}$
$\displaystyle+\mathrm{ln}(p({\mathcal{I}}^{i}_{k}|\text{\sout{${\mathbf{x}_{k}},\hat{\bm{\mathcal{I}}}^{i-}_{k},\mathbf{y}_{1:{k-1}}$}}))+const.$
(20)

where $const.$ is some constant and
$p(\mathbf{x}|\mathbf{y},\text{\sout{$\mathbf{z}$}})$ denotes the conditional
independence of $\mathbf{x}$ and $\mathbf{z}$ given $\mathbf{y}$. We can
further write

$\displaystyle{\hat{\mathcal{I}}}^{i}_{k}=$ $\displaystyle\
\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\Big{\\{}{-}\mfrac{1}{2}\mathrm{tr}\big{(}\mathbf{W}_{k}{\mathbf{R}}_{k}^{-1}({\mathcal{I}^{i}_{k}},\hat{\bm{\mathcal{I}}}^{i-}_{k})\big{)}{-}\mfrac{1}{2}\ln|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}},\hat{\bm{\mathcal{I}}}^{i-}_{k})|$
$\displaystyle+\ln\big{(}(1-\theta_{k}^{i})\delta(\mathcal{I}_{k}^{i}-\epsilon)+\theta_{k}^{i}\delta({\mathcal{I}}_{k}^{i}-1)\big{)}\Big{\\}}$
(21)

where
${\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}},\hat{\bm{\mathcal{I}}}^{i-}_{k})$
denotes ${\mathbf{R}}_{k}({\bm{\mathcal{I}}}_{k})$ evaluated at
${\bm{\mathcal{I}}}_{k}$ with it $i$th element as ${\mathcal{I}^{i}_{k}}$ and
remaining entries $\hat{\bm{\mathcal{I}}}^{i-}_{k}$ and

$\displaystyle\mathbf{W}_{k}=\int\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\mathbf{m}^{+}_{k},\mathbf{P}^{+}_{k})d\mathbf{x}_{k}$

Resultingly, $\hat{\mathcal{I}}^{i}_{k}$ can be determined as

$\displaystyle{\hat{\mathcal{I}}}^{i}_{k}$
$\displaystyle=\begin{cases}1&\text{if }\hat{\tau}^{i}_{k}\leq 0,\\\
0&\text{if }\hat{\tau}^{i}_{k}>0\end{cases}$ (22)

with

$\displaystyle\hat{\tau}^{i}_{k}=$ $\displaystyle\
\Big{\\{}\mathrm{tr}\big{(}\mathbf{W}_{k}\triangle{\hat{\mathbf{R}}}_{k}^{-1}\big{)}+\ln\Big{(}\frac{|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}=1},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}{|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}=\epsilon},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}\Big{)}$
$\displaystyle+2\ln\Big{(}\frac{1}{\theta^{i}_{k}}-1\Big{)}\Big{\\}}$ (23)

where

$\displaystyle\triangle{\hat{\mathbf{R}}}_{k}^{-1}=({\mathbf{R}}_{k}^{-1}({\mathcal{I}^{i}_{k}=1},\hat{\bm{\mathcal{I}}}^{i-}_{k})-{\mathbf{R}}_{k}^{-1}({\mathcal{I}^{i}_{k}=\epsilon},\hat{\bm{\mathcal{I}}}^{i-}_{k}))$
(24)

Using the steps outlined in Appendix B, we can further simplify
$\hat{\tau}^{i}_{k}$ as

$\displaystyle\hat{\tau}^{i}_{k}=$ $\displaystyle\
\Big{\\{}\mathrm{tr}\big{(}\mathbf{W}_{k}\triangle{\hat{\mathbf{R}}}_{k}^{-1}\big{)}+\ln\Bigg{|}\textbf{I}-\frac{\mathbf{R}^{-i,i}_{k}\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}}{{R}^{i,i}_{k}}\Bigg{|}$
$\displaystyle+\ln(\epsilon)+2\ln\Big{(}\frac{1}{\theta^{i}_{k}}-1\Big{)}\Big{\\}}$
(25)

where $\hat{\mathbf{R}}^{-i,-i}_{k}$ is the submatrix of
${\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$ corresponding to entries of
$\hat{\bm{\mathcal{I}}}^{i-}_{k}$. $\mathbf{R}^{i,-i}_{k}$ and
$\mathbf{R}^{-i,i}_{k}$ contain the measurement covariances between ith and
rest of the dimensions.

Initialize $\mathbf{m}^{+}_{0},\mathbf{P}^{+}_{0}$;

for  _$k=1,2...K$_ do

Initialize
$\theta^{i}_{k},{\bm{\hat{\mathcal{I}}}_{k}},\mathbf{Q}_{k},\mathbf{R}_{k}$;

Prediction:

Evaluate $\mathbf{m}^{-}_{k},\mathbf{P}^{-}_{k}$ with (13) and (14);

Update:

while  _not converged_ do

Update ${\mathbf{m}^{+}_{k}}$ and ${\mathbf{P}^{+}_{k}}$ with (17)-(18);

Update ${\hat{\mathcal{I}}^{i}_{k}}$ $\forall\ i$ with (22);

end while

end for

Algorithm 1 The proposed filter: EMORF

Though we can directly evaluate $\triangle{\hat{\bm{\mathbf{R}}}}_{k}^{-1}$ in
(24), we can save computations by avoiding repetitive calculations. To this
end, we first need to compute

$\displaystyle\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}=\begin{bmatrix}\bm{\Xi}^{i,i}&\bm{\Xi}^{i,-i}\\\
\bm{\Xi}^{-i,i}&\bm{\Xi}^{-i,-i}\end{bmatrix}$ (26)

with

$\displaystyle\bm{\Xi}^{i,i}$
$\displaystyle=\frac{1}{R^{i,i}_{k}-\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}}-\frac{\epsilon}{R^{i,i}_{k}}$
(27) $\displaystyle\bm{\Xi}^{i,-i}$
$\displaystyle=-\frac{\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}}{R^{i,i}_{k}-\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}}$
(28) $\displaystyle\bm{\Xi}^{-i,i}$
$\displaystyle=-\frac{(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}}{R^{i,i}_{k}-\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}}$
(29) $\displaystyle\bm{\Xi}^{-i,-i}$
$\displaystyle=\frac{(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}}{R^{i,i}_{k}-\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}}$
(30)

where the $i$th row/column entries in
$\triangle{\hat{\bm{\mathbf{R}}}}_{k}^{-1}$ have been conveniently swapped
with the first row/column elements to obtain
$\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}$. By swapping the first
row/column entries of $\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}$ to the
$i$th row/column positions we can reclaim
$\triangle{\hat{\bm{\mathbf{R}}}}_{k}^{-1}$. Appendix C provides further
details in this regard.

The resulting EM-based outlier-robust filter (EMORF) is outlined as Algorithm
1. For the convergence criterion, we suggest using the ratio of the L2 norm of
the difference of the state estimates from the current and previous VB
iterations and the L2 norm of the estimate from the previous iteration. This
criterion has been commonly chosen in similar robust filters [18, 33].

###  III-C VB factorization of
$p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}})$ and the associated
computational overhead

Note that for better accuracy we can factorize
$p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:{k}})\approx
q^{f}(\mathbf{x}_{k})q^{f}(\bm{\mathcal{I}}_{k})$ instead of forcing
independence between all ${\mathcal{I}}^{i}_{k}$. In this case, expression for
evaluating $q^{f}(\mathbf{x}_{k})$ remains same as in (8). However, this
choice leads to the following VB marginal distribution of
$\bm{\mathcal{I}}_{k}$

$\displaystyle q^{f}(\bm{\mathcal{I}}_{k})$
$\displaystyle\propto\exp\big{(}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},\bm{\mathcal{I}}_{k}|\mathbf{y}_{1:k}))\big{\rangle}_{q^{f}(\mathbf{x}_{k})}\
\big{)}$ (31)

This results in a modified M-Step in the EM algorithm given as

#### M-Step

${\hat{\bm{\mathcal{I}}}}_{k}=\underset{\bm{{\mathcal{I}}}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{k},{\bm{\mathcal{I}}}_{k}|\mathbf{y}_{1:{k}}))\big{\rangle}_{q^{f}(\mathbf{x}_{k})}$
(32)

Proceeding further, with the prediction and update steps during inference, we
can arrive at

$\displaystyle{\hat{\bm{\mathcal{I}}}}_{k}=$ $\displaystyle\
\underset{{\bm{\mathcal{I}}}_{k}}{\mathrm{argmax}}\Big{\\{}{-}\mfrac{1}{2}\mathrm{tr}\big{(}\mathbf{W}_{k}{\mathbf{R}}_{k}^{-1}({\bm{\mathcal{I}}}_{k})\big{)}{-}\mfrac{1}{2}\ln|{\mathbf{R}}_{k}({\bm{\mathcal{I}}}_{k})|$
$\displaystyle+\ln\big{(}(1-\theta_{k}^{i})\delta(\mathcal{I}_{k}^{i}-\epsilon)+\theta_{k}^{i}\delta({\mathcal{I}}_{k}^{i}-1)\big{)}\Big{\\}}$
(33)

It is not hard to notice that determining ${\hat{\bm{\mathcal{I}}}}_{k}$ using
(33) involves tedious calculations. In fact, we run into the same
computational difficulty that we have been dodging. To arrive at the result,
we need to evaluate the inverses and determinants, for each of the $2^{m}$
combinations corresponding to the entries of ${\bm{\mathcal{I}}}_{k}$. This
entails the bothering complexity level of $\mathcal{O}(m^{3}2^{m})$.

With the proposed factorization in (7), we obtain a more practical and
scalable algorithm. The resulting complexity is $\mathcal{O}(m^{4})$ following
from the evaluation of matrix inverses and determinants for calculating each
of the $\hat{\mathcal{I}}^{i}_{k}$ $\forall\ i=1\cdots m$ in (22).

###  III-D Connection between EMORF and SORF [18]

Since the construction of EMORF is motivated by selective observations
rejecting filter (SORF) [18], it is insightful to discuss their connection. We
derived SORF considering a diagonal measurement covariance matrix
$\mathbf{R}_{k}$ [18]. In SORF, we used distributional estimates for
$\bm{\mathcal{I}}_{k}$ since it did not induce any significant computational
strain. It is instructive to remark here that the if we use point estimates of
$\bm{\mathcal{I}}_{k}$ in SORF it becomes a special case of EMORF. In
particular, the point estimates for ${\mathcal{I}}^{i}_{k}\ \forall\ i$ in
[18] can be obtained with the following criterion

$\displaystyle{\hat{\mathcal{I}}}^{i}_{k}$
$\displaystyle=\begin{cases}1&\text{if }\bar{\tau}^{i}_{k}\leq 0,\\\
0&\text{if }\bar{\tau}^{i}_{k}>0\end{cases}$ (34)

To deduce ${\hat{\mathcal{I}}}^{i}_{k}=1$, the following should hold

$\displaystyle\Omega^{i}_{k}$ $\displaystyle\geq 1-\Omega^{i}_{k}$ (35)
$\displaystyle\Omega^{i}_{k}$ $\displaystyle\geq 0.5$ (36)

where $\Omega^{i}_{k}$ denotes the posterior probability of
${{\mathcal{I}}}^{i}_{k}=1$. Using the expression of $\Omega^{i}_{k}$ from
[18] we can write (36) as

$\displaystyle\frac{1}{1+{\sqrt{\epsilon}}\left({\frac{1}{\theta^{i}_{k}}-1}\right){\exp\left({\frac{W^{i,i}_{k}}{2R^{i,i}_{k}}(1-\epsilon)}\right)}}\geq
0.5$ (37)

leading to

$\displaystyle\bar{\tau}^{i}_{k}\leq 0$ (38)

where

$\displaystyle\bar{\tau}^{i}_{k}=$
$\displaystyle{{\frac{W^{i,i}_{k}}{R^{i,i}_{k}}(1-\epsilon)}}+\ln({\epsilon})+2\ln\Big{(}{\frac{1}{\theta^{i}_{k}}-1}\Big{)}$
(39)

which can be recognized as the particular case of (25) given
$\mathbf{R}^{i,-i}_{k}$ and $\mathbf{R}^{-i,i}_{k}$ vanish as $\mathbf{R}_{k}$
is considered to be diagonal. The first term in $\eqref{eqn_fl_up10_i}$
reduces to ${{\frac{W^{i,i}_{k}}{R^{i,i}_{k}}(1-\epsilon)}}$ given
$\bm{\Xi}^{i,i}=\frac{1-\epsilon}{R^{i,i}_{k}}$ and
$\bm{\Xi}^{i,-i},\bm{\Xi}^{-i,i}$ and $\bm{\Xi}^{-i,-i}$ all reduce to zero.
Moreover, the second term in (25) also disappears resulting in (39).

###  III-E Choice of the parameters $\theta^{i}_{k}$ and $\epsilon$

For EMORF we propose setting the parameters $\theta^{i}_{k}$ and $\epsilon$
the same as in SORF. Specifically, we suggest choosing a neutral value of 0.5
or an uninformative prior for $\theta^{i}_{k}$. The Bayes-Laplace and the
maximum entropy approaches for obtaining uninformative prior for a finite-
valued parameter lead to the choice of the uniform prior distribution [38,
39]. Moreover, the selection has been justified in the design of outlier-
resistant filters assuming no prior information about the outliers statistics
is available [18, 32]. For $\epsilon$ we recommend its value to be close to
zero since the exact value of $0$ denies the VB/EM updates as in [18].

##  IV Proposed Robust Smoothing

In smoothing, our interest lies in determining the posterior distribution of
all the states $\mathbf{x}_{1:K}$ conditioned on the batch of all the
observations $\mathbf{y}_{1:{K}}$. With that goal, we take a similar approach
to filtering and approximate the joint posterior distribution as a product of
marginals

$p(\mathbf{x}_{1:K},\bm{\mathcal{I}}_{1:K}|\mathbf{y}_{1:{K}})\approx
q^{s}(\mathbf{x}_{1:K})\prod_{k}\prod_{i}q^{s}({\mathcal{I}}^{i}_{k})$ (40)

where the true marginals are approximated as
$p(\mathbf{x}_{1:K}|\mathbf{y}_{1:{K}})\approx{q^{s}(\mathbf{x}_{1:K})}$ and
$p(\bm{\mathcal{I}}_{1:K}|\mathbf{y}_{1:{K}})\approx{\prod_{k}\prod_{i}q^{s}({\mathcal{I}}^{i}_{k})}$.
Let us assume
$q^{s}({\mathcal{I}}^{i}_{k})=\delta({{\mathcal{I}}}^{i}_{k}-{\breve{{\mathcal{I}}}}^{i}_{k})$
where ${\breve{{\mathcal{I}}}}^{i}_{k}$ denotes the point approximation of
${{{\mathcal{I}}}}^{i}_{k}$. Consequently, the EM steps are given as

#### E-Step

$\displaystyle{q^{s}(\bm{\mathbf{x}}{1:K})}$
$\displaystyle=p(\mathbf{x}_{1:K}|\mathbf{y}_{1:{K}},\breve{\bm{\mathcal{I}}}_{1:K})\propto
p(\mathbf{x}_{1:K},\breve{\bm{\mathcal{I}}}_{1:K}|\mathbf{y}_{1:{K}})$ (41)

#### M-Step

${\breve{\mathcal{I}}}^{i}_{k}=\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{1:K},{\mathcal{I}}^{i}_{k},\breve{\bm{\mathcal{I}}}^{i-}_{k},\breve{\bm{\mathcal{I}}}_{k{-}}|\mathbf{y}_{1:{K}})\big{\rangle}_{q^{s}(\mathbf{x}_{1:K})}$
(42)

where all ${\hat{\mathcal{I}}}^{i}_{k}$ in the M-Step are sequentially updated
using the latest estimates.

####  IV-1 Derivation of $q^{s}\left(\mathbf{x}_{1:K}\right)$

With the E-Step in (41) we can write

$\displaystyle{q^{s}(\bm{\mathbf{x}}{1:K})}$ $\displaystyle\propto
p(\mathbf{y}_{{1:K}}|\mathbf{x}_{{1:K}},\breve{\bm{\mathcal{I}}}_{1:K})p(\mathbf{x}_{{1:K}})$
(43) $\displaystyle{q^{s}(\bm{\mathbf{x}}{1:K})}$
$\displaystyle\propto\prod_{k}p(\mathbf{y}_{{k}}|\mathbf{x}_{{k}},\breve{\bm{\mathcal{I}}}_{k})p(\mathbf{x}_{{k}}|\mathbf{x}_{{k-1}})$
(44)

We can identify that ${q^{s}(\bm{\mathbf{x}}{1:K})}$ can be approximated as a
Gaussian distribution from the results of general Gaussian RTS smoothing [17].
Using the forward and backward passes, we can determine the parameters of
${q^{s}(\bm{\mathbf{x}}{k})}\sim\mathcal{N}(\bm{\mathfrak{m}}^{s}_{k},\bm{\mathcal{P}}^{s}_{k})$,
which denotes the marginalized densities of ${q^{s}(\bm{\mathbf{x}}{1:K})}$.

#### Forward pass

The forward pass essentially involves the filtering equations given as

$\displaystyle\bm{\mathfrak{m}}^{-}_{k}=$
$\displaystyle\int\mathbf{f}(\mathbf{x}_{k-1})\
{\mathcal{N}}\left(\mathbf{x}_{k-1}|\bm{\mathfrak{m}}^{+}_{k-1},\bm{\mathcal{P}}^{+}_{k-1}\right)d\mathbf{x}_{k-1}$
(45) $\displaystyle\bm{\mathcal{P}}^{-}_{k}=$
$\displaystyle\int\Big{\\{}(\mathbf{f}(\mathbf{x}_{k-1})-\bm{\mathfrak{m}}^{-}_{k})(\mathbf{f}(\mathbf{x}_{k-1})-\bm{\mathfrak{m}}^{-}_{k})^{\top}$
(46)
$\displaystyle{\mathcal{N}}({\mathbf{x}_{k-1}}|\bm{\mathfrak{m}}^{+}_{k-1},\bm{\mathcal{P}}^{+}_{k-1})\Big{\\}}d\mathbf{x}_{k-1}+\mathbf{Q}_{k-1}$
(47) $\displaystyle\bm{\mathfrak{m}}^{+}_{k}=$ $\displaystyle\
\bm{\mathfrak{m}}^{-}_{k}+\bm{\mathcal{K}}_{k}(\mathbf{y}_{k}-\bm{\nu}_{k})$
(48) $\displaystyle\bm{\mathcal{P}}^{+}_{k}=$ $\displaystyle\
\bm{\mathcal{P}}^{-}_{k}-\bm{\mathcal{C}}_{k}\bm{\mathcal{K}}^{\top}_{k}$ (49)

where

$\displaystyle\bm{\mathcal{K}}_{k}=$ $\displaystyle\
\bm{\mathcal{C}}_{k}(\bm{\mathcal{U}}_{k}+{\mathbf{R}}_{k}({\bm{\breve{\mathcal{I}}}_{k}}))^{-1}=\\{\bm{\mathcal{C}}_{k}({\mathbf{R}}_{k}^{-1}({\bm{\breve{\mathcal{I}}}_{k}})$
$\displaystyle{-}{\mathbf{R}}_{k}^{-1}({\bm{\breve{\mathcal{I}}}_{k}})(\mathbf{I}+\bm{\mathcal{U}}_{k}{\mathbf{R}}_{k}^{-1}({\bm{\breve{\mathcal{I}}}_{k}}))^{-1}\bm{\mathcal{U}}_{k}{\mathbf{R}}_{k}^{-1}({\bm{\breve{\mathcal{I}}}_{k}}))\\}$
$\displaystyle\bm{\nu}_{k}=$ $\displaystyle\int\mathbf{h}(\mathbf{x}_{k})\
{\mathcal{N}}\left(\mathbf{x}_{k}|\bm{\mathfrak{m}}^{-}_{k},\bm{\mathcal{P}}^{-}_{k}\right)d\mathbf{x}_{k}$
$\displaystyle\bm{\mathcal{U}}_{k}=$
$\displaystyle\int(\mathbf{h}(\mathbf{x}_{k})-\bm{\nu}_{k})(\mathbf{h}(\mathbf{x}_{k})-\bm{\nu}_{k})^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\bm{\mathfrak{m}}^{-}_{k},\bm{\mathcal{P}}^{-}_{k})d\mathbf{x}_{k}$
$\displaystyle\bm{\mathcal{C}}_{k}=$
$\displaystyle\int(\mathbf{x}_{k}-\bm{\mathfrak{m}}^{-}_{k})(\mathbf{h}(\mathbf{x}_{k})-\bm{\nu}_{k})^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\bm{\mathfrak{m}}^{-}_{k},\bm{\mathcal{P}}^{-}_{k})d\mathbf{x}_{k}$

Note that ${\mathbf{R}}_{k}({\bm{\breve{\mathcal{I}}}_{k}})$ and
${\mathbf{R}}_{k}^{-1}({\bm{\breve{\mathcal{I}}}_{k}})$ can be evaluated
similar to ${\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$ and
${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$.

#### Backward pass

The backward pass can be completed as

$\displaystyle\bm{\mathcal{L}}_{k+1}=$
$\displaystyle\int(\mathbf{x}_{k}-\bm{\mathfrak{m}}^{+}_{k})(\mathbf{f}\left(\mathbf{x}_{k}\right)-\bm{\mathfrak{m}}_{k+1}^{-})^{\top}{\mathcal{N}}\left(\mathbf{x}_{k}|\bm{\mathfrak{m}}^{+}_{k},\bm{\mathcal{P}}^{+}_{k}\right)d\mathbf{x}_{k}$
(50) $\displaystyle\bm{\mathcal{G}}_{k}=$ $\displaystyle\
\bm{\mathcal{L}}_{k+1}\left(\bm{\mathcal{P}}_{k+1}^{-}\right)^{-1}$ (51)
$\displaystyle\bm{\mathfrak{m}}_{k}^{{s}}=$ $\displaystyle\
\bm{\mathfrak{m}}^{+}_{k}+\bm{\mathcal{G}}_{k}\left(\bm{\mathfrak{m}}_{k+1}^{{s}}-\bm{\mathfrak{m}}_{k+1}^{-}\right)$
(52) $\displaystyle\bm{\mathcal{P}}_{k}^{{s}}=$ $\displaystyle\
\bm{\mathcal{P}}^{+}_{k}+\bm{\mathcal{G}}_{k}\left(\bm{\mathcal{P}}_{k+1}^{{s}}-\bm{\mathcal{P}}_{k+1}^{-}\right)\bm{\mathcal{G}}_{k}^{\top}$
(53)

####  IV-2 Derivation of $\breve{\mathcal{I}}^{i}_{k}$

With the M-Step in (42) we can write

$\displaystyle{\breve{\mathcal{I}}}^{i}_{k}=\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{x}_{1:K},{\mathcal{I}}^{i}_{k},\breve{\bm{\mathcal{I}}}^{i-}_{k},\breve{\bm{\mathcal{I}}}_{k{-}}|\mathbf{y}_{1:{k}})\big{\rangle}_{q^{s}(\mathbf{x}_{1:K})}$
(54)

Using the Bayes rule we can proceed as

$\displaystyle{\breve{\mathcal{I}}}^{i}_{k}=$ $\displaystyle\
\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{y}_{{k}}|\mathbf{x}_{k},{\mathcal{I}}^{i}_{k},\breve{\bm{\mathcal{I}}}^{i-}_{k},\text{\sout{$\breve{\bm{\mathcal{I}}}_{k{-}},{\mathbf{y}_{k-}}$}}))$
$\displaystyle+\mathrm{ln}(p(\mathbf{x}_{1:K},{\mathcal{I}}^{i}_{k},\breve{\bm{\mathcal{I}}}^{i-}_{k},\breve{\bm{\mathcal{I}}}_{k{-}}|\mathbf{y}_{k-}))\big{\rangle}_{q^{s}(\mathbf{x}_{1:K})}+const.$
(55)

which leads to

$\displaystyle{\breve{\mathcal{I}}}^{i}_{k}=$ $\displaystyle\
\underset{{{\mathcal{I}}}^{i}_{k}}{\mathrm{argmax}}\big{\langle}\mathrm{ln}(p(\mathbf{y}_{{k}}|\mathbf{x}_{k},{\mathcal{I}}^{i}_{k},\breve{\bm{\mathcal{I}}}^{i-}_{k}))\big{\rangle}_{q^{s}(\mathbf{x}_{k})}$
$\displaystyle+\mathrm{ln}(p({\mathcal{I}}^{i}_{k}|\text{\sout{${\mathbf{x}_{1:K}},{\breve{\bm{\mathcal{I}}}^{i-}_{k}},{\hat{\bm{\mathcal{I}}}_{k{-}}},{\mathbf{y}_{k-}}$}}))+const.$
(56)

which is similar to (20) except that the expectation is taken with respect to
the marginal smoothing distribution ${q^{s}(\bm{\mathbf{x}}{k})}$.
Consequently, $\breve{\mathcal{I}}^{i}_{k}$ can be determined as

$\displaystyle{\breve{\mathcal{I}}}^{i}_{k}$
$\displaystyle=\begin{cases}1&\text{if }\breve{\tau}^{i}_{k}\leq 0,\\\
0&\text{if }\breve{\tau}^{i}_{k}>0\end{cases}$ (57)

where

$\displaystyle\breve{\tau}^{i}_{k}=$ $\displaystyle\
\Big{\\{}\mathrm{tr}\big{(}\bm{\mathcal{W}}_{k}\triangle{\breve{\mathbf{R}}}_{k}^{-1}\big{)}+\ln\Bigg{|}\textbf{I}-\frac{\mathbf{R}^{-i,i}_{k}\mathbf{R}^{i,-i}_{k}(\breve{\mathbf{R}}^{-i,-i}_{k})^{-1}}{{R}^{i,i}_{k}}\Bigg{|}$
$\displaystyle+\ln(\epsilon)+2\ln\Big{(}\frac{1}{\theta^{i}_{k}}-1\Big{)}\Big{\\}}$
(58)

with

$\displaystyle\triangle{\breve{\mathbf{R}}}_{k}^{-1}=({\mathbf{R}}_{k}^{-1}({\mathcal{I}^{i}_{k}=1},\breve{\bm{\mathcal{I}}}^{i-}_{k})-{\mathbf{R}}_{k}^{-1}({\mathcal{I}^{i}_{k}=\epsilon},\breve{\bm{\mathcal{I}}}^{i-}_{k}))$
(59)

Initialize $\bm{\mathfrak{m}}^{+}_{0},\bm{\mathcal{P}}^{+}_{0}$;

while  _not converged_ do

for  _$k=1,2...K$_ do

Initialize
$\theta^{i}_{k},{\bm{\breve{\mathcal{I}}}_{k}},\mathbf{Q}_{k},\mathbf{R}_{k}$;

Evaluate $\bm{\mathfrak{m}}^{-}_{k},\bm{\mathcal{P}}^{-}_{k}$ with (45)-(46);

Evaluate $\bm{\mathfrak{m}}^{+}_{k},\bm{\mathcal{P}}^{+}_{k}$ with (48)-(49);

end for

$\bm{\mathfrak{m}}^{s}_{K}=\bm{\mathfrak{m}}^{+}_{K}$;

$\bm{\mathcal{P}}^{s}_{K}=\bm{\mathcal{P}}^{+}_{K}$ ;

for  _$k=K-1,...1$_ do

Evaluate $\bm{\mathfrak{m}}^{s}_{k},\bm{\mathcal{P}}^{s}_{k}$ with (52)-(53);

end for

for  _$k=1,2...K$_ do

Update ${\breve{\mathcal{I}}^{i}_{k}}$ $\forall\ i$ with (57);

end for

end while

Algorithm 2 The proposed smoother: EMORS

which can be be calculated similar to $\triangle{\hat{\mathbf{R}}}_{k}^{-1}$
in (24). $\breve{\mathbf{R}}^{-i,-i}_{k}$ denotes the submatrix of
${\mathbf{R}}_{k}({\bm{\breve{\mathcal{I}}}_{k}})$ corresponding to entries of
$\breve{\bm{\mathcal{I}}}^{i-}_{k}$ and

$\displaystyle\bm{\mathcal{W}}_{k}=\int\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)\left(\mathbf{y}_{k}-\mathbf{h}\left(\mathbf{x}_{k}\right)\right)^{\top}{\mathcal{N}}({\mathbf{x}_{k}}|\bm{\mathfrak{m}}^{s}_{k},\bm{\mathcal{P}}^{s}_{k})d\mathbf{x}_{k}$

The resulting EM-based outlier-robust smoother (EMORS) is outlined as
Algorithm 2. We suggest using the same convergence criterion and parameters as
for robust filtering.

##  V Performance Bounds

It is useful to determine the performance bounds of outlier-discarding state
estimators considering correlated measurement noise. We evaluate the
estimation bounds of filtering and smoothing approaches that are perfect
outlier rejectors, having complete knowledge of the instances of outlier
occurrences. In particular, we assume that the measurement covariance matrix
is a function of perfectly known values of $\bm{\mathcal{I}}_{k}$ given as
$\mathbf{R}_{k}(\bm{\mathcal{I}}_{k})$. In this case,
${\mathcal{I}}^{i}_{k}=0$ means rejection of the $i$th corrupted dimension,
whereas ${\mathcal{I}}^{i}_{k}=1$ denotes inclusion of the $i$th measurement.
Resultingly, $\mathbf{R}^{-1}_{k}(\bm{\mathcal{I}}_{k})$ has zeros at the
diagonals, rows, and columns corresponding to dimensions for which
${\mathcal{I}}^{i}_{k}=0$. Remaining submatrix of
$\mathbf{R}^{-1}_{k}(\bm{\mathcal{I}}_{k})$ can be evaluated as the inverse of
submatrix of $\mathbf{R}_{k}$ considering the dimensions with
${\mathcal{I}}^{i}_{k}=1$.

Note that we set ${\mathcal{I}}^{i}_{k}=\epsilon$, not exactly as $0$, for
outlier rejection in the proposed state estimators as it declines inference.
However, for evaluating performance bounds this choice is appropriate
resulting in perfect outlier rejection. Also note that during robust state
estimation, we do not exactly know $\bm{\mathcal{I}}_{k}$ apriori and model it
statistically for subsequent inference. The use of perfectly known
$\bm{\mathcal{I}}_{k}$ for estimation bounds gives us an idea of how well we
can estimate the state if outliers are somehow perfectly detected and
rejected.

We evaluate BCRBs for the perfect rejector for the model in (1)-(2) that have
been corrupted with measurement outliers for both filtering and smoothing.

###  V-A Filtering

For the estimation error of $\mathbf{x}_{k}$ during filtering, the BCRB matrix
can be written as [40]

$\text{BCRB}^{f}_{k}\triangleq({\mathbf{J}^{+}_{k}})^{-1}$ (49)

where the corresponding filtering Fisher information matrix (FIM) denoted as
$\mathbf{J}^{+}_{k}$ can be evaluated recursively as

$\displaystyle\mathbf{J}^{-}_{k}=\mathbf{D}_{k-1}^{22}(1)-\mathbf{D}_{k-1}^{21}\left(\mathbf{J}^{+}_{k-1}+\mathbf{D}_{k-1}^{11}\right)^{-1}\mathbf{D}_{k-1}^{12}$
(60)
$\displaystyle\mathbf{J}^{+}_{k}=\mathbf{J}^{-}_{k}+\mathbf{D}_{k-1}^{22}(2)$
(61)

where $\mathbf{J}^{+}_{0}=\langle-\Delta_{\mathbf{x}_{0}}^{\mathbf{x}_{0}}\ln
p(\mathbf{x}_{0})\rangle_{p(\mathbf{x}_{0})}$ and

$\displaystyle\Delta_{\Psi}^{\Theta}=\nabla_{\Psi}\nabla_{\Theta}^{\top}$ (62)
$\displaystyle\nabla_{\Theta}=\left[\frac{\partial}{\partial\Theta_{1}},\ldots,\frac{\partial}{\partial\Theta_{r}}\right]^{\top}$
(63)

$\displaystyle\mathbf{D}_{k}^{11}$
$\displaystyle=\langle-\Delta_{\mathbf{x}_{k}}^{\mathbf{x}_{k}}\ln
p\left(\mathbf{x}_{k+1}\mid\mathbf{x}_{k}\right)\rangle_{p(\mathbf{x}_{k+1},\mathbf{x}_{k})}$
(64) $\displaystyle\mathbf{D}_{k}^{12}$
$\displaystyle=\langle-\Delta_{\mathbf{x}_{k}}^{\mathbf{x}_{k+1}}\ln
p\left(\mathbf{x}_{k+1}\mid\mathbf{x}_{k}\right)\rangle_{p(\mathbf{x}_{k+1},\mathbf{x}_{k})}$
(65) $\displaystyle\mathbf{D}_{k}^{21}$
$\displaystyle=\langle-\Delta_{\mathbf{x}_{k+1}}^{\mathbf{x}_{k}}\ln
p\left(\mathbf{x}_{k+1}\mid\mathbf{x}_{k}\right)\rangle_{p(\mathbf{x}_{k+1},\mathbf{x}_{k})}=\left(\mathbf{D}_{k}^{12}\right)^{\top}$
(66) $\displaystyle\mathbf{D}_{k}^{22}$
$\displaystyle=\mathbf{D}_{k}^{22}(1)+\mathbf{D}_{k}^{22}(2)$ (67)
$\displaystyle\mathbf{D}_{k}^{22}(1)$
$\displaystyle=\langle-\Delta_{\mathbf{x}_{k+1}}^{\mathbf{x}_{k+1}}\ln
p\left(\mathbf{x}_{k+1}\mid\mathbf{x}_{k}\right)\rangle_{p(\mathbf{x}_{k+1},\mathbf{x}_{k})}$
(68) $\displaystyle\mathbf{D}_{k}^{22}(2)$
$\displaystyle=\langle-\Delta_{\mathbf{x}_{k+1}}^{\mathbf{x}_{k+1}}\ln
p\left(\mathbf{y}_{k+1}\mid\mathbf{x}_{k+1}\right)\rangle_{p(\mathbf{y}_{k+1},\mathbf{x}_{k+1})}$
(69)

The bound is valid given the existence of the following derivatives and
expectations terms for an asymptotically unbiased estimator [41]. For the
perfect rejector considering the system model in (1)-(2) that is infested with
observation outliers we can write

$\displaystyle\mathbf{D}_{k}^{11}$
$\displaystyle=\langle\tilde{\mathbf{F}}^{\top}(\mathbf{x}_{k})\mathbf{Q}_{k}^{-1}\tilde{\mathbf{F}}(\mathbf{x}_{k})\rangle_{p(\mathbf{x}_{k})}$
(70) $\displaystyle\mathbf{D}_{k}^{12}$
$\displaystyle=-\langle\tilde{\mathbf{F}}^{\top}(\mathbf{x}_{k})\rangle_{p(\mathbf{x}_{k})}\mathbf{Q}_{k}^{-1}$
(71) $\displaystyle\mathbf{D}_{k}^{22}(1)$ $\displaystyle=\mathbf{Q}_{k}^{-1}$
(72) $\displaystyle\mathbf{D}_{k}^{22}(2)$
$\displaystyle=\langle\tilde{\mathbf{H}}^{\top}(\mathbf{x}_{k+1})\mathbf{R}_{k+1}^{-1}(\bm{\mathcal{I}}_{k+1})\tilde{\mathbf{H}}(\mathbf{x}_{k+1})\rangle_{p(\mathbf{x}_{k+1})}$
(73)

where $\tilde{\mathbf{F}}(.)$ and $\tilde{\mathbf{H}}(.)$ are the Jacobians of
the transformations $\mathbf{f}(.)$ and $\mathbf{h}(.)$ respectively.

###  V-B Smoothing

Similarly, for the estimation error of $\mathbf{x}_{k}$ during smoothing, the
BCRB matrix can be written as [42]

$\text{BCRB}^{s}_{k}\triangleq({\mathbf{J}^{s}_{k}})^{-1}$ (74)

where $\mathbf{J}^{s}_{K}=\mathbf{J}^{+}_{K}$. We can compute the associated
smoothing FIM denoted as $\mathbf{J}^{s}_{k}$ recursively as

$\displaystyle\mathbf{J}^{s}_{k}=\mathbf{J}^{+}_{k}+\mathbf{D}_{k}^{11}-\mathbf{D}_{k}^{12}\left(\mathbf{D}_{k}^{22}(1)+\mathbf{J}^{s}_{k+1}+\mathbf{J}^{-}_{k+1}\right)^{-1}\mathbf{D}_{k}^{21}$
(75)

##  VI Numerical Experiments

To test the performance of the proposed outlier-resilient state estimators, we
carry out several numerical experiments. We use Matlab on a computer powered
by an Intel i7-8550U processor. All the experiments were conducted while
considering SI units.

Figure 3: Target tracking test example setup

For performance evaluation, we resort to a target tracking problem with TDOA-
based range measurements inspired by [21]. Fig. 3 shows the setup of the
considered example. Owing to the use of a common reference sensor to obtain
the TDOA observations, from the difference of the time of arrival (TOA)
measurements, the resulting covariance matrix becomes fully populated.

We consider the process equation for the target assuming an unknown turning
rate as [33]

$\displaystyle\mathbf{x}_{k}=\mathbf{f}(\mathbf{x}_{k-1})+\mathbf{q}_{k-1}$
(76)

with

$\displaystyle\mathbf{f}(\mathbf{x}_{k-1})$
$\displaystyle=\begin{bmatrix}\text{1}&\frac{\text{sin}(\omega_{k-1}{{\zeta}})}{\omega_{k-1}}&\text{0}&\frac{\text{cos}(\omega_{k-1}{{\zeta}})-\text{1}}{\omega_{k-1}}&\text{0}\\\
\text{0}&\text{cos}(\omega_{k-1}{{\zeta}})&\text{0}&-\text{sin}(\omega_{k-1}{{\zeta}})&\text{0}\\\
\text{0}&\frac{\text{1}-\text{cos}(\omega_{k-1}{{\zeta}})}{\omega_{k-1}}&\text{1}&\frac{\text{sin}(\omega_{k-1}{{\zeta}})}{\omega_{k-1}}&\text{0}\\\
\text{0}&\text{\text{sin}}(\omega_{k-1}{{\zeta}})&\text{0}&\text{cos}(\omega_{k-1}{{\zeta}})&\text{0}\\\
\text{0}&\text{0}&\text{0}&\text{0}&\text{1}\end{bmatrix}\mathbf{x}_{k-\text{1}}$
(77)

where the state vector
$\mathbf{x}_{k}=[a_{k},\dot{{a_{k}}},b_{k},\dot{{b_{k}}},\omega_{k}]^{\top}$
is composed of the 2D position coordinates $({a_{k}},{b_{k}})$, the
corresponding velocities $(\dot{{a_{k}}},\dot{{b_{k}}})$, the angular velocity
$\omega_{k}$ of the target at time instant $k$, ${{\zeta}}$ denotes the
sampling period, and $\mathbf{q}_{k-\text{1}}\sim
N\left(0,\mathbf{Q}_{k-\text{1}}\right)$. $\mathbf{Q}_{k-\text{1}}$ is given
in terms of scaling parameters $\eta_{1}$ and $\eta_{2}$ as [33]

$\mathbf{Q}_{k-\text{1}}=\begin{bmatrix}\eta_{1}\mathbf{M}&0&0\\\
0&\eta_{1}\mathbf{M}&0\\\
0&0&\eta_{2}\end{bmatrix},\mathbf{M}=\begin{bmatrix}{{\zeta}}^{3}/3&{{\zeta}}^{2}/2\\\
{{\zeta}}^{2}/2&{{\zeta}}\end{bmatrix}$

Range readings are obtained using $m$ sensors installed in a zig-zag fashion
as depicted in Fig. 3. The $i$th sensor is located at
$\big{(}a^{\rho_{i}}=350(i-1),b^{\rho_{i}}=350\ ((i-1)\mod 2)\big{)}$ for
$i=1\cdots m$. We assume the first sensor as the common sensor for reference
resulting in $m-1$ TDOA-based measurements. The nominal measurement equation
can be expressed as

$\mathbf{y}_{k}=\mathbf{h}(\mathbf{x}_{k})+\mathbf{r}_{k}$ (78)

with

$\displaystyle{h^{j}(\mathbf{x}_{k})}=$
$\displaystyle\Big{\\{}\sqrt{(a_{k}-a^{\rho_{1}})^{2}+(b_{k}-b^{\rho_{1}})^{2}}$
$\displaystyle-\sqrt{(a_{k}-a^{\rho_{j+1}})^{2}+(b_{k}-b^{\rho_{j+1}})^{2}}\Big{\\}}$
(79)

for $j=1\cdots m-1$. The corresponding nominal covariance measurement matrix
is fully populated given as [21]

$\displaystyle\mathbf{R}_{k}=\begin{bmatrix}{\sigma^{2}_{1}}+{\sigma^{2}_{2}}&\dots&{\sigma^{2}_{1}}\\\
\vdots&\ddots&\vdots\\\
{\sigma^{2}_{1}}&\cdots&{\sigma^{2}_{1}}+{\sigma^{2}_{m}}\end{bmatrix}$ (80)

where $\sigma^{2}_{i}$ is the variance contribution of the $i$th sensor in the
resulting covariance matrix. To consider the effect of outliers the
measurement equation can be modified as

$\mathbf{y}_{k}=\mathbf{h}(\mathbf{x}_{k})+\mathbf{r}_{k}+\mathbf{o}_{k}$ (81)

where $\mathbf{o}_{k}$ produces the effect of outliers in the measurements and
is assumed to obey the following distribution

$\displaystyle p(\mathbf{o}_{k})$
$\displaystyle=\prod_{j=1}^{m-1}\mathcal{J}^{j}_{k}{\mathcal{N}}(o^{j}_{k}|0,\gamma({\sigma^{2}_{1}}+{\sigma^{2}_{j}}))$

where $\mathcal{J}^{j}_{k}$ is a Bernoulli random variable, with values $0$
and $1$, that controls whether an outlier in the $j$th dimension occurs. Let
$\lambda$ denote the probability that a sensor’s TOA measurement is affected.
Therefore, the probability that no outlier appears in the $j$th dimension,
corresponding to $\mathcal{J}^{j}_{k}=0$, is $(1-\lambda)^{2}$ since the first
sensor is a common reference for the TDOA-based measurements. We assume that
the TOA measurements are independently affected and the corruption of the
first TOA observation affects all the measurements. Similarly, the parameter
$\gamma$ controls the variance of an outlier in each dimension respectively.
Using the proposed model we generate the effect of outliers in the data.

For filtering performance comparisons, we choose a hypothetical Gaussian
filter that is a perfect rejector having apriori knowledge of all outlier
instances. We also consider the generalized and independent VBKF estimators
[19], referred from hereon as Gen. VBKF and Ind. VBKF, for comparisons. In
VBKFs, we set the design parameter as $N=1$ (essentially resorting to the EM
method) since higher $N$ results in more computational strain. Lastly, we use
the derived BCRB-based filtering lower bounds to benchmark the performance of
all the filters. Similarly, for smoothing we use the counterparts of all the
considered filters i.e. a perfect outlier-rejecting general Gaussian RTS
smoother and the generalized/independent VBKF-based RTS smoothers denoted as
Gen. VBKS and Ind. VBKS.

For simulations the following values of parameters are used: the initial state
$\mathbf{x}_{0}=[0,1,0,-1,-0.0524]^{\text{T}}$, $\zeta=1$, $\eta_{1}=0.1$,
$\eta_{2}=1.75\times 10^{-4}$, ${\sigma^{2}_{j}}=10$. The initialization
parameters of estimators are:
${\mathbf{m}}^{+}_{0}\sim\mathcal{N}(\mathbf{x}_{0},\mathbf{P}^{+}_{0}$),
$\mathbf{P}^{+}_{0}=\mathbf{Q}_{k}$, $\epsilon=10^{-6}$ and
$\theta^{i}_{k}=0.5~{}\forall~{}i$. For each method we use the Unscented
transform (UT) for approximating the Gaussian integrals [43], in all the
considered methods. Resultingly, the Unscented Kalman Filter (UKF) becomes the
core inferential engine for all the techniques.

In each method, UT parameters are set as $\alpha=1$, $\beta=2$, and
$\kappa=0$. Moreover, we use the same threshold of $10^{-4}$ for the
convergence criterion in each algorithm. Other parameters for VBKFs/VBKSs are
assigned values as originally documented. All the simulations are repeated
with a total time duration $K=400$ and $100$ independent MC runs. Moreover, we
use box and whisker plots to visualize all the results.

###  VI-A Filtering Performance

We assess the relative filtering performance under different scenarios.

Figure 4: MSE vs $\lambda$ $(m=10,\gamma=1000)$

First, we choose $10$ number of sensors with $\gamma=1000$ and increase the
TOA contamination probability $\lambda$. Fig. 4 shows the mean squared error
(MSE) of the state estimate of each filter as $\lambda$ is increased. For
$\lambda=0$ all the filters essentially work as the standard UKF having
similar performance. As $\lambda$ increases, MSE of each method and the lower
bound value are seen to increase. The hypothetical ideal UKF exhibits the best
performance followed by the proposed EMORF, Gen. VBKF, and Ind. VBKF
respectively. The trend remains the same for each $\lambda$. Similar patterns
have been observed for other combinations of $m$ and $\gamma$. Performance
degradation of Ind. VBKF as compared to EMORF and Gen. VBKF is expectable as
it ignores the measurement correlations during filtering. We find EMORF to be
generally more robust in comparison to Gen. VBKF. Our results are not
surprising given that we found the modified selective observation rejecting
(mSOR)-UKF to be more resilient to outliers as compared to the modified
outlier-detecting (mOD)-UKF [18], which are designed for independent
measurements having similar structures to EMORF and VBKF respectively.

Figure 5: MSE vs $m$ ($\lambda=0.3$ and $\gamma=200$)

Next, we vary the number of sensors and assess the estimation performance of
the filters. Fig. 5 shows the MSE of each method as the number of sensors is
increased with $\lambda=0.3$ and $\gamma=200$. As expected, the error bound
and MSE of each filter decrease with increasing number of sensors since more
sources of information become available. We see a pattern similar to the
previous case with the best performance exhibited by the hypothetical ideal
UKF followed by EMORF, Gen. VBKF, and Ind. VBKF respectively. Moreover, we
have observed similar trends for other values of $\lambda$ and $\gamma$ as
well.

Figure 6: Time vs $m$ ($\lambda=0.3$ and $\gamma=200$)

Subsequently, we evaluate the processing overhead of each algorithm by varying
the number of sensors. Fig. 6 shows the execution time taken by each algorithm
as the number of sensors is increased with $\lambda=0.2$ and $\gamma=100$. We
observe that the ideal UKF and Ind. VBKF take lesser time for execution having
a complexity of $\mathcal{O}(m^{3})$. However, EMORF and Gen. VBKF induce more
computational overhead, having a complexity of $\mathcal{O}(m^{4})$, due to
utilization of matrix inverses and determinants for evaluating each of the
$\mathcal{I}^{i}_{k}$ and $\mathbf{z}^{(i)}_{t}$ $\forall\ i=1\cdots m$ in
EMORF and Gen. VBKF respectively. This is the cost we pay for achieving
robustness with correlated measurement noise. Nevertheless, we find that EMORF
generally takes less processing time as compared to Gen. VBKF as shown in Fig.
6. Moreover, similar performance has been observed for other combinations of
$\lambda$ and $\gamma$. This can be attributed to a simpler model being
employed in EMORF resulting in reduced computations.

###  VI-B Smoothing Performance

For smoothing we perform analogous experiments and observe similar
performance.

Figure 7: MSE vs $\lambda$ $(m=15,\gamma=100)$

First, we choose $15$ number of sensors and increase the TOA contamination
probability $\lambda$ with $\gamma=100$. Fig. 7 shows how MSE of the state
estimate of each smoother changes as $\lambda$ is increased. Similar to
filtering, we observe that MSE of each estimator grows with increasing
$\lambda$ including the BCRB-based smoothing lower bound. The hypothetical RTS
smoother performs the best followed by the proposed EMORS, Gen. VBKS, and Ind.
VBKS respectively. The trend remains the same for each $\lambda$. Similar
patterns have been seen for other combinations of $m$ and $\gamma$.

Figure 8: MSE vs $m$ ($\lambda=0.2$ and $\gamma=100$)

Thereafter, we assess the estimation performance of the filters by varying the
number of sensors. Fig. 8 depicts MSE of each estimator as the number of
sensors increase with $\lambda=0.2$ and $\gamma=100$. MSE for each smoother
decreases with growing number of sensors including the BRCB-based lower bound.
The hypothetical RTS smoother is the best performing followed by EMORS, Gen.
VBKS, and Ind. VBKS respectively. We have observed similar trends for other
values of $\lambda$ and $\gamma$ as well.

Figure 9: Time vs $m$ ($\lambda=0.2$ and $\gamma=100$)

Lastly, we evaluate the computational overhead of each algorithm by varying
the number of sensors. Fig. 9 shows the time each method takes as the number
of sensors is increased with $\lambda=0.2$ and $\gamma=100$. We observe
similar patterns as for filtering. The ideal RTS smoother and Ind. VBKS having
a complexity of $\mathcal{O}(m^{3})$ take lesser execution time. EMORS and
Gen. VBKS having a complexity $\mathcal{O}(m^{4})$ are more time consuming.
EMORS generally induces lesser computing overhead as compared to Gen. VBKS as
depicted in Fig. 9. We have observed similar patterns for different
combinations of $\lambda$ and $\gamma$ .

##  VII Conclusion

We consider the problem of outlier-robust state estimation assuming the
existence of measurement noise correlation. Given their advantages, resorting
to tuning-free learning-based approaches is an attractive option in this
regard. Identifying the shortcomings of such existing VB-based tractable
methods, we propose EMORF and EMORS. Since the standard VB approach entails
significant processing complexity, we adopt EM in our algorithmic
constructions. We can conclude that the presented methods are simpler and
hence more practicable as compared to the state-of-the-art Gen. VBKF/VBKS,
devised for the same conditions. This is possible due to the reduction of
inference parameters resulting from the proposal of an uncomplicated model.
Also, the need of the specialized digamma function during implementation is
obviated. In addition, numerical experiments in an illustrative TDOA-based
target tracking example suggest further merits of the proposed methods. We
find that EMORF/EMORS generally exhibit lesser errors as compared to Gen. and
Ind. VBKF/VBKS in different scenarios of the example. Moreover, though the
complexity order of EMORF/EMORS and Gen. VBKF/VBKS is the same, the proposed
estimators are found to be computationally more competitive in general for
different test conditions. These merits make the proposed state estimators
worthy candidates for implementation in relevant scenarios.

##  Appendix A Evaluating
${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$

For evaluating ${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$ we
consider that ${\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$ can easily be
rearranged by swapping rows/columns depending on $\hat{\bm{\mathcal{I}}}_{k}$
as

$\displaystyle{{\bm{\mathfrak{R}}}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})=\begin{bmatrix}\bm{\mathtt{R}}_{k}/{\epsilon}&\mathbf{0}\\\
\mathbf{0}&\hat{\mathbf{R}}_{k}\end{bmatrix}$ (82)

where $\bm{\mathtt{R}}_{k}$ is a sub-matrix with diagonal entries of
${\mathbf{R}}_{k}$. $\hat{\mathbf{R}}_{k}$ contains the rest of the fully
populated submatrix of ${\mathbf{R}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$
corresponding to entries of $\hat{\bm{\mathcal{I}}}_{k}=1$. Inversion of
${\bm{\mathfrak{R}}}_{k}({\bm{\hat{\mathcal{I}}}_{k}})$ results in

$\displaystyle{\bm{\mathfrak{R}}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})=$
$\displaystyle\begin{bmatrix}{\epsilon}\bm{\mathtt{R}}^{-1}_{k}&\mathbf{0}\\\
\mathbf{0}&\hat{\mathbf{R}}^{-1}_{k}\end{bmatrix}$ (83)

Finally, ${\bm{\mathfrak{R}}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$ can be
swapped accordingly to obtain the required matrix
${\mathbf{R}}_{k}^{-1}({\bm{\hat{\mathcal{I}}}_{k}})$ .

##  Appendix B Simplifying $\hat{\tau}^{i}_{k}$

We can swap the ith row/column entries of
${\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}}=1,\hat{\bm{\mathcal{I}}}^{i-}_{k})$
and
${\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}}=1,\hat{\bm{\mathcal{I}}}^{i-}_{k})$
with the first row/column elements to obtain

$\displaystyle|{\bm{\mathfrak{R}}}_{k}({\mathcal{I}^{i}_{k}}=1,\hat{\bm{\mathcal{I}}}^{i-}_{k})|$
$\displaystyle=\begin{vmatrix}R^{i,i}_{k}&\mathbf{R}^{i,-i}_{k}\\\
\mathbf{R}^{-i,i}_{k}&\hat{\mathbf{R}}^{-i,-i}_{k}\end{vmatrix}$ (84)
$\displaystyle|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}}=\epsilon,\hat{\bm{\mathcal{I}}}^{i-}_{k})|$
$\displaystyle=\begin{vmatrix}R^{i,i}_{k}/\epsilon&\mathbf{0}\\\
\mathbf{0}&\hat{\mathbf{R}}^{-i,-i}_{k}\end{vmatrix}$ (85)

Consequently, we can write

$\displaystyle\ln\Big{(}\frac{|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}=1},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}{|{\mathbf{R}}_{k}({\mathcal{I}^{i}_{k}=\epsilon},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}\Big{)}=\ln\Big{(}\frac{|{\bm{\mathfrak{R}}}_{k}({\mathcal{I}^{i}_{k}=1},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}{|{\bm{\mathfrak{R}}}_{k}({\mathcal{I}^{i}_{k}=\epsilon},\hat{\bm{\mathcal{I}}}^{i-}_{k})|}\Big{)}$
$\displaystyle=\ln\Bigg{|}\textbf{I}-\frac{\mathbf{R}^{-i,i}_{k}\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}}{{R}^{i,i}_{k}}\Bigg{|}+\ln(\epsilon)$
(86)

where we have used the following property from matrix algebra [44]

$\displaystyle\begin{vmatrix}\bm{\mathsf{A}}&\bm{\mathsf{B}}\\\
\bm{\mathsf{C}}&\bm{\mathsf{D}}\end{vmatrix}=|\bm{\mathsf{A}}||\bm{\mathsf{D}}-\bm{\mathsf{C}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{B}}|$

Resultingly, we can simplify (23) to (25).

##  Appendix C Evaluating $\triangle{\hat{\mathbf{R}}}_{k}^{-1}$

To avoid redundant calculations during the evaluation of
$\triangle{\hat{\mathbf{R}}}_{k}^{-1}$, we can first swap the ith row/column
elements of matrices with the first row/column entries in (24) to obtain

$\displaystyle\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}$
$\displaystyle=({\bm{\mathfrak{R}}}_{k}^{-1}({\mathcal{I}^{i}_{k}=1},\hat{\bm{\mathcal{I}}}^{i-}_{k})-{\bm{\mathfrak{R}}}_{k}^{-1}({\mathcal{I}^{i}_{k}=\epsilon},\hat{\bm{\mathcal{I}}}^{i-}_{k}))$
$\displaystyle=\begin{bmatrix}R^{i,i}_{k}&\mathbf{R}^{i,-i}_{k}\\\
\mathbf{R}^{-i,i}_{k}&\hat{\mathbf{R}}^{-i,-i}_{k}\end{bmatrix}^{-1}-\begin{bmatrix}R^{i,i}_{k}/\epsilon&\mathbf{0}\\\
\mathbf{0}&\hat{\mathbf{R}}^{-i,-i}_{k}\end{bmatrix}^{-1}$ (87)

To simplify (87), we use the following property from matrix algebra [44]

$\displaystyle\begin{bmatrix}\bm{\mathsf{A}}&\bm{\mathsf{B}}\\\
\bm{\mathsf{C}}&\bm{\mathsf{D}}\end{bmatrix}^{-1}=\begin{bmatrix}\bm{\mathsf{S}}^{-1}&-\bm{\mathsf{S}^{-1}\bm{\mathsf{B}}\bm{\mathsf{D}}^{-1}}\\\
-\bm{\mathsf{D}}^{-1}\bm{\mathsf{C}}\bm{\mathsf{S}^{-1}}&\bm{\mathsf{D}}^{-1}+\bm{\mathsf{D}}^{-1}\bm{\mathsf{C}}\bm{\mathsf{S}}^{-1}\bm{\mathsf{B}}\bm{\mathsf{D}}^{-1}\end{bmatrix}$

where $\bm{\mathsf{S}}$ is the Schur’s complement of $\bm{\mathsf{D}}$ given
as
$\bm{\mathsf{S}}=\bm{\mathsf{A}}-\bm{\mathsf{B}}\bm{\mathsf{D}}^{-1}\bm{\mathsf{C}}$.
As a result, we obtain

$\displaystyle\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}$
$\displaystyle=\begin{bmatrix}\bm{\Xi}^{i,i}&\bm{\Xi}^{i,-i}\\\
\bm{\Xi}^{-i,i}&\bm{\Xi}^{-i,-i}\end{bmatrix}$ (88)

where the expressions of $\bm{\Xi}^{i,i},\bm{\Xi}^{i,-i},\bm{\Xi}^{-i,i}$ and
$\bm{\Xi}^{-i,-i}$ are given in (27)-(30). The redundant calculations in
(27)-(30) can be computed once and stored for further computations e.g.
$({R^{i,i}_{k}-\mathbf{R}^{i,-i}_{k}(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}\mathbf{R}^{-i,i}_{k}})^{-1}$
and $(\hat{\mathbf{R}}^{-i,-i}_{k})^{-1}$. Lastly, the first row/column
entries of $\triangle{\hat{\bm{\mathfrak{R}}}}_{k}^{-1}$ are interchanged to
the actual $i$th row/column positions to obtain the required
$\triangle{\hat{\mathbf{R}}}_{k}^{-1}$.

## References

  * [1] J. Kuti, I. J. Rudas, H. Gao, and P. Galambos, “Computationally relaxed Unscented Kalman filter,” _IEEE Transactions on Cybernetics_ , vol. 53, no. 3, pp. 1557–1565, 2023.
  * [2] B. Lian, Y. Wan, Y. Zhang, M. Liu, F. L. Lewis, and T. Chai, “Distributed Kalman consensus filter for estimation with moving targets,” _IEEE Transactions on Cybernetics_ , vol. 52, no. 6, pp. 5242–5254, 2022.
  * [3] H. Zhang, X. Zhou, Z. Wang, and H. Yan, “Maneuvering target tracking with event-based mixture Kalman filter in mobile sensor networks,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 10, pp. 4346–4357, 2020.
  * [4] A. Donder and F. R. y. Baena, “Kalman-filter-based, dynamic 3-D shape reconstruction for steerable needles with fiber bragg gratings in multicore fibers,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2262–2275, 2022.
  * [5] M. Gardner and Y.-B. Jia, “Pose and motion estimation of free-flying objects: Aerodynamics, constrained filtering, and graph-based feature tracking,” _IEEE Transactions on Robotics_ , vol. 38, no. 5, pp. 3187–3202, 2022.
  * [6] C. Mishra, L. Vanfretti, and K. D. Jones, “Synchrophasor phase angle data unwrapping using an Unscented Kalman filter,” _IEEE Transactions on Power Systems_ , vol. 36, no. 5, pp. 4868–4871, 2021.
  * [7] B. L. Boada, M. J. L. Boada, and H. Zhang, “Sensor fusion based on a dual Kalman filter for estimation of road irregularities and vehicle mass under static and dynamic conditions,” _IEEE/ASME Transactions on Mechatronics_ , vol. 24, no. 3, pp. 1075–1086, 2019.
  * [8] M. Jahja, D. Farrow, R. Rosenfeld, and R. J. Tibshirani, “Kalman filter, sensor fusion, and constrained regression: Equivalences and insights,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [9] A. H. Chughtai, M. Tahir, and M. Uppal, “A robust Bayesian approach for online filtering in the presence of contaminated observations,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–15, 2021.
  * [10] A. H. Chughtai, U. Akram, M. Tahir, and M. Uppal, “Dynamic state estimation in the presence of sensor outliers using MAP-based EKF,” _IEEE Sensors Letters_ , vol. 4, no. 4, pp. 1–4, 2020.
  * [11] R. E. Kalman, “A new approach to linear filtering and prediction problems,” 1960\.
  * [12] M. S. Grewal and A. P. Andrews, _Kalman filtering: Theory and Practice with MATLAB_. John Wiley & Sons, 2014\.
  * [13] S. J. Julier and J. K. Uhlmann, “New extension of the Kalman filter to nonlinear systems,” in _Signal processing, sensor fusion, and target recognition VI_ , vol. 3068. Spie, 1997, pp. 182–193.
  * [14] N. J. Gordon, D. J. Salmond, and A. F. Smith, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation,” in _IEE proceedings F (radar and signal processing)_ , vol. 140, no. 2. IET, 1993, pp. 107–113.
  * [15] G. Evensen, “The ensemble Kalman filter: Theoretical formulation and practical implementation,” _Ocean dynamics_ , vol. 53, pp. 343–367, 2003\.
  * [16] H. E. Rauch, F. Tung, and C. T. Striebel, “Maximum likelihood estimates of linear dynamic systems,” _AIAA journal_ , vol. 3, no. 8, pp. 1445–1450, 1965\.
  * [17] S. Särkkä and L. Svensson, _Bayesian filtering and smoothing_. Cambridge university press, 2023, vol. 17.
  * [18] A. H. Chughtai, M. Tahir, and M. Uppal, “Outlier-robust filtering for nonlinear systems with selective observations rejection,” _IEEE Sensors Journal_ , vol. 22, no. 7, pp. 6887–6897, 2022.
  * [19] H. Li, D. Medina, J. Vilà-Valls, and P. Closas, “Robust variational-based Kalman filter for outlier rejection with correlated measurements,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 357–369, 2021.
  * [20] Y. Bar-Shalom, “Negative correlation and optimal tracking with Doppler measurements,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 37, no. 3, pp. 1117–1120, 2001.
  * [21] R. Kaune, J. Hörst, and W. Koch, “Accuracy analysis for TDOA localization in sensor networks,” in _14th International Conference on Information Fusion_ , 2011, pp. 1–8.
  * [22] Y. Zhu, R. Blum, Z.-Q. Luo, and K. M. Wong, “Unexpected properties and optimum-distributed sensor detectors for dependent observation cases,” _IEEE Transactions on Automatic Control_ , vol. 45, no. 1, pp. 62–72, 2000\.
  * [23] H. Ma and Y. Liu, “Correlation based video processing in video sensor networks,” in _2005 International Conference on Wireless Networks, Communications and Mobile Computing_ , vol. 2, 2005, pp. 987–992 vol.2.
  * [24] K. Yuen, B. Liang, and B. Li, “A distributed framework for correlated data gathering in sensor networks,” _IEEE Transactions on Vehicular Technology_ , vol. 57, no. 1, pp. 578–593, 2008.
  * [25] Y. Huang, Y. Zhang, N. Li, and J. Chambers, “Robust Student’s t based nonlinear filter and smoother,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 52, no. 5, pp. 2586–2596, 2016.
  * [26] H. Wang, H. Li, W. Zhang, and H. Wang, “Laplace l1 robust Kalman filter based on majorization minimization,” in _2017 20th International Conference on Information Fusion (Fusion)_ , 2017, pp. 1–5.
  * [27] C. D. Karlgaard, “Nonlinear regression Huber–Kalman filtering and fixed-interval smoothing,” _Journal of guidance, control, and dynamics_ , vol. 38, no. 2, pp. 322–330, 2015.
  * [28] L. Chang and K. Li, “Unified form for the robust Gaussian information filtering based on M-estimate,” _IEEE Signal Processing Letters_ , vol. 24, no. 4, pp. 412–416, 2017.
  * [29] L. Chang, B. Hu, G. Chang, and A. Li, “Multiple outliers suppression derivative-free filter based on Unscented transformation,” _Journal of guidance, control, and dynamics_ , vol. 35, no. 6, pp. 1902–1906, 2012.
  * [30] S. Wang, W. Gao, and A. P. S. Meliopoulos, “An alternative method for power system dynamic state estimation based on Unscented transform,” _IEEE Transactions on Power Systems_ , vol. 27, no. 2, pp. 942–950, 2012\.
  * [31] A. K. Singh and B. C. Pal, “Decentralized dynamic state estimation in power systems using Unscented transformation,” _IEEE Transactions on Power Systems_ , vol. 29, no. 2, pp. 794–804, 2014.
  * [32] A. Nakabayashi and G. Ueno, “Nonlinear filtering method using a switching error model for outlier-contaminated observations,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 7, pp. 3150–3156, 2020.
  * [33] H. Wang, H. Li, J. Fang, and H. Wang, “Robust Gaussian Kalman filter with outlier detection,” _IEEE Signal Processing Letters_ , vol. 25, no. 8, pp. 1236–1240, 2018.
  * [34] R. Piché, S. Särkkä, and J. Hartikainen, “Recursive outlier-robust filtering and smoothing for nonlinear systems using the multivariate Student-t distribution,” in _2012 IEEE International Workshop on Machine Learning for Signal Processing_ , 2012, pp. 1–6.
  * [35] T. Zhang, S. Zhao, X. Luan, and F. Liu, “Bayesian inference for State-space models with Student-t mixture distributions,” _IEEE Transactions on Cybernetics_ , vol. 53, no. 7, pp. 4435–4445, 2023.
  * [36] K. P. Murphy, _Machine learning: a probabilistic perspective_ , 2012.
  * [37] V. Šmídl and A. Quinn, _The variational Bayes method in signal processing_. Springer Science & Business Media, 2006.
  * [38] A. H. Chughtai, A. Majal, M. Tahir, and M. Uppal, “Variational-based nonlinear Bayesian filtering with biased observations,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 5295–5307, 2022.
  * [39] M. A. A. Turkman, C. D. Paulino, and P. Müller, _Computational Bayesian statistics: an introduction_. Cambridge University Press, 2019, vol. 11.
  * [40] H. L. Van Trees, _Detection, estimation, and modulation theory, part I: detection, estimation, and linear modulation theory_. John Wiley & Sons, 2004.
  * [41] P. Tichavsky, C. Muravchik, and A. Nehorai, “Posterior Cramer-Rao bounds for discrete-time nonlinear filtering,” _IEEE Transactions on Signal Processing_ , vol. 46, no. 5, pp. 1386–1396, 1998.
  * [42] M. Šimandl, J. Královec, and P. Tichavský, “Filtering, predictive, and smoothing Cramér–Rao bounds for discrete-time nonlinear dynamic systems,” _Automatica_ , vol. 37, no. 11, pp. 1703–1716, 2001.
  * [43] E. A. Wan, R. Van Der Merwe, and S. Haykin, “The Unscented Kalman filter,” _Kalman filtering and neural networks_ , vol. 5, no. 2007, pp. 221–280, 2001.
  * [44] F. Zhang, _The Schur complement and its applications_. Springer Science & Business Media, 2006, vol. 4.