#  Outlier-robust Kalman Filtering through Generalised Bayes

Gerardo Duran-Martin  Matias Altamirano  Alexander Y. Shestopaloff  Leandro
Sánchez-Betancourt  Jeremias Knoblauch  Matt Jones  François-Xavier Briol
Kevin Murphy

###### Abstract

We derive a novel, provably robust, and closed-form Bayesian update rule for
online filtering in state-space models in the presence of outliers and
misspecified measurement models. Our method combines generalised Bayesian
inference with filtering methods such as the extended and ensemble Kalman
filter. We use the former to show robustness and the latter to ensure
computational efficiency in the case of nonlinear models. Our method matches
or outperforms other robust filtering methods (such as those based on
variational Bayes) at a much lower computational cost. We show this
empirically on a range of filtering problems with outlier measurements, such
as object tracking, state estimation in high-dimensional chaotic systems, and
online learning of neural networks.

Machine Learning, ICML, filtering, online learning, robust

##  1 Introduction

Probabilistic state-space models (SSMs) are widely used to address problems in
time-series forecasting, online learning, tracking problems, and signal
processing. A key challenge in SSMs is to perform online (sequential)
posterior inference, also known as Bayesian filtering. If the model is linear
and the and the dynamics are known, then the Kalman filter (KF) algorithm is
the optimal filter in terms of minimal mean squared error (Morris, 1976). If
the model is not linear, approximate inference methods must be used. Although
there are a number of techniques available, in this paper we focus on the
extended Kalman filter (EKF) (see e.g., Sarkka & Svensson, 2023) and the
ensemble Kalman filter (EnKF) (see e.g., Roth et al., 2017a), because they
both admit closed-form Bayesian updates to Gaussian posterior approximations
and scale to high dimensions.

Figure 1:  First state component of the SSM (20). The grey dots are
measurements sampled from (20) and the red crosses are measurements sampled
from an outlier measurement process. The dotted blue line shows the KF
posterior mean estimate and the solid orange line shows our proposed WoLF
posterior mean estimate. The regions around the posterior mean cover two
standard deviations. For comparison, the dashed black line shows the true
sampled state process.

In a Bayesian setting, the KF and KF-like methods typically use a Gaussian
measurement model for computational convenience. However, a weakness of this
assumption is that extreme observations, such as outliers, are not well
modelled by Gaussians, leading to model misspecification (Grubbs, 1969; Hampel
et al., 2011). Moreover, in practice, the true dynamics are often different
from the one assumed by the filter. This could be because the measurements are
corrupted, or because the user does not have access to the measurement
function. For example, in tracking problems, sensor measurement errors record
the wrong position of an object; similarly, in economics and finance, there
are instances of data providers sending time series with erroneous values (see
e.g. Liu, 2020, for a survey).

There is a large literature that studies the filtering problem in the presence
of outliers and misspecified measurement models111Alternatively, a measurement
model can be referred to as an observation model.. Some approaches extend the
state space by introducing hierarchical priors on the measurement model. Since
closed-form estimation of the state is not generally feasible, a factorised
variational Bayes (VB) approach can be employed to obtain a fixed-point
solution (see e.g., Ting et al., 2007; Agamennoni et al., 2012; Huang et al.,
2016; Nurminen et al., 2015; Wang et al., 2018; Piché et al., 2012).
Alternatively, Huberised (E)KF-type methods are optimisation-based approaches
to filtering, which minimise the Huber loss; see Boncelet & Dickinson (1983);
Karlgaard (2015); Das (2023) and the references therein. Recently, Boustati et
al. (2020) proposed a particle filter (PF) approach that uses generalised
Bayes (GB) to handle model misspecification from outliers; this builds on work
showing that GB updates are theoretically sound and robust to outliers and
model misspecification (see, e.g., Bissiri et al., 2016; Jewson et al., 2018;
Knoblauch et al., 2022; Fong et al., 2021; Jewson & Rossell, 2022; Husain &
Knoblauch, 2022; Matsubara et al., 2023, 2022). However, these approaches can
be slow, due to the need to perform iterative optimisation per step (for VB)
(Knoblauch et al., 2018), or the need to use a large number of samples (for
PF) (Boustati et al., 2020). See Appendix A for an overview of a number of
these methods.

In this paper, we propose a novel approach that tackles the filtering problem
in the presence of outliers and measurement model misspecification. Our method
is based on the GB approach where one replaces the log-likelihood of the
measurement process with a loss function. We call our method the _weighted
observation likelihood filter (WoLF)_ because it uses a weighted log-
likelihood as loss. A key advantage of this choice is that we can have closed-
form (conjugate) update equations that compute an approximate Gaussian
posterior. We derive WoLF variants for the KF, the EKF, and the ensemble
Kalman Filter (EnKF).

Our approach has several key advantages over prior work: (i) it is fast and
has similar computational cost to the KF thanks to closed-form updates, (ii)
it is flexible to the form of misspecification, (iii) it is provably robust to
outliers, and (iv) it is easy to implement and straightforward to apply to
other filtering methods. See Figure 1 for an illustrative example.

The remainder of the paper proceeds as follows: in Section 2 we introduce the
filtering problem and some common algorithms to approximate a Gaussian
posterior. In Section 3 we present our method and derive WoLF variants to the
KF, the EKF, and the EnKF. Next, in Section 3.4 we show that for certain
choices of weighting functions, our method is provably robust. Finally, in
Section 4, we show empirically that our method matches and in some cases
outperforms previous robust filtering methods at a lower running time. We
consider a variety of filtering problems, including 2d-tracking, state
estimation in high-dimensional chaotic systems, and online learning of neural
networks. Our code can be found at https://github.com/gerdm/weighted-
likelihood-filter.

##  2 Background: Filtering in SSMs

We briefly review Kalman filtering and the main extensions we consider in this
work. In what follows, $m,d\in\mathbb{N}$ are the dimensions of the state and
measurement processes, and $p(\cdot)$ is used for densities. Given an initial
state ${\bm{\theta}}_{0}\in{}^{m}$, a Markovian state-space model (SSM) is
defined as

$\displaystyle{\bm{\theta}}_{t}$
$\displaystyle=f_{t}({\bm{\theta}}_{t-1})+\bm{\phi}_{t},$ (1)
$\displaystyle{\bm{y}}_{t}$
$\displaystyle=h_{t}({\bm{\theta}}_{t})+\bm{\varphi}_{t},$ (2)

for $t\in\\{1,\ldots,T\\}$. Here, ${\bm{\theta}}_{t}\in{}^{m}$ is the (latent)
state vector, ${\bm{y}}_{t}\in{}^{d}$ is the (observed) measurement vector,
$f_{t}:{}^{m}\to{}^{m}$ is the dynamics function, $h_{t}:{}^{m}\to{}^{d}$ is
the measurement function,222 In many applications, the function $h_{t}$ is
modulated by an exogenous feature vector ${\bm{x}}_{t}$. See e.g., Section
4.2.  $\bm{\phi}_{t}$ is a zero-mean Gaussian-distributed random vector with
known covariance matrix $\mathbf{Q}_{t}$, and $\bm{\varphi}_{t}$ is any zero-
mean random vector representing the measurement noise. The meaning of the
state depends on the application; for example, it can be the position of an
object, the state of the atmosphere, or the weights of a neural network, as we
will see in the results section.

The filtering of an SSM consists of two steps: the predict step and the update
step. Given the result of a previous update step (i.e., the iterative
posterior) $p({\bm{\theta}}_{t-1}|{\bm{y}}_{1:t-1})$, the predict step
estimates the prior predictive distribution

$p({\bm{\theta}}_{t}|{\bm{y}}_{1:t-1})=\int
p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})p({\bm{\theta}}_{t-1}|{\bm{y}}_{1:t-1}){\rm
d}{\bm{\theta}}_{t-1},$ (3)

and the update step estimates the new posterior distribution

$p({\bm{\theta}}_{t}|{\bm{y}}_{1:t})\propto
p({\bm{y}}_{t}|{\bm{\theta}}_{t})p({\bm{\theta}}_{t}|{\bm{y}}_{1:t-1}).$ (4)

where $\propto$ indicates equality up to a multiplicative normalisation
constant.

###  2.1 The Kalman filter

Suppose the SSM is linear and Gaussian, that is,

$\displaystyle{\bm{\theta}}_{t}$
$\displaystyle=\mathbf{F}_{t}\,{\bm{\theta}}_{t-1}+\bm{\phi}_{t},$ (5)
$\displaystyle{\bm{y}}_{t}$
$\displaystyle=\mathbf{H}_{t}\,{\bm{\theta}}_{t}+\bm{\varphi}_{t},$

with $\bm{\varphi}_{t}$ a zero-mean Gaussian with known covariance matrix
$\mathbf{R}_{t}$, $\mathbf{F}_{t}\in{}^{m\times m}$, and
$\mathbf{H}_{t}\in{}^{d\times m}$. Then, given the initial condition
$p({\bm{\theta}}_{0})={\cal
N}({\bm{\theta}}_{0}\,|\,\bm{\mu}_{0},\,\bm{\Sigma}_{0})$ — a Gaussian density
with known mean $\bm{\mu}_{0}\in\mathbb{R}^{m}$ and covariance
$\bm{\Sigma}_{0}\in\mathbb{R}^{m\times m}$ — the exact Bayesian predict and
update steps are given by

$\displaystyle p({\bm{\theta}}_{t}|{\bm{y}}_{1:t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\bm{\mu}_{t|t-1},\,\bm{\Sigma}_{t|t-1}),$ (6)
$\displaystyle p({\bm{\theta}}_{t}|{\bm{y}}_{1:t})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\bm{\mu}_{t},\,\bm{\Sigma}_{t}),$ (7)

with prior predictive mean and covariance given by

$\displaystyle\bm{\mu}_{t|t-1}$
$\displaystyle=\mathbf{F}_{t}\,\bm{\mu}_{t-1},$ (8)
$\displaystyle\bm{\Sigma}_{t|t-1}$
$\displaystyle=\mathbf{F}_{t}\,\bm{\Sigma}_{t-1}\,\mathbf{F}_{t}^{\intercal}+\mathbf{Q}_{t},$

and posterior mean and covariance given by

$\displaystyle\bm{\Sigma}_{t}^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\intercal}\,\mathbf{R}_{t}^{-1}\,\mathbf{H}_{t},$
(9) $\displaystyle\mathbf{K}_{t}$
$\displaystyle=\bm{\Sigma}_{t}\,\mathbf{H}_{t}^{\intercal}\,\mathbf{R}_{t}^{-1},$
$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}\,({\bm{y}}_{t}-\hat{\bm{y}}_{t}),$

where $\hat{{\bm{y}}}_{t}=\mathbf{H}_{t}\,\bm{\mu}_{t|t-1}$ is the predicted
observation and $\mathbf{K}_{t}$ is the Kalman gain matrix used to map the
error (residual) vector in observation space to an update in latent state
space.

###  2.2 The Extended Kalman filter

When the state and measurement functions are non-linear, a common approach is
to introduce a Gaussian posterior density
$q({\bm{\theta}}_{t}|{\bm{y}}_{1:t})={\cal
N}({\bm{\theta}}_{t}\,|\,\bm{\mu}_{t},\,\bm{\Sigma}_{t})$ that approximates
the posterior density $p({\bm{\theta}}_{t}|{\bm{y}}_{1:t})$ through a
Kullback-Leibler projection. This is done using Gaussian approximations of the
joint densities $q({\bm{\theta}}_{t},{\bm{\theta}}_{t-1}|{\bm{y}}_{1:t-1})$
and $q({\bm{\theta}}_{t},{\bm{y}}_{t}|{\bm{y}}_{1:t-1})$ and then performing a
closed-form Bayesian update; see Section 8.4 in Sarkka & Svensson (2023) for
details.

The EKF is a special case of the above in which one linearises $f_{t}$ in (1)
and $h_{t}$ in (2), and assumes a Gaussian measurement noise
$\bm{\varphi}_{t}$. As a consequence, the predict and update equations
resemble those of the standard KF. More precisely, the EKF replaces the mean
of the state transition density $p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ in
(1) with

$\mathbb{E}[{\bm{\theta}}_{t}|{\bm{\theta}}_{t-1}]\approx\mathbf{F}_{t}({\bm{\theta}}_{t-1}-\bm{\mu}_{t-1})+f_{t}(\bm{\mu}_{t-1})=:\bar{\bm{\mu}}_{t|t-1},$
(10)

and the measurement mean of $p({\bm{y}}_{t}|{\bm{\theta}}_{t})$ in (2) with

$\mathbb{E}[{\bm{y}}_{t}|{\bm{\theta}}_{t}]\approx\mathbf{H}_{t}({\bm{\theta}}_{t}-\bm{\mu}_{t|t-1})+h_{t}(\bm{\mu}_{t|t-1})=:\bar{{\bm{y}}}_{t},$
(11)

where
$\bm{\mu}_{t|t-1}=\mathbb{E}[\bar{\bm{\mu}}_{t|t-1}|\bm{\mu}_{t-1}]=f_{t}(\bm{\mu}_{t-1})$,
$\mathbf{F}_{t}$ is the Jacobian of $f_{t}$ evaluated at $\bm{\mu}_{t-1}$,
$\mathbf{H}_{t}$ is the Jacobian of $h_{t}$ evaluated at $\bm{\mu}_{t|t-1}$,
and
$\hat{{\bm{y}}}_{t}=\mathbb{E}[\bar{{\bm{y}}}_{t}]=h_{t}(\bm{\mu}_{t|t-1})$.
The linearisation of the transition and measurement around the respective
previous means allows the resulting predict and update equations to closely
resemble (8) and (9). Therefore, the algorithm remains relatively scalable and
more efficient low-rank extensions can also be derived (see e.g. Chang et al.,
2023; Lambert et al., 2023; Cartea et al., 2023).

###  2.3 The ensemble Kalman filter

The ensemble Kalman filter (EnKF) was developed as an alternative to the
extended Kalman filter for high-dimensional and highly non-linear state
dynamics (see e.g., Evensen, 1994; Burgers et al., 1998; Roth et al., 2017a).
This method avoids the need to compute Jacobians and the storage of an
explicit $m\times m$ posterior covariance matrix, by instead representing the
belief state with an ensemble of $N\in\mathbb{N}$ particles
$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}\in{}^{m}$ for $i=1,\ldots,N$ evolved
through time. This enables the method to scale to high-dimensional state
spaces with complex nonlinear dynamics, such as those frequently arising in
data assimilation for weather forecasting (Evensen, 2009). In the EnKF, the
predict step samples (1) to obtain
$\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}$. Then, the update step samples
predictions $\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}\in{}^{d}$, for each
particle, according to

$\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}\sim\mathcal{N}\left(h_{t}\left(\hat{{\bm{\theta}}}_{t|t-1}^{(i)}\right),\mathbf{R}_{t}\right).$
(12)

The particles are then updated according to

$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+\bar{\mathbf{K}}_{t}\big{(}{\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\big{)},$
(13)

with gain matrix $\bar{\mathbf{K}}_{t}$ calculated from the ensemble. See
Appendix D.3 for details.

##  3 The weighted observation likelihood filter

Our method is based on the GB approach where one modifies the update step in
(4) to use a loss function $\ell_{t}:{}^{m}\to\real$ in place of the negative
log-likelihood of the measurement process. This gives the generalised
posterior

$q({\bm{\theta}}_{t}|{\bm{y}}_{1:t})\propto\exp(-\ell_{t}({\bm{\theta}}_{t}))q({\bm{\theta}}_{t}|{\bm{y}}_{1:t-1}).$
(14)

We propose to gain robustness to outliers in observation space by taking the
loss function to be the model’s negative log-likelihood scaled by a data-
dependent weighting term

$\ell_{t}({\bm{\theta}}_{t})=-W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})\,\log
q({\bm{y}}_{t}|{\bm{\theta}}_{t}),$ (15)

with $W:{}^{d}\times{}^{d}\to\real$ the weighting function and
$q({\bm{y}}_{t}|{\bm{\theta}}_{t})$ the modelled measurement process. We call
our method the _weighted observation likelihood filter (WoLF)_. To specify an
instance of our method, ones needs to define the likelihood
$q({\bm{y}}_{t}|{\bm{\theta}}_{t})$ and the weighting function $W$. In the
next subsections, we show the flexibility of WoLF and derive weighted-
likelihood-based KF, EKF, and EnKF algorithms. Setting
$W({\bm{y}}_{t},\bar{{\bm{y}}}_{t})=1$ trivially recovers existing methods,
but we will instead use non-constant weighting functions inspired by the work
of (Barp et al., 2019; Matsubara et al., 2022; Altamirano et al., 2023a, b).

###  3.1 Linear weighted observation likelihood filter

The following proposition gives a closed-form solution for the update step of
WoLF under a linear measurement function and a Gaussian likelihood (see
Section C.1 for the proof).

######  Proposition 3.1.

Consider the linear-Gaussian SSM (5) with weighting function
$W:{}^{d}\times{}^{d}\to\real$. Then, the update step of WoLF with loss
function (15) is given by (9) with $\mathbf{R}_{t}^{-1}$ replaced by
$\bar{\mathbf{R}}_{t}^{-1}=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})\,\mathbf{R}_{t}^{-1}$.

The resulting predict and update steps for WoLF under linear dynamics and
zero-mean Gaussians for the state and measurement process are shown in
Algorithm 1.

0: $\mathbf{F}_{t}$, $\mathbf{Q}_{t}$ // predict step

$\bm{\mu}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\mu}_{t-1}$

$\bm{\Sigma}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\Sigma}_{t-1}\,\mathbf{F}_{t}^{\intercal}+\mathbf{Q}_{t}$

${\bm{y}}_{t}$, $\mathbf{H}_{t}$, $\mathbf{R}_{t}$ // update step

$\hat{{\bm{y}}}_{t}\leftarrow\mathbf{H}_{t}\,\bm{\mu}_{t|t-1}$

$w_{t}\leftarrow W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$

$\bm{\Sigma}_{t}^{-1}\leftarrow\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{2}\,\mathbf{H}_{t}^{\intercal}\,\mathbf{R}_{t}^{-1}\,\mathbf{H}_{t}$

$\mathbf{K}_{t}\leftarrow
w_{t}^{2}\,\bm{\Sigma}_{t}\,\mathbf{H}_{t}^{\intercal}\,\mathbf{R}_{t}^{-1}$

$\bm{\mu}_{t}\leftarrow\bm{\mu}_{t|t-1}+\mathbf{K}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\right)$

Algorithm 1 WoLF predict and update step

The computational complexity of WoLF under linear dynamics matches that of the
KF, i.e., $O(m^{3})$. Alternative robust filtering algorithms require multiple
iterations per measurement to achieve robustness and stability, making them
significantly slower; see Table 1 for the computational complexity for the
methods we consider, and Table 2 and Figure 5 for empirical comparisons.

###  3.2 Nonlinear weighted observation likelihood filter

Our method readily extends to other nonlinear filtering algorithms. For
example, a WoLF version of the EKF is obtained by introducing a weighting
function to $\eqref{eq:linearised-measurement-mean}$ yielding the approximate
log-likelihood

$\log
q({\bm{y}}_{t}|{\bm{\theta}}_{t})=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})\log{\cal
N}({\bm{y}}_{t}\,|\,\bar{{\bm{y}}}_{t},\,\mathbf{R}_{t}).$ (16)

Similarly, we can derive a weighted ensemble Kalman filter by weighting error
terms in the update step (13); see Appendix D.3 for details. Finally, our
method can be easily extended to handle measurement processes modelled as an
exponential family, such as in classification; see Appendix D.1 for details.

###  3.3 The choice of weighting function

Weighted likelihoods have a well-established history in Bayesian inference and
have demonstrated their efficacy in improving robustness (Grünwald, 2012;
Holmes & Walker, 2017; Grünwald & van Ommen, 2017; Miller & Dunson, 2018;
Bhattacharya et al., 2019; Alquier & Ridgway, 2020; Dewaskar et al., 2023). In
this context, the corresponding posteriors are often referred to as
fractional, tempered, or power posteriors. In most existing work, the
determination of weights relies on heuristics and the assigned weights remain
constant across all data points so that
$W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})=w\in\real$ for all $t$. In contrast, we
dynamically incorporate information from the most recent observations without
incurring additional computational costs by defining the weight as a function
of the current observation ${\bm{y}}_{t}$ and its prediction
$\hat{{\bm{y}}}_{t}=h_{t}(\bm{\mu}_{t|t-1})$, which is based on all of the
past observations.

To define the weighting function, we take inspiration from previous work for
dealing with outliers. In particular, Wang et al. (2018) proposed classifying
robust filtering algorithms into two main types: compensation-based
algorithms, which incorporate information from tail events into the model in a
robust way (see, e.g., Huang et al., 2016; Agamennoni et al., 2012), and
detect-and-reject algorithms, which assume that outlier observations bear no
useful information and thus are ignored (see, e.g., Wang et al., 2018; Mu &
Yuen, 2015). Below we show how both of these strategies can be implemented
using our WoLF method by merely changing the weighting function.

##### Inverse multi-quadratic weighting function:

As an example of a compensation-based method, we follow Altamirano et al.
(2023b) and use the Inverse Multi-Quadratic (IMQ) weighting, which in our SSM
setting is

$\displaystyle W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$
$\displaystyle=\left(1+\frac{||{\bm{y}}_{t}-\hat{{\bm{y}}}_{t}||_{2}^{2}}{c^{2}}\right)^{-1/2},$
(17)

where $c>0$ is the soft threshold and $\|\cdot\|$ denotes the $l_{2}$ norm. We
call WoLF with IMQ weighting “WoLF-IMQ”.

##### Mahalanobis-based weighting function:

The $l_{2}$ norm in the IMQ can be modified to account for the covariance
structure of the measurement process by replacing it with the Mahalanobis
distance between ${\bm{y}}_{t}$ and $\hat{{\bm{y}}}_{t}$:

$\displaystyle W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$
$\displaystyle=\left(1+\frac{\|\mathbf{R}_{t}^{-1/2}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{2}^{2}}{c^{2}}\right)^{-1/2}.$
(18)

We call WoLF with this weighting function the WoLF-MD method. This type of
weighted IMQ function has been used extensively in the kernel literature (see
e.g. Chen et al., 2019; Detommaso et al., 2018; Riabiz et al., 2022).

##### Threshold Mahalanobis-based weighting function:

As an example of a detect-and-reject method, we consider

$W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})=\begin{cases}1&\text{if
}\|\mathbf{R}_{t}^{-1/2}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{2}^{2}\leq c\\\
0&\text{otherwise}\end{cases}$ (19)

with $c>0$ the fixed threshold. The weighting function (19) corresponds to
ignoring information from estimated measurements whose Mahalanobis distance to
the true measurement is larger than some predefined threshold $c$. In the
linear setting, this weighting function is related to the benchmark method
employed in Ting et al. (2007). We refer to WoLF with this weighting function
as “WoLF-TMD”.

For ${\bm{y}}_{t}\in{}^{d}$ with $d=m\gg 1$ and diagonal measurement
covariance $\mathbf{R}_{t}=\text{diag}\left(r_{t,1},\ldots,r_{t,d}\right)$,
the WoLF-TMD function can be modified to weight individual observations so
that $W:{}^{d}\times{}^{d}\to{}^{d}$. See Section 4.3 for an example and
Appendix D.2 for a discussion.

The proposed weighting functions — the IMQ, the MD, and the TMD — are defined
such that $W:{}^{d}\times{}^{d}\to[0,1]$ and therefore can only down-weight
observations. This means that our updates are always conservative, i.e., our
posteriors will be wider in the presence of outliers (see Figure 1 for an
example).

###  3.4 Theoretical properties

In this section, we prove the outlier-robustness for WoLF-type methods. We use
the classical framework of Huber (1981). Consider measurements
${\bm{y}}_{1:t}$. We measure the influence of a contamination
${\bm{y}}^{c}_{t}$ by examining the divergence between the posterior with the
original observation ${\bm{y}}_{t}$ and the posterior with the contamination
${\bm{y}}_{t}^{c}$, which is allowed to be arbitrarily large. As a function of
${\bm{y}}_{t}^{c}$, this divergence is called the _posterior influence
function_ (PIF) and was studied in Ghosh et al. (2016); Matsubara et al.
(2022); Altamirano et al. (2023a, b). Following Altamirano et al. (2023b), we
consider the Kullback-Leibler (KL) divergence, which allows us to obtain
closed-form expressions for Gaussians. The PIF is given by

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})=\operatorname{KL}\left(p({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})\|p({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})\right).$

If
$\sup_{{\bm{y}}_{t}^{c}\in\mathbb{R}^{d}}|\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})|<\infty$,
then the posterior is called outlier-robust, which indicates that as
$\|{\bm{y}}_{t}-{\bm{y}}_{t}^{c}\|_{2}\to\infty$, the contamination’s effect
on the posterior is bounded. This is the Bayesian equivalent to bias-
robustness in frequentist statistics.

######  Theorem 3.2.

Consider the linear Gaussian SSM (5) or its linearised (EKF) approximation.
The standard (E)KF posterior has an unbounded PIF and is not outlier robust.

In contrast, the generalised posterior presented in Proposition 3.1 has
bounded PIF and is, therefore, outlier robust for any weighting function $W$
such that
$\sup_{{\bm{y}}_{t}\in{}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})<\infty$ and
$\sup_{{\bm{y}}_{t}\in{}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\,\|{\bm{y}}_{t}\|_{2}<\infty$.

The proof is in Section C.2. In particular, the conditions are satisfied when
$W$ is (17), (18), or (19), which are the focus of our paper. Figure 2 shows
an empirical validation of this proposition for the 2D tracking problem
detailed in Section 4.1. Here, it is clear that the PIF for the standard KF is
unbounded, whereas our methods exhibit a bounded influence. In Appendix C.3 we
show a version of the theorem above for the EnKF.

Figure 2:  PIF for the 2d tracking problem of Section 4.1. The last
measurement ${\bm{y}}_{t}$ is replaced with
${\bm{y}}^{c}_{t}={\bm{y}}_{t}+\bm{\epsilon}$, where
$\bm{\epsilon}\in[-5,5]\times[-5,5]$. We observe that the PIF is asymetric for
the weighted methods; this is because the weighting term is a function of the
prior predictive and the measurement at time $t$. See Appendix E.1 for a more
detailed explanation.

##  4 Experiments

In this section, we study the performance of the WoLF methods in multiple
filtering settings. Each experiment employs a dataset (or samples data from an
SSM), a collection of benchmark methods, and a metric to compare the methods.

For our robust baselines, we make use of three methods that are representative
of recent state-of-the-art approaches to robust filtering: the Bernoulli KF of
Wang et al. (2018) ( KF-B), which is an example of a detect-and-reject
strategy; the inverse-Wishart filter of Agamennoni et al. (2012) ( KF-IW),
which is an example of a compensation-based strategy; and the Huberised EnKF
of Roh et al. (2013) ( Hub-EnKF), which is an example of a Huberised
algorithm. The  KF-B and  KF-IW are deterministic and optimise a VB objective
to compute a Gaussian approximation to the state posterior (see Appendix A for
details). We do not compare against sophisticated hierarchical methods nor
methods based on particle filtering, because these do not scale well to high-
dimensional state spaces. For the neural network fitting problem, we also
consider a variant of online gradient descent ( OGD) based on Adam (Kingma &
Ba, 2017), which uses multiple inner iterations per step (measurement). This
method does scale to high-dimensional state spaces, but sadly only gives a
maximum a posteriori (MAP) estimate and is not as sample efficient as a robust
Bayesian filter.

For experiments where KF or EKF is used as the baseline, we consider the
following WoLF variants: (i) the WoLF version with inverse multi-quadratic
weighting function ( WoLF-IMQ), (ii) the thresholded WoLF with Mahalanobis-
based weighting function ( WoLF-TMD). When using the EKF variants, we
linearise the state mean (10) and measurement mean (11). For experiments where
the ensemble KF ( EnKF) is taken as the baseline algorithm, we benchmark the
performance of the weighted likelihood EnKF with (i) weighting with averaged-
particles (AP-EnKF) and (ii) the per-particle weightings ( PP-EnKF). See
Appendix D.3 for a detailed description of the robust EnKF methods.

Method | Cost | #HP | Ref
---|---|---|---
KF | $O(m^{3})$ | 0 | Kalman
KF-B | $O(I\,m^{3})$ | 3 | Wang et al.
KF-IW | $O(I\,m^{3})$ | 2 | Agamennoni et al.
OGD | $O(I\,m^{2})$ | 2 | Bencomo et al.
WoLF-IMQ | $O(m^{3})$ | 1 | (Ours)
WoLF-TMD | $O(m^{3})$ | 1 | (Ours)
Table 1:  Computational complexity of the update step, assuming $d\leq m$ and
assuming linear dynamics. Here, $I$ is the number of inner iterations, #HP
refers to the number of hyperparameters we tune, and ”Cost” refers to the
computational complexity.

In each experiment, and unless otherwise specified, we run $100$ trials to
evaluate each method. The hyperparameters of each method are chosen on the
first trial using the Bayesian optimisation (BO) package of Nogueira (2014).
BO is a popular derivative-free approach to function maximisation (see e.g.
Frazier, 2018). Specifically, we optimise the hyperparameters that minimise
the chosen metric on the first run of each experiment. Where a multi-output
metric is specified, the minimisation is taken over the maximum of the output.
The hyperparameters for KF/EKF-like methods are: the noise scaling and number
of inner iterations for the  KF-IW, the two shape parameters and the number of
inner iterations for the  KF-B, the learning rate and number of inner
iterations for the  OGD, the thresholding value for the  WoLF-TMD, and the
soft threshold hyperparameter for the  WoLF-IMQ. See Table 1 for a summary.

The results we obtain can be summarised as follows: WoLF-based methods either
outperform or match the performance of their counterparts in the metrics we
specify below, but typically at a much lower running time.

###  4.1 Robust KF for tracking a 2D object

We consider the classical problem of estimating the position of an object
moving in 2D with constant velocity, which is commonly used to benchmark
tracking problems (see e.g., Example 8.2.1.1 in Murphy (2023) or Example 4.5
in Sarkka & Svensson (2023)). The SSM takes the form

$\displaystyle p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\mathbf{F}_{t}{\bm{\theta}}_{t-1},\,\mathbf{Q}_{t}),$
(20) $\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t})$ $\displaystyle={\cal
N}({\bm{y}}_{t}\,|\,\mathbf{H}_{t}{\bm{\theta}}_{t},\,\mathbf{R}_{t}),$

where $\mathbf{Q}_{t}=q\,{\bf I}_{4}$, $\mathbf{R}_{t}=r\,{\bf I}_{2}$,
$({\bm{\theta}}_{0,t},{\bm{\theta}}_{1,t})$ is the position,
$({\bm{\theta}}_{2,t},{\bm{\theta}}_{3,t})$ is the velocity,

$\displaystyle\mathbf{F}_{t}$ $\displaystyle=\begin{pmatrix}1&0&\Delta&0\\\
0&1&0&\Delta\\\ 0&0&1&0\\\ 0&0&0&1\end{pmatrix},$
$\displaystyle\mathbf{H}_{t}$ $\displaystyle=\begin{pmatrix}1&0&0&0\\\
0&1&0&0\end{pmatrix},$

$\\!\\!\Delta=0.1$ is the sampling rate, $q=0.10$ is the system noise, $r=10$
is the measurement noise, and ${\bf I}_{K}$ is a $K\times K$ identity matrix.
We simulate 500 trials, each with 1,000 steps. For each method, we compute the
scaled RMSE metric
$J_{T,i}=\sqrt{\sum_{t=1}^{T}({\bm{\theta}}_{t,i}-\bm{\mu}_{t,i})^{2}}$ for
$i\in\\{0,1,2,3\\}$ as well as the total running time (relative to the  KF).

In our experiments, the true data generating process is one of two variants of
(20). The first variant (which we call Student observations) corresponds to a
system whose measurement process comes from the Student-t likelihood:

$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t})$
$\displaystyle=\text{St}({{\bm{y}}_{t}\,|\,\mathbf{H}_{t}{\bm{\theta}}_{t},\,\mathbf{R}_{t},\nu_{t}})$
(21) $\displaystyle=\int_{0}^{\infty}{\cal
N}\left({\bm{y}}_{t}\,|\,\mathbf{H}_{t}{\bm{\theta}}_{t},\frac{\mathbf{R}_{t}}{\tau}\right)\text{Gam}\left(\tau|\frac{\nu_{t}}{2},\frac{\nu_{t}}{2}\right)d\tau,$

with $\text{Gam}(\cdot|a,b)$ the density of a Gamma distribution with shape
$a$ and rate $b$, and $\nu_{t}=2.01$. The second variant (which we call
mixture observations) corresponds to a system where the mean of the
observations changes sporadically. Instances of this variant can occur as a
form of human error or a software bug in a data-entry program. To emulate this
scenario, we modify (20) by using the following mixture model for the
observation process:

$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t})$ $\displaystyle={\cal
N}({\bm{y}}_{t}\,|\,{\bm{m}}_{t},\,\mathbf{R}_{t}),$ (22)
$\displaystyle{\bm{m}}_{t}$
$\displaystyle=\begin{cases}\mathbf{H}_{t}\,{\bm{\theta}}_{t}&\text{w.p.}\
1-p_{\epsilon},\\\ 2\,\mathbf{H}_{t}\,{\bm{\theta}}_{t}&\text{w.p.}\
p_{\epsilon},\end{cases}$

where $p_{\epsilon}=0.05$.

##### Results

Figure 3:  The left panel shows a sample path using the Student variant and
the right panel shows a sample path using the mixture variant. The top left
figure on each panel shows the true underlying state in black, and the
measurements as grey dots.

Figure 3 shows a sample of each variant along with the filtered state for each
method. For the Student variant (left panel), the  WoLF-IMQ and the  WoLF-TMD
estimate the true state more closely than the competing methods. Both the  KF-
IW and the  KF-B look comparable to the  KF, which are not robust to outliers.
For the mixture variant (right panel),333 The top left figure in the right
panel is cropped, see Figure 13 for the uncropped version.  we see that the
WoLF-IMQ, the  WoLF-TMD, and the  KF-B filter the true state correctly. In
contrast, the  KF-IW and the  KF are not robust to outliers.444 KF-B removes
outliers that bear no information according to some criterion, but in the
Student-t case, it fails.  KF-IW, on the other hand, estimates a measurement
covariance rather than the dispersion of such a measurement covariance. In
this sense, it is misspecified in both cases.

Figure 4:  Distribution (across 500 2d tracking trials) of RMSE for first
component of the state vector, $J_{T,0}$. Left panel: Student observation
model. Right panel: Mixture observation model.

The results in Figure 3 hold for multiple trials as shown in Figure 4, which
plots the distribution of the errors in the first component of the state
vector. The heavy tail of the  WoLF-TMD on the mixture observations variant,
and the distribution of $J_{T,k}$ for all $k$ are studied in Appendix E.2.1.

Method | Student | Mixture
---|---|---
KF-B | 2.0x | 3.7x
KF-IW | 1.2x | 5.3x
WoLF-IMQ (ours) | 1.0x | 1.0x
WoLF-TMD (ours) | 1.0x | 1.0x
Table 2: Mean slowdown rate over  KF.

Table 2 shows the median slowdown (in running time) to process the
measurements relative to the  KF. The slowdown for method X is obtained diving
the running time of method X over the running time of the  KF. Under the
Student variant, the  WoLF-IMQ, the  WoLF-TMD, and the  KF-IW have similar
running time to the  KF. In contrast, the  KF-B takes twice the amount of
time. Under the mixture variant, the  KF-B and the  KF-IW are almost four
times and five times slower than the  KF respectively. The changes in slowdown
rate are due the number of inner iterations that were chosen during the first
trial.

###  4.2 Robust EKF for online MLP regression (UCI)

In this section, we benchmark the methods using a corrupted version of the
tabular UCI regression datasets.555The dataset is available at
https://github.com/yaringal/DropoutUncertaintyExps. Similar to other papers
that deal with non-linear online learning (see, e.g. Chang et al., 2023), we
consider a single-hidden-layer multi-layered perceptron (MLP) with twenty
hidden units and a real-valued output unit. In this experiment, the state
dimension (number of parameters in the MLP) is $m=(n_{\text{in}}\times
20+20)+(20\times 1+1)$, where $n_{\text{in}}$ is the dimension of the feature
${\bm{x}}_{t}$.666 See Table 3 for the values that $n_{\text{in}}$ takes for
each dataset.  One of the main advantages of using a Bayesian filtering method
for fitting neural networks (compared to using  OGD) is the ability to handle
non-stationary distributions (see e.g. Chang et al., 2023). Below, we take a
static state (see e.g. Lambert et al., 2021), so that the prior predictive
mean is $\bm{\mu}_{t|t-1}=\bm{\mu}_{t-1}$. In Appendix E.3 we study online
learning with non-stationary environments.

Each trial is carried out as follows: first, we randomly shuffle the rows in
the dataset; second, we divide the dataset into a warmup dataset (10% of rows)
and a corrupted dataset (remaining 90% of rows); third, we normalise the
corrupted dataset using min-max normalisation from the warmup dataset; fourth,
with probability $p_{\epsilon}=0.1$, we replace a measurement
${\bm{y}}_{t}\in\real$ with a corrupted data point ${\bm{u}}_{t}\sim{\cal
U}[-50,50]$; and fifth, we run each method on the corrupted dataset.

For each dataset and for each method, we evaluate the prior predictive
$\text{RMedSE}=\sqrt{\text{median}\\{({\bm{y}}_{t}-h_{t}(\bm{\mu}_{t|t-1}))^{2}\\}_{t=1}^{T}}$,
which is the squared root of the median squared error between the measurement
${\bm{y}}_{t}$ and the prior predictive
$h_{t}(\bm{\mu}_{t|t-1})=h(\bm{\mu}_{t|t-1},{\bm{x}}_{t})$.777We use median
instead of mean because we have outliers in measurement space. Here, $h$ is
the MLP. We also evaluate the average time step of each method, i.e., we run
each method and divide the total running time by the number of samples in the
corrupted dataset.

##### Results

Figure 5 shows the percentage change of the RMedSE and the percentage change
of running time with respect to those of the  OGD for all corrupted UCI
datsets.

Figure 5:  RMedSE versus time per step (relative to the  OGD minus $1$) across
the corrupted UCI datasets.

Given the computational complexity of the remaining methods, ideally, a robust
Bayesian alternative to the  OGD should be as much to the left as possible on
the $x$-axis (rel. time step) and as low as possible on the $y$-axis (rel.
RMedSE). We observe that the  WoLF-IMQ and the  WoLF-TMD have both of these
traits. In particular, we observe that the only two points in the third
quadrant are those of the  WoLF-IMQ and the  WoLF-TMD. Note that the  EKF-IW
and the  EKF-B have much higher relative running time and the  EKF has much
higher relative RMedSE.

###  4.3 Robust EnKF for Lorenz96 model

We consider a modified version of the Lorenz96 model that is commonly used to
simulate the atmosphere (see e.g. Lorenz, 2006; Arnold et al., 2013)

For a fixed $\Delta>0$, the SSM is given by

$\displaystyle\frac{{\bm{\theta}}_{t+\Delta,i}-{\bm{\theta}}_{t,i}}{\Delta}$
$\displaystyle=\Big{(}{\bm{\theta}}_{t,i+1}-{\bm{\theta}}_{t,i-2}\Big{)}{\bm{\theta}}_{t,i-1}-{\bm{\theta}}_{t,i}+\bm{\phi}_{t,i},$
(23) $\displaystyle{\bm{y}}_{t,i}$
$\displaystyle=\begin{cases}{\bm{\theta}}_{t,i}+\bm{\varphi}_{t,i}&\text{w.p.
}1-p_{\epsilon},\\\ 100&\text{w.p. }p_{\epsilon}.\end{cases}$

Here, ${\bm{\theta}}_{t,k}$ is the value of the state component $k$ at step
$t$, $\bm{\phi}_{t,i}\sim{\cal N}(8,1)$, $\bm{\varphi}_{t,i}\sim{\cal
N}(0,1)$, $p_{\epsilon}=0.001$, $i=1,\ldots,d$, $t=1,\ldots,T$, with $T\gg 1$
the number of steps, and we use the convention
${\bm{\theta}}_{t,d+k}={\bm{\theta}}_{t,k}$,
${\bm{\theta}}_{t,-k}={\bm{\theta}}_{t,d-k}$. Similar to Roth et al. (2017a),
we integrate the state process in (23) to match the formulation in (1) using
the Runge-Kutta-4 (RK4) procedure with discretisation step $\Delta=0.05$,
integrated over $T=10^{3}$ steps, $N=1,000$ number of particles, and
$d=m=100$. A run of the state process is shown in Figure 6. Note that the
probability of an outlier happening on any state component at any timestep is
$p_{\epsilon}\times d=0.1$.

Figure 6:  Sample of the Lorenz process with $d=40$. Here, the “waves” move
westward. The red dots represent where measurement outliers occur. We take
$p_{\epsilon}=0.01$ and $\Delta=0.05$ integrated over $T=60$ steps.

In this experiment,  EnKF is the baseline. As in Roth et al. (2017b), we use
the metric
$L_{t}=\sqrt{\frac{1}{d}({\bm{\theta}}_{t}-\bm{\mu}_{t})^{\intercal}({\bm{\theta}}_{t}-\bm{\mu}_{t})}$
to measure the in-state RMSE.

##### Results

An evaluation of $L_{t}$ for the  EnKF, the AP-EnKF, the  PP-EnKF, and the
Hub-EnKF is shown in the top row of Figure 7. The grey vertical lines denote
timesteps where an outlier event happened, i.e., at least one entry of
${\bm{y}}_{t}$ is $100$. The top plot shows that the  Hub-EnKF and AP-EnKF
have an almost-identical behaviour. However, the bottom row of Figure 7 shows
that the AP-EnKF and the  PP-EnKF are more robust to the choice of threshold
$c$ compared to the  Hub-EnKF. This is because the  Hub-EnKF makes updates
with outlier observations clipped at $c$, whereas WoLF-like methods disregard
error measurements above $c$.

Figure 7:  The top row shows a sample run of the  EnKF, the AP-EnKF, the  PP-
EnKF, and the  Hub-EnKF; outlier events are shown in grey vertical bars. The
bottom row shows the bootstrap estimate of $L_{T}$ over 20 runs and 500
bootstrapped samples as a function of the $c$ hyperparameter.

We present the results of the EnKF when the number of particles is less than
the number of state components in Appendix E.2.3. In the experiment, we modify
the algorithms to make use of the covariance inflation correction first
proposed by Anderson (1999). The conclusions from this section extend to the
covariance inflation case.

##  5 Conclusion

We introduced a provably robust filtering algorithm based on generalised Bayes
which we call the weighted observation likelihood filter or WoLF. Our
algorithm is as fast as the KF, has closed-form update equations, and is
straightforward to apply to various filtering methods. The superior
performance of the WoLF is shown on a wide range of filtering problems. In
contrast, alternative robust methods either have higher computational
complexity than the WoLF, or similar computational complexity but not higher
performance.

Future work will investigate how to overcome the limitations of our approach.
For example, (i) not being robust to outliers in the state-process, (ii) the
assumption of a known covariance $\mathbf{R}_{t}$, (iii) the assumption of
known dynamics for the covariance $\mathbf{Q}_{t}$, and (iv) the assumption of
a unimodal posterior.

## Acknowledgements

JK and FXB were supported by the EPSRC grant [EP/Y011805/1].

## Impact statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none
which we feel must be specifically highlighted here.

## References

  * Agamennoni et al. (2012) Agamennoni, G., Nieto, J. I., and Nebot, E. M.  Approximate inference in state-space models with heavy-tailed noise.  _IEEE Transactions on Signal Processing_ , 2012.
  * Agee et al. (1980) Agee, W. S., Dunn, B. A., and NM, W. S. M. R.  Robust filtering and smoothing via gaussian mixtures.  _NTIS, SPRINGFIELD, VA, 1980, 34_ , 1980.
  * Alquier & Ridgway (2020) Alquier, P. and Ridgway, J.  Concentration of tempered posteriors and of their variational approximations.  2020\.
  * Altamirano et al. (2023a) Altamirano, M., Briol, F.-X., and Knoblauch, J.  Robust and scalable Bayesian online changepoint detection.  In _International Conference on Machine Learning_ , 2023a.
  * Altamirano et al. (2023b) Altamirano, M., Briol, F.-X., and Knoblauch, J.  Robust and conjugate Gaussian process regression.  _arXiv:2311.00463_ , 2023b.
  * Anderson (1999) Anderson, J. L.  A monte carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts.  _Monthly weather review_ , 127(12):2741–2758, 1999.
  * Arnold et al. (2013) Arnold, H., Moroz, I., and Palmer, T.  Stochastic parametrizations and model uncertainty in the Lorenz’96 system.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 371(1991):20110479, 2013.
  * Arroyo et al. (2024) Arroyo, A., Cartea, A., Moreno-Pino, F., and Zohren, S.  Deep attentive survival analysis in limit order books: Estimating fill probabilities with convolutional-transformers.  _Quantitative Finance_ , pp. 1–23, 2024.
  * Barp et al. (2019) Barp, A., Briol, F.-X., Duncan, A. B., Girolami, M., and Mackey, L.  Minimum Stein discrepancy estimators.  In _Neural Information Processing Systems_ , pp. 12964–12976, 2019.
  * Bencomo et al. (2023) Bencomo, G. M., Snell, J. C., and Griffiths, T. L.  Implicit maximum a posteriori filtering via adaptive optimization.  _arXiv:2311.10580_ , 2023.
  * Bhattacharya et al. (2019) Bhattacharya, A., Pati, D., and Yang, Y.  Bayesian fractional posteriors.  2019\.
  * Bishop (2006) Bishop, C. M.  _Pattern Recognition and Machine Learning (Information Science and Statistics)_.  Springer-Verlag, 2006.
  * Bissiri et al. (2016) Bissiri, P. G., Holmes, C. C., and Walker, S. G.  A general framework for updating belief distributions.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 78(5):1103–1130, 2016.
  * Boncelet & Dickinson (1983) Boncelet, C. G. and Dickinson, B. W.  An approach to robust Kalman filtering.  In _IEEE Conference on Decision and Control_ , pp. 304–305, 1983.
  * Boustati et al. (2020) Boustati, A., Akyildiz, O. D., Damoulas, T., and Johansen, A.  Generalised Bayesian filtering via sequential monte carlo.  _Advances in neural information processing systems_ , 33:418–429, 2020.
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q.  JAX: composable transformations of Python+NumPy programs, 2018.  URL http://github.com/google/jax.
  * Burgers et al. (1998) Burgers, G., Van Leeuwen, P. J., and Evensen, G.  Analysis scheme in the ensemble Kalman filter.  _Monthly weather review_ , 126(6):1719–1724, 1998.
  * Calvet et al. (2015) Calvet, L. E., Czellar, V., and Ronchetti, E.  Robust filtering.  _Journal of the American Statistical Association_ , 110(512):1591–1606, 2015.
  * Cao et al. (2022) Cao, W., Liu, C., Lan, Z., Li, S. E., Pan, W., and Alessandri, A.  Robust bayesian inference for moving horizon estimation.  _arXiv preprint arXiv:2210.02166_ , 2022.
  * Cartea et al. (2023) Cartea, Á., Duran-Martin, G., and Sánchez-Betancourt, L.  Detecting toxic flow.  _arXiv:2312.05827_ , 2023.
  * Chang et al. (2023) Chang, P. G., Duran-Martin, G., Shestopaloff, A. Y., Jones, M., and Murphy, K.  Low-rank extended Kalman filtering for online learning of neural networks from streaming data.  In _Conference on Lifelong Learning Agents_ , 2023.
  * Chen et al. (2019) Chen, W. Y., Barp, A., Briol, F.-X., Gorham, J., Girolami, M., Mackey, L., and Oates, C.  Stein point markov chain Monte Carlo.  In _International Conference on Machine Learning_ , pp. 1011–1021, 2019.
  * Das (2023) Das, S.  _Robust state estimation methods for robotics applications_.  PhD thesis, West Virginia University, 2023.
  * Detommaso et al. (2018) Detommaso, G., Cui, T., Marzouk, Y., Spantini, A., and Scheichl, R.  A Stein variational newton method.  _Advances in Neural Information Processing Systems_ , 2018.
  * Dewaskar et al. (2023) Dewaskar, M., Tosh, C., Knoblauch, J., and Dunson, D. B.  Robustifying likelihoods by optimistically re-weighting data.  _arXiv:2303.10525_ , 2023.
  * Duran-Martin et al. (2022) Duran-Martin, G., Kara, A., and Murphy, K.  Efficient online Bayesian inference for neural bandits.  In _International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Evensen (1994) Evensen, G.  Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics.  _Journal of Geophysical Research: Oceans_ , 99(C5):10143–10162, 1994.
  * Evensen (2009) Evensen, G.  _Data Assimilation: The Ensemble Kalman Filter_.  Springer, 2nd ed. 2009 edition edition, 2009.
  * Fong et al. (2021) Fong, E., Holmes, C., and Walker, S. G.  Martingale posterior distributions.  _arXiv:2103.15671_ , 2021.
  * Frazier (2018) Frazier, P. I.  A tutorial on Bayesian optimization.  _arXiv:1807.02811_ , 2018.
  * Ghosh et al. (2016) Ghosh, S., Delle Fave, F. M., and Yedidia, J.  Assumed density filtering methods for learning Bayesian neural networks.  In _Conference on Artificial Intelligence and Interactive Digital Entertainment_ , 2016.
  * Grubbs (1969) Grubbs, F. E.  Procedures for detecting outlying observations in samples.  _Technometrics_ , pp. 1–21, 1969.
  * Grünwald (2012) Grünwald, P.  The safe Bayesian.  In _International Conference on Algorithmic Learning Theory_ , pp. 169–183, 2012.
  * Grünwald & van Ommen (2017) Grünwald, P. and van Ommen, T.  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  _Bayesian Analysis_ , 12(4):1069 – 1103, 2017.
  * Hampel et al. (2011) Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., and Stahel, W. A.  _Robust statistics: the approach based on influence functions_.  John Wiley & Sons, 2011.
  * Holmes & Walker (2017) Holmes, C. C. and Walker, S. G.  Assigning a value to a power likelihood in a general Bayesian model.  _Biometrika_ , 104(2):497–503, 2017.
  * Horn & Johnson (2012) Horn, R. A. and Johnson, C. R.  _Matrix analysis_.  Cambridge university press, 2012.
  * Huang et al. (2016) Huang, Y., Zhang, Y., Li, N., and Chambers, J.  A robust Gaussian approximate filter for nonlinear systems with heavy tailed measurement noises.  In _IEEE International Conference on Acoustics, Speech and Signal Processing_ , pp. 4209–4213, 2016.
  * Huber (1981) Huber, P. J.  Robust statistics.  _Wiley Series in Probability and Mathematical Statistics_ , 1981.
  * Husain & Knoblauch (2022) Husain, H. and Knoblauch, J.  Adversarial interpretation of Bayesian inference.  In _International Conference on Algorithmic Learning Theory_ , pp. 553–572. PMLR, 2022.
  * Jewson & Rossell (2022) Jewson, J. and Rossell, D.  General Bayesian loss function selection and the use of improper models.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , pp. 1640–1665, 2022.
  * Jewson et al. (2018) Jewson, J., Smith, J. Q., and Holmes, C.  Principled Bayesian minimum divergence inference.  _Entropy_ , pp. 442, 2018.
  * Kalman (1960) Kalman, R. E.  A new approach to linear filtering and prediction problems.  _Journal of Basic Engineering_ , 82(Series D):35–45, 1960.
  * Karlgaard (2015) Karlgaard, C. D.  Nonlinear regression Huber–Kalman filtering and fixed-interval smoothing.  _Journal of guidance, control, and dynamics_ , 38(2):322–330, 2015.
  * Khan & Rue (2023) Khan, M. E. and Rue, H.  The bayesian learning rule.  _Journal of Machine Learning Research_ , 24(281):1–46, 2023.  URL http://jmlr.org/papers/v24/22-0291.html.
  * Kingma & Ba (2017) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization, 2017.
  * Knoblauch et al. (2018) Knoblauch, J., Jewson, J. E., and Damoulas, T.  Doubly robust bayesian inference for non-stationary streaming data with $beta$-divergences.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Knoblauch et al. (2022) Knoblauch, J., Jewson, J., and Damoulas, T.  An optimization-centric view on Bayes’ rule: Reviewing and generalizing variational inference.  _Journal of Machine Learning Research_ , 23(132):1–109, 2022.
  * Lambert et al. (2021) Lambert, M., Bonnabel, S., and Bach, F.  The recursive variational Gaussian approximation (R-VGA).  _Statistics and Computing._ , 32(1):10, 2021.
  * Lambert et al. (2023) Lambert, M., Bonnabel, S., and Bach, F.  The limited-memory recursive variational Gaussian approximation (L-RVGA).  _Statistics and Computing_ , 33(70), 2023.  doi: 10.1007/s11222-023-10239-x.  URL https://inria.hal.science/hal-03501920.
  * Lin et al. (2019) Lin, W., Khan, M. E., and Schmidt, M.  Fast and simple Natural-Gradient variational inference with mixture of exponential-family approximations.  In _ICML_ , June 2019.  URL http://arxiv.org/abs/1906.02914.
  * Liu (2020) Liu, G.  Data quality problems troubling business and financial researchers: A literature review and synthetic analysis.  _Journal of Business & Finance Librarianship_, 25(3-4):315–371, 2020.
  * Lorenz (2006) Lorenz, E. N.  _Predictability – a problem partly solved_ , pp. 40–58.  Cambridge University Press, 2006.
  * Masreliez (1975) Masreliez, C.  Approximate non-gaussian filtering with linear state and observation relations.  _IEEE Transactions on Automatic Control_ , 20(1):107–110, 1975.
  * Masreliez & Martin (1977) Masreliez, C. and Martin, R.  Robust bayesian estimation for the linear model and robustifying the kalman filter.  _IEEE transactions on Automatic Control_ , 22(3):361–371, 1977.
  * Matsubara et al. (2022) Matsubara, T., Knoblauch, J., Briol, F.-X., and Oates, C. J.  Robust generalised Bayesian inference for intractable likelihoods.  _Journal of the Royal Statistical Society: Series B_ , 84(3):997–1022, 2022.
  * Matsubara et al. (2023) Matsubara, T., Knoblauch, J., Briol, F.-X., and Oates, C. J.  Generalized Bayesian inference for discrete intractable likelihood.  _Journal of the American Statistical Association_ , pp. 1–11, 2023.
  * Meyr & Spies (1984) Meyr, H. and Spies, G.  The structure and performance of estimators for real-time estimation of randomly varying time delay.  _IEEE transactions on acoustics, speech, and signal processing_ , 32(1):81–94, 1984.
  * Miller & Dunson (2018) Miller, J. W. and Dunson, D. B.  Robust Bayesian inference via coarsening.  _Journal of the American Statistical Association_ , 2018.
  * Morris (1976) Morris, J.  The Kalman filter: A robust estimator for some classes of linear quadratic problems.  _IEEE Transactions on Information Theory_ , 22(5):526–534, 1976.
  * Mu & Yuen (2015) Mu, H.-Q. and Yuen, K.-V.  Novel outlier-resistant extended Kalman filter for robust online structural identification.  _Journal of Engineering Mechanics_ , 141(1):04014100, 2015.
  * Murphy (2023) Murphy, K. P.  _Probabilistic Machine Learning: Advanced Topics_.  MIT Press, 2023.
  * Nogueira (2014) Nogueira, F.  Bayesian Optimization: Open source constrained global optimization tool for Python, 2014.  URL https://github.com/fmfn/BayesianOptimization.
  * Nurminen et al. (2015) Nurminen, H., Ardeshiri, T., Piché, R., and Gustafsson, F.  Robust inference for state-space models with skewed measurement noise.  _IEEE Signal Processing Letters_ , 22(11):1898–1902, 2015.
  * Ollivier (2018) Ollivier, Y.  Online natural gradient as a Kalman filter.  _Electronic Journal of Statistics_ , 12(2):2930–2961, 2018.
  * Piché et al. (2012) Piché, R., Särkkä, S., and Hartikainen, J.  Recursive outlier-robust filtering and smoothing for nonlinear systems using the multivariate student-t distribution.  In _IEEE International Workshop on Machine Learning for Signal Processing_ , pp. 1–6, 2012.
  * Riabiz et al. (2022) Riabiz, M., Chen, W., Cockayne, J., Swietach, P., Niederer, S. A., Mackey, L., and Oates, C. J.  Optimal thinning of MCMC output.  _arXiv:2005.03952_ , 2022.
  * Roh et al. (2013) Roh, S., Genton, M. G., Jun, M., Szunyogh, I., and Hoteit, I.  Observation quality control with a robust ensemble Kalman filter.  _Monthly Weather Review_ , 141(12):4414–4428, 2013.
  * Roth et al. (2017a) Roth, M., Hendeby, G., Fritsche, C., and Gustafsson, F.  The ensemble Kalman filter: a signal processing perspective.  _Journal on Advances in Signal Processing_ , 2017:1–16, 2017a.
  * Roth et al. (2017b) Roth, M., Hendeby, G., Fritsche, C., and Gustafsson, F.  The ensemble Kalman filter: a signal processing perspective.  _EURASIP Journal on Advances in Signal Processing_ , pp. 1–16, 2017b.
  * Sarkka & Svensson (2023) Sarkka, S. and Svensson, L.  _Bayesian Filtering and Smoothing (2nd edition)_.  Cambridge University Press, 2023.
  * Schick & Mitter (1994) Schick, I. C. and Mitter, S. K.  Robust recursive estimation in the presence of heavy-tailed observation noise.  _The annals of statistics_ , 22(2):1045–1080, 1994.
  * Ting et al. (2007) Ting, J.-A., Theodorou, E., and Schaal, S.  Learning an outlier-robust kalman filter.  In _European Conference on Machine Learning_ , pp. 748–756. Springer, 2007.
  * Wang et al. (2018) Wang, H., Li, H., Fang, J., and Wang, H.  Robust Gaussian Kalman filter with outlier detection.  _IEEE Signal Processing Letters_ , 25(8):1236–1240, 2018.
  * West (1981) West, M.  Robust sequential approximate bayesian estimation.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 43(2):157–166, 1981.

## Supplementary Materials

The Appendix is structured as follows. Appendix A provides an overview of
outlier-robust filtering methods. Appendix B derives the WoLF-MD method from a
maximum-a-posteriori (MAP) estimate perspective. Appendix C collects proofs;
more precisely, Appendix C.1 proves Proposition 3.1 and Appendix C.2 proves
Theorem 3.2. We also show robustness for the IMQ, MD, and TMD weighting
functions. Next, in Appendix D, we discuss the exponential-family, multi-
output weighting function, and EnKF extensions to the WoLF. Finally, in
Appendix E, we investigate the 2d tracking problem of Section 4.1 in more
detail; we conduct robustness checks and we introduce an additional experiment
for corrupted non-linear non-stationary learning.

##  Appendix A  Background on existing robust filters and related methods

Robust Bayesian filtering and the minimum variance estimator dates back to
Masreliez (1975) and Masreliez & Martin (1977). These methods propose a
modified KF-recursion for the linear SSM written in terms of the score
function of the measurement prior predictive. In these early works, inference
relies on Monte Carlo.

West (1981) follows these earlier works and proposes KF-like updates for non-
normal measurement models. That paper also studies whether several popular
likelihood functions are robust in the sense of “ignoring outliers” — the
analysis includes the Student-t, power exponential, Huber, logistic, and
stable-law likelihoods. For the case of Student-t likelihood with one degree
of freedom (i.e., Cauchy) and linear dynamics, the update equation for the
posterior mean derived in West (1981) is equivalent to ours when using IMQ
weights. However, the approach taken in West (1981) cannot recover our TMD
scheme, due to being tied to a given choice of measurement model. Furthermore,
West (1981) does not provide a theoretical foundation to use their approach in
non-linear measurement models.

The work in Meyr & Spies (1984) proposes a scheme to eliminate observations
that a KF procedure labels as outliers. Their methodology relies on a
“secondary decision system” which checks for discrepancies between the
predicted mean and the observation, eliminating observations with high
discrepancies. This scheme is analogous to the TMD scheme, however, it is not
shown to be provably robust.

Another alternative to robustify measurement models against outliers is the
work in Agee et al. (1980). Their paper introduces Gaussian-mixture models for
robust filtering and smoothing. Their inference method is based on particle
filtering.

To the best of our knowledge, the first work that proposes robust filtering in
the context of robust statistics is Calvet et al. (2015). Their robust filter
follows Masreliez & Martin (1977) and is based on a “Huberisation” of the
derivative of the log-measurement density (score function), which then they
integrate. See Schick & Mitter (1994) for a comprehensive review of classical
robust-KF methods for linear SSMs. Recently proposed provably-robust methods
include the work by Boustati et al. (2020) and Cao et al. (2022).

Finally, we discuss similarities and differences with the Bayesian learning
rule (BLR) line of work by Khan & Rue (2023). BLR also allows the replacing of
the likelihood with a loss function, but that loss depends only on the
parameter and the observation, i.e., BLR uses
$\ell_{t}({\bm{\theta}})=L({\bm{y}}_{t},f_{t}({\bm{\theta}}))$, with
$L:{}^{d}\times{}^{d}$ a loss function. By contrast, our loss depends on the
weight $W_{t}$, which in turn depends on the belief state, via
$\hat{{\bm{y}}}_{t}=f_{t}(\bm{\mu}_{t|t-1})$, which in turn is a function of
all the past data, i.e., our loss has the form
$\ell_{t}({\bm{\theta}})=L({\bm{y}}_{1:t},f_{t}({\bm{\theta}}))$. In addition,
WoLF is designed for the online inference setting, whereas BLR, at least in
its standard form, is designed for offline inference. Finally, note that the
exponential family extension of BLR Lin et al. (2019) has also been proposed
as a way to get robustness by using scale mixture posteriors. However, this
induces robustness in parameter space, whereas we focus on robustness in
observation space. In particular, the weighting term for the BLR mixture model
depends on the distance in parameter space, ${\bm{\theta}}_{t}-\bm{\mu}_{t}$,
whereas ours depends on the difference in observation space,
${\bm{y}}_{t}-\hat{{\bm{y}}_{t}}$.

###  A.1 Variational-based methods

In this section, we provide an overview of variational-Bayes (VB) robust
filtering methods. As above, ${\bm{\theta}}_{t}$ is the state vector of
interest and $\bm{\Psi}_{t}$ are additional state parameters. Given the SSM
(1) and measurement model $p({\bm{y}}_{t}|{\bm{\theta}}_{t},\bm{\Psi}_{t})$,
VB-based methods seek an approximate posterior distribution over the extended
state process $\bm{\Phi}_{t}=({\bm{\theta}}_{t},\bm{\Psi}_{t})$ that
factorises as

$q(\bm{\Phi}_{t})=\prod_{k=1}^{K}q(\bm{\Phi}_{t,k}),$ (24)

with $\bm{\Phi}_{t}=(\bm{\Phi}_{t,1},\ldots,\bm{\Phi}_{t,K})$, and $K$ the
number of collections. It can be shown that the log-density $q^{*}$ that
minimises the KL divergence between the true posterior distribution and the
variational distribution is given by

$\log q^{*}(\bm{\Phi}_{t,k})=\mathbb{E}_{\neg k}\left[\log
p({\bm{y}}_{t},\bm{\Phi}_{t})\right]+C,$ (25)

where $C$ is the normalising constant of $q^{*}$, and the notation
$\mathbb{E}_{\neg k}\left[\cdot\right]$ denotes the conditional expectation
given all elements in $\bm{\Phi}_{t}$ except from $\bm{\Phi}_{t,k}$. See
Section 10.1.1 in Bishop (2006) for details. Below, we discuss the robust VB-
based filtering variants we use in the paper.

###  A.2 KF-IW method of Agamennoni et al. (2012)

Agamennoni et al. (2012) extend the state-space to be
$\bm{\Psi}_{t}=({\bm{\theta}}_{t},\mathbf{R}_{t})$, where $\mathbf{R}_{t}$ is
the measurement covariance. Note that the classical KF setting,
$\mathbf{R}_{t}$ is known. The SSM is of the form

$\displaystyle p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\mathbf{F}_{t}{\bm{\theta}}_{t-1},\,\mathbf{Q}_{t}),$
(26) $\displaystyle p(\mathbf{R}_{t})$ $\displaystyle={\cal
W}^{-1}(\mathbf{R}_{t}\,|\,\nu\bm{\Lambda},\nu),$ $\displaystyle
p({\bm{y}}_{t}|{\bm{\theta}}_{t},\mathbf{R}_{t})$ $\displaystyle={\cal
N}({\bm{y}}_{t}\,|\,\mathbf{H}_{t}{\bm{\theta}}_{t},\,\mathbf{R}_{t}),$

where ${\cal W}^{-1}(\cdot|{\bf P},\eta)$ is the density of an inverse Wishart
distribution with positive-definite scale matrix ${\bf P}\in{}^{m\times m}$,
$\eta>m-1$ degrees of freedom, and $\nu>0$ is the noise-scaling
hyperparameter. They consider the class of variational distributions

$q(\bar{{\bm{\theta}}},\bar{\mathbf{R}})=q({\bm{\theta}}_{0})q(\mathbf{R}_{0})\prod_{t=1}^{T}q({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})q(\mathbf{R}_{t}),$
(27)

with $\bar{{\bm{\theta}}}=({\bm{\theta}}_{0},\ldots,{\bm{\theta}}_{T})$ and
$\bar{\mathbf{R}}=(\mathbf{R}_{0},\ldots,\mathbf{R}_{T})$. They show that the
class of VB posteriors (25), under the model in (26) and (27), take the form

$\displaystyle q({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\bm{\mu}_{t},\,\bm{\Sigma}_{t}),$ (28) $\displaystyle
q(\mathbf{R}_{t})$ $\displaystyle={\cal
W}^{-1}(\mathbf{R}_{t}|\nu_{t}\bm{\Lambda}_{t},\nu_{t}),$ (29)

with $\bm{\mu}_{t}$, $\bm{\Sigma}_{t}$, $\bm{\Lambda}_{t}$, and $\nu_{t}$
specified in Algorithm 2.

0: $\mathbf{F}_{t}$, $\mathbf{Q}_{t}$, $\bm{\mu}_{t-1}$, $\bm{\Sigma}_{t-1}$
// predict step

$\bm{\mu}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\mu}_{t-1}$

$\bm{\Sigma}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\Sigma}_{t-1}\,\mathbf{F}_{t}^{\intercal}+\mathbf{Q}_{t}$

$\bm{\mu}_{t},\bm{\Sigma}_{t}\leftarrow\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1}$

${\bm{y}}_{t}$, $\mathbf{H}_{t}$, $\mathbf{R}_{t}$, $\ell\in{}_{+}$ // update
step

for $i=1,\ldots,I$ do

$\mathbf{S}_{t}\leftarrow({\bm{y}}_{t}-\mathbf{H}_{t}\,\bm{\mu}_{t})({\bm{y}}_{t}-\mathbf{H}_{t}\,\bm{\mu}_{t})^{\intercal}+\mathbf{H}_{t}^{\intercal}\,\bm{\Sigma}_{t}\,\mathbf{H}_{t}$

$\bm{\Lambda}_{t}\leftarrow(\ell+1)^{-1}(\ell\,\mathbf{R}_{0}+\mathbf{S}_{t})$

$\mathbf{K}_{t}\leftarrow(\mathbf{H}_{t}\,\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}+\bm{\Lambda}_{t})^{-1}\mathbf{H}_{t}^{\intercal}\,\bm{\Sigma}_{t|t-1}$

$\bm{\mu}_{t}\leftarrow\bm{\mu}_{t|t-1}+\mathbf{K}_{t}^{\intercal}({\bm{y}}_{t}-\mathbf{H}_{t}^{\intercal}\,\bm{\mu}_{t|t-1})$

$\bm{\Sigma}_{t}\leftarrow\mathbf{K}_{t}^{\intercal}\,\bm{\Lambda}_{t}\,\mathbf{K}_{t}+(\mathbf{I}-\mathbf{H}_{t}\,\mathbf{K}_{t})^{\intercal}\,\bm{\Sigma}_{t|t-1}(\mathbf{I}-\mathbf{H}_{t}\,\mathbf{K}_{t})$

end for

Algorithm 2 Agamennoni et al. (2012) predict and update step for i.i.d. noise
with $I\geq 1$ inner iterations.

In this method the hyperparameters are the number of iterations $I$ and the
scaling term $\ell$. The prior measurement covariance is $\mathbf{R}_{0}$.

###  A.3 KF-B method of Wang et al. (2018)

Wang et al. (2018) extend the state-space to be
$\bm{\Psi}_{t}=({\bm{\theta}}_{t},\xi_{t},w_{t})$, where $w_{t}$ is an outlier
event and $\xi_{t}$ is its probability. The SSM is of the form:

$\displaystyle p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\mathbf{F}_{t}{\bm{\theta}}_{t-1},\,\mathbf{Q}_{t}),$
(30) $\displaystyle p(\xi_{t})$
$\displaystyle=\text{Beta}(\xi_{t}|\alpha_{0},\beta_{0}),$ $\displaystyle
p(\rho_{t}|\xi_{t})$ $\displaystyle=\text{Bern}(\rho_{t}|\xi_{t}),$
$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t},\mathbf{R}_{t}\rho_{t})$
$\displaystyle=\begin{cases}{\cal
N}({\bm{y}}_{t}\,|\,h_{t}({\bm{\theta}}_{t}),\,\mathbf{R}_{t})&\text{if
}\rho_{t}=1,\\\ 1&\text{if }\rho_{t}=0.\end{cases}$

In the above equations, $\text{Beta}(\cdot|a,b)$ is the density of a Beta
distribution with shape parameters $a$ and $b$, and $\text{Bern}(\cdot|\pi)$
is the mass of a Bernoulli random variable with parameter $\pi\in[0,1]$. They
consider the class of variational distributions

$q(\bar{{\bm{\theta}}},\bar{\rho},\bar{\xi})=q({\bm{\theta}}_{0})q(\rho_{0})q(\xi_{0})\prod_{t=1}^{T}q({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})q(\rho_{t})q(\xi_{t}),$
(31)

with $\bar{{\bm{\theta}}}=({\bm{\theta}}_{0},\ldots,{\bm{\theta}}_{T})$,
$\bar{\rho}=(\rho_{0},\ldots,\rho_{T})$, $\bar{\xi}=(\xi_{0},\ldots,\xi_{T})$.
We provide the the predict and update equations in Algorithm 3.

0: $\alpha_{0},\beta_{0}$, $\bm{\mu}_{t-1}$, $\bm{\Sigma}_{t-1}$

0: $\mathbf{F}_{t}$, $\mathbf{Q}_{t}$ // predict step

$\bm{\mu}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\mu}_{t-1}$

$\bm{\Sigma}_{t|t-1}\leftarrow\mathbf{F}_{t}\,\bm{\Sigma}_{t-1}\,\mathbf{F}_{t}^{\intercal}+\mathbf{Q}_{t}$

$\bm{\mu}_{t},\bm{\Sigma}_{t}\leftarrow\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1}$

${\bm{y}}_{t}$, $\mathbf{H}_{t}$, $\mathbf{R}_{t}$, $\text{tol.}\ll 1$ //
update step

$\rho_{t},\alpha^{\prime},\beta^{\prime}\leftarrow 1,\alpha_{0},\beta_{0}$

for $i=1,\ldots,I$ do

if $\rho_{\epsilon}<\text{tol.}$ then

$\bm{\mu}_{t}\leftarrow\bm{\mu}_{t|t-1}$

$\bm{\Sigma}_{t}\leftarrow\bm{\Sigma}_{t|-1}$

else

$\bar{\mathbf{R}}_{t}\leftarrow\mathbf{R}/\rho_{t}$

$\hat{{\bm{y}}}_{t}\leftarrow\mathbf{H}_{t}\,\bm{\mu}_{t|t-1}$

$\bm{\Sigma}_{t}^{-1}\leftarrow\bm{\Sigma}_{t|t-1}^{-1}+\,\mathbf{H}_{t}^{\top}\,\bar{\mathbf{R}}_{t}^{-1}\,\mathbf{H}_{t}$

$\mathbf{K}_{t}\leftarrow\,\bm{\Sigma}_{t}\,\mathbf{H}_{t}^{\intercal}\,\bar{\mathbf{R}}_{t}^{-1}$

$\bm{\mu}_{t}\leftarrow\bm{\mu}_{t|t-1}+\mathbf{K}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\right)$

end if

$\mathbf{B}_{t}\leftarrow\mathbb{E}_{{\bm{\theta}}\sim{\cal
N}(\bm{\mu}_{t},\bm{\Sigma}_{t})}\left[({\bm{y}}_{t}-h_{t}({\bm{\theta}})({\bm{y}}_{t}-h_{t}({\bm{\theta}}))^{\intercal}\right]$

$\log\bar{\pi}_{t}\leftarrow\Psi(\alpha^{\prime})-\Psi(\alpha^{\prime}+\beta^{\prime}+1)$

$\log(1-\bar{\pi}_{t})\leftarrow\Psi(\beta^{\prime}+1)-\Psi(\alpha^{\prime}+\beta^{\prime}+1)$

$\rho_{t}\leftarrow\frac{\exp(\log\bar{\pi}_{t}-\text{Tr}(\mathbf{B}_{t}\mathbf{R}_{t}^{-1})/2)}{\exp(\log\bar{\pi}_{t}-\text{Tr}(\mathbf{B}_{t}\mathbf{R}_{t}^{-1})/2)+\exp(\log(1-\bar{\pi}_{t}))}$

$\alpha^{\prime}\leftarrow\alpha_{0}+\rho_{t}$

$\beta^{\prime}\leftarrow\beta_{0}+1-\rho_{t}$

end for

Algorithm 3 Wang et al. (2018) predict and update step with $I\geq 1$ inner
iterations.

Here, the hyperparameters are the prior rates $\alpha_{0}$ and $\beta_{0}$,
and the number of inner iterations $I$, $\Psi(\cdot)$ is the digamma function,
and $\mathbf{B}_{t}$ is of closed form after linearising the measurement
function.

###  A.4 Ting et al. (2007)

Ting et al. (2007) extend the state-space to be
$\bm{\Psi}_{t}=({\bm{\theta}}_{t},w_{t})$, where $w_{t}$ is a weighting term
for the observation covariance $\mathbf{R}$. In their method, $\mathbf{R}$ is
known and fixed. The SSM is of the form:

$\displaystyle p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,\mathbf{F}{\bm{\theta}}_{t-1},\,\mathbf{Q}),$ (32)
$\displaystyle p(w_{t})$ $\displaystyle={\rm Gam}(w_{t}|a_{w},b_{w}),$
$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t})$ $\displaystyle={\cal
N}({\bm{y}}_{t}\,|\,\mathbf{H}{\bm{\theta}}_{t},\,\mathbf{R}/w_{t}),$

for a diagonal dynamics covariance $\mathbf{Q}$, $a_{w},b_{w}>0$, and diagonal
observation covariance $\mathbf{R}$. They consider the class of variational
distributions

$q(\bar{w},\bar{{\bm{\theta}}})=q({\bm{\theta}}_{0})\prod_{t=1}^{T}q({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})q(w_{t}),$
(33)

with $\bar{{\bm{\theta}}}=({\bm{\theta}}_{0},\ldots,{\bm{\theta}}_{T})$ and
$\bar{w}=(w_{0},\ldots,w_{T})$. They show, for known $\mathbf{F}$,
$\mathbf{H}$, $\mathbf{Q}$, and $\mathbf{R}$, that the variational
distributions are of the form

$\displaystyle q({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$
$\displaystyle=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t}),$
(34) $\displaystyle q(w_{t})$
$\displaystyle=\text{Gam}(w_{t}|a_{w,t},b_{w,t}),$ (35)

where

$\displaystyle a_{w,t}$ $\displaystyle=a_{w}+\frac{1}{2},$ (36) $\displaystyle
b_{w,t}$ $\displaystyle=b_{w}+\mathbb{E}_{{\bm{\theta}}\sim{\cal
N}(\bm{\mu}_{t},\bm{\Sigma}_{t})}[({\bm{y}}_{t}-\mathbf{H}{\bm{\theta}})^{\intercal}\mathbf{R}^{-1}({\bm{y}}_{t}-\mathbf{H}{\bm{\theta}})],$
$\displaystyle\bm{\Sigma}_{t}^{-1}$
$\displaystyle=\mathbf{Q}^{-1}+v_{t}\mathbf{H}^{\intercal}\mathbf{R}^{-1}\mathbf{H},$
$\displaystyle\mathbf{K}_{t}$
$\displaystyle=\bm{\Sigma}_{t}\mathbf{H}^{\intercal}\mathbf{R}^{-1},$
$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\mathbf{F}\bm{\mu}_{t-1}v_{t}\,\mathbf{K}_{t}({\bm{x}}_{t}-\mathbf{H}\mathbf{F}\bm{\mu}_{t-1}),$
$\displaystyle v_{t}$
$\displaystyle=\frac{a_{w,t}+\frac{1}{2}}{b_{w,t}+\mathbb{E}_{{\bm{\theta}}\sim{\cal
N}(\bm{\mu}_{t},\bm{\Sigma}_{t})}[({\bm{y}}_{t}-\mathbf{H}{\bm{\theta}})^{\intercal}\mathbf{R}^{-1}({\bm{y}}_{t}-\mathbf{H}{\bm{\theta}})]}.$

Their method assumes no prior knowledge of either the measurement matrix
$\mathbf{H}$ or the projection matrix $\mathbf{F}$. These are estimated using
the EM algorithm. Assuming known $\mathbf{H}$ and $\mathbf{F}$ — as we do in
this paper — allows us to bypass the M-step. However, this is detrimental to
their approach since no information about the posterior covariance is
propagated forward.

###  A.5 Huang et al. (2016)

Huang et al. (2016) extend the state-space to be
$\bm{\Psi}_{t}=({\bm{\theta}}_{t},\mathbf{R}_{t},\nu_{t},w_{t})$, where
$\mathbf{R}_{t}$ is the measurement covariance, $w_{t}$ is a weighting term
for the measurement covariance, and $\nu_{t}$ are the degrees of freedom for
the weighting term. The SSM takes the form:

$\displaystyle p({\bm{\theta}}_{t}|{\bm{\theta}}_{t-1})$ $\displaystyle={\cal
N}({\bm{\theta}}_{t}\,|\,f_{t}({\bm{\theta}}_{t-1}),\,\mathbf{Q}_{t}),$ (37)
$\displaystyle p(\nu_{t})$ $\displaystyle={\rm Gam}(\nu_{t}\,|\,a_{t},b_{t}),$
$\displaystyle p(w_{t})$ $\displaystyle={\rm
Gam}(\lambda_{t}|\nu_{t}/2,\nu_{t}/2),$ $\displaystyle p(\mathbf{R}_{t})$
$\displaystyle={\cal W}^{-1}(\mathbf{R}_{t}\,|\,\bm{\Lambda}_{t},u_{t}),$
$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t},\mathbf{R}_{t})$
$\displaystyle={\cal
N}({\bm{y}}_{t}\,|\,h_{t}({\bm{\theta}}_{t}),\,\mathbf{R}_{t}/w_{t}),$

with $\mathbf{Q}_{t}=\text{diag}(q_{t,1},\ldots,q_{t,D})$. Note that their
method combines the SSMs in Ting et al. (2007) and Agamennoni et al. (2012).

##  Appendix B WoLF-MD as a MAP estimator

In this section we show how to derive WoLF-MD as a MAP estimator. This is an
alternative derivation that circumvents the use of generalised Bayes.

###  B.1 Deriving the Mahalanobis IMQ term

Consider the modified observation model of (32):

$\displaystyle p(w_{t})$ $\displaystyle={\rm Gam}(w_{t}|\alpha,\beta),$ (38)
$\displaystyle p({\bm{y}}_{t}|\hat{{\bm{y}}}_{t})$
$\displaystyle=\mathcal{N}({\bm{y}}_{t}|\hat{{\bm{y}}}_{t},w_{t}^{-1}\mathbf{R}_{t}),$
(39)

with $\hat{{\bm{y}}}_{t}=\mathbf{H}_{t}\bm{\mu}_{t|t-1}$ known at time $t$,
and $\alpha,\beta>0$. The posterior on $w_{t}$ is

$\displaystyle p(w_{t}|{\bm{y}}_{t})$ $\displaystyle\propto
w_{t}^{\alpha-1}e^{-\beta
w_{t}}\left|w_{t}\mathbf{R}_{t}^{-1}\right|^{1/2}\exp\left(-\frac{1}{2}{\bm{e}}_{t}^{\top}w_{t}\mathbf{R}_{t}^{-1}{\bm{e}}_{t}\right)$
(40) $\displaystyle\propto{\rm
Gam}\left(w_{t}|\alpha+\frac{n_{y}}{2},\beta+\frac{1}{2}\left\|{\bm{e}}_{t}\right\|_{\mathbf{R}_{t}^{-1}}^{2}\right),$
(41)

where ${\bm{e}}_{t}={\bm{y}}_{t}-\hat{{\bm{y}}}_{t}$ and
$\left\|{\bm{v}}\right\|_{\mathbf{A}}=\left\|\mathbf{A}^{1/2}{\bm{v}}\right\|_{2}$
is the Mahalanobis distance. The maximum-a-posteriori (MAP) estimate for
$w_{t}$ is

$W_{t}=\operatornamewithlimits{argmax}_{w_{t}\in{}^{+}}p(w_{t}|{\bm{y}}_{t})=\frac{\alpha+\frac{n_{y}}{2}-1}{\beta+\frac{1}{2}\left\|{\bm{e}}_{t}\right\|_{\mathbf{R}_{t}^{-1}}^{2}}.$
(42)

For a given $c\in\real$, take the hyperparameters $\alpha$ and $\beta$ to be

$\displaystyle\alpha=\frac{c^{2}-n_{y}+2}{2}\,,$
$\displaystyle\quad\beta=\frac{c^{2}}{2}\,,$ (43)

where $n_{y}$ is the number of measurements. We obtain

$\displaystyle w_{t}$
$\displaystyle=\left(1+\frac{\left\|{\bm{e}}_{t}\right\|_{\mathbf{R}_{t}^{-1}}^{2}}{c^{2}}\right)^{-1/2}.$
(44)

This is the Mahalanobis-based IMQ weighting function (18). Substituting the
MAP estimate back into the observation model yields the weighted loglikelihood
approximation

$\log p\left({\bm{y}}_{t}|\hat{{\bm{y}}}_{t}\right)\approx
w_{t}^{2}\log\mathcal{N}\left({\bm{y}}_{t}|\hat{{\bm{y}}}_{t},\mathbf{R}_{t}\right).$
(45)

###  B.2 Prior Uncertainty

In this section, we take the measurement mean to be the output of a predictive
model $\bar{{\bm{y}}}_{t}$ with unknown parameter ${\bm{\theta}}_{t}$. We let

$\displaystyle p\left({\bm{\theta}}_{t}|{\bm{y}}_{1:t-1}\right)$
$\displaystyle=\mathcal{N}\left({\bm{\theta}}_{t}|\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1}\right),$
(46) $\displaystyle p\left({\bm{y}}_{t}|{\bm{\theta}}_{t},w_{t}\right)$
$\displaystyle=\mathcal{N}\left({\bm{y}}_{t}|\bar{{\bm{y}}}_{t},w_{t}^{-1}\mathbf{R}_{t}\right),$
(47)

where $\bar{{\bm{y}}}_{t}$ is given by (11). The joint posterior is

$\displaystyle p\left({\bm{\theta}}_{t},w_{t}|{\bm{y}}_{1:t}\right)$
$\displaystyle\propto w_{t}^{\alpha-1}e^{-\beta
w_{t}}\exp\left(-\frac{1}{2}\left({\bm{\theta}}_{t}-\bm{\mu}_{t|t-1}\right)^{\top}\bm{\Sigma}_{t|t-1}^{-1}\left({\bm{\theta}}_{t}-\bm{\mu}_{t|t-1}\right)\right)$
$\displaystyle\qquad\times\left|w_{t}\mathbf{R}_{t}^{-1}\right|^{1/2}\exp\left(-\frac{1}{2}\left({\bm{y}}_{t}-\bar{{\bm{y}}}_{t}\right)^{\top}w_{t}\mathbf{R}_{t}^{-1}\left({\bm{y}}_{t}-\bar{{\bm{y}}}_{t}\right)\right)$
$\displaystyle\propto
w_{t}^{\alpha-1+n_{y}/2}\exp\left(-\frac{1}{2}\left\|{\bm{\theta}}_{t}-\bm{\mu}_{t|t-1}-\left(\bm{\Sigma}_{t|t-1}^{-1}+w_{t}\mathbf{H}_{t}^{\top}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)^{-1}w_{t}\mathbf{H}_{t}^{\top}\mathbf{R}_{t}^{-1}{\bm{e}}_{t}\right\|_{\bm{\Sigma}_{t|t-1}^{-1}+w_{t}\mathbf{H}_{t}^{\top}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}}^{2}\right)$
$\displaystyle\qquad\times\exp\left(-\beta
w_{t}-\frac{1}{2}{\bm{e}}_{t}^{\top}\left(w_{t}^{-1}\mathbf{R}_{t}+\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}\right)^{-1}{\bm{e}}_{t}\right).$

where the notation $||{{\bm{x}}}||_{\mathbf{A}}$ means
${\bm{x}}^{\intercal}\,\mathbf{A}^{-1}\,{\bm{x}}$. Then, the marginal for
$w_{t}$ can be written as

$\displaystyle p\left(w_{t}|{\bm{y}}_{1:t}\right)$ $\displaystyle\propto
w_{t}^{\alpha-1+n_{y}/2}\left|\bm{\Sigma}_{t|t-1}^{-1}+w_{t}\mathbf{H}_{t}^{\top}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right|^{-1/2}\exp\left(-\beta
w_{t}-\frac{1}{2}{\bm{e}}_{t}^{\top}\left(w_{t}^{-1}\mathbf{R}_{t}+\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}\right)^{-1}{\bm{e}}_{t}\right).$
(48)

By taking the limit $\bm{\Sigma}_{t|t-1}\to\bm{0}$, equation (48) becomes

$\lim_{\bm{\Sigma}_{t|t-1}\to\bm{0}}p\left(w_{t}|{\bm{y}}_{1:t}\right)\propto
w_{t}^{\alpha-1+n_{y}/2}\exp\left(-\beta
w_{t}-\frac{1}{2}{\bm{e}}_{t}^{\top}\mathbf{R}_{t}^{-1}{\bm{e}}_{t}w_{t}\right),$
(49)

with maximum at

$w_{t}^{*}=\operatornamewithlimits{argmax}_{w_{t}}\lim_{\bm{\Sigma}_{t|t-1}\to\bm{0}}=\frac{\alpha-1+\frac{n_{y}}{2}}{\beta+\frac{1}{2}{\bm{e}}_{t}^{\top}\mathbf{R}_{t}^{-1}{\bm{e}}_{t}},$
(50)

that matches (42). Therefore in the SSM setting, the Mahalanobis IMQ weighting
is the MAP estimate for $w_{t}$ after ignoring the prior uncertainty
$\bm{\Sigma}_{t|t-1}$.

##  Appendix C Proofs of theoretical results

###  C.1 Proof of Proposition 3.1

###### Proof.

Let $w_{t}^{2}:=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$. The loss function
takes the form

$\displaystyle\ell_{t}({\bm{\theta}}_{t})$ $\displaystyle=-w_{t}^{2}\log{\cal
N}({\bm{y}}_{t}\,|\,\mathbf{H}_{t}{\bm{\theta}}_{t},\,\mathbf{R}_{t})$ (51)
$\displaystyle=\frac{1}{2}\left({\bm{y}}_{t}-\mathbf{H}_{t}{\bm{\theta}}_{t}\right)^{\intercal}(\mathbf{R}_{t}/w_{t}^{2})^{-1}\left({\bm{y}}_{t}-\mathbf{H}_{t}{\bm{\theta}}_{t}\right)-\frac{w_{t}^{2}\,d}{2}\log\pi-\frac{w_{t}^{2}}{2}\log|\mathbf{R}_{t}|$
$\displaystyle=\frac{1}{2}\left({\bm{y}}_{t}-\mathbf{H}_{t}{\bm{\theta}}_{t}\right)^{\intercal}\bar{\mathbf{R}}_{t}^{-1}\left({\bm{y}}_{t}-\mathbf{H}_{t}{\bm{\theta}}_{t}\right)+C,$

with $\bar{\mathbf{R}}_{t}=\mathbf{R}_{t}/w_{t}^{2}$, and where
$C=-\frac{w_{t}^{2}\,d}{2}\log\pi-\frac{w_{t}^{2}}{2}\log|\mathbf{R}_{t}|$ is
a term that does not depend on ${\bm{\theta}}_{t}$. The remaining follows from
the standard KF derivation. Note that the loss function does not correspond to
the log-likelihood for a homoskedastic Gaussian model since
$\bar{\mathbf{R}}_{t}$ may depend on all data, including ${\bm{y}}_{t}$. ∎

Figure 8 shows the weighted log-likelihood (51) for a univariate ${\cal
N}(0,1)$ Gaussian density as a function of the weighting term
$w_{t}^{2}\in(0,1]$.

Figure 8: Weighted likelihood (unnormalised) for a standard Gaussian.

We observe that a weighting log-likelihood resembles a heavy-tailed likelihood
for $w_{t}<1$.

###  C.2 Proof of Theorem 3.2

Theorem 3.2 comprises several sub-claims, each of which we will prove as a
separate lemma in this section. At the end of the section, we will integrate
all the lemmas to prove the main result. The first results are for the linear
Gaussian SSM case.

######  Lemma C.1.

Consider the linear Gaussian SSM. The standard KF posterior has an unbounded
PIF and is not outlier robust.

###### Proof.

Let
$p({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$,
and
$p({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t}^{c},\bm{\Sigma}_{t}^{c})$,
the uncontaminated and contaminated standard Kalman filter posterior. Here

$\displaystyle\bm{\mu}_{t}=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\right),$
$\displaystyle\bm{\mu}_{t}^{c}=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}^{c}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}\right),$
$\displaystyle\bm{\Sigma}_{t}^{-1}=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t},$
$\displaystyle(\bm{\Sigma}_{t}^{c})^{-1}=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}.$

Here
$\bm{\mu}_{t|t-1},\mathbf{K}_{t},\mathbf{H}_{t},\bm{\Sigma}^{-1}_{t|t-1},\mathbf{R}_{t}^{-1}$
do not depend on the contamination ${\bm{y}}_{t}^{c}$. It is clear that
$\bm{\Sigma}_{t}^{-1}=(\bm{\Sigma}_{t}^{c})^{-1}$, hence
$\mathbf{K}_{t}^{c}=\bm{\Sigma}_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}=\bm{\Sigma}_{t}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}=\mathbf{K}_{t}$.
Using the fact that given two $d$-dimensional Gaussians
$\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$ and
$\mathcal{N}(\bm{\mu}_{1},\bm{\Sigma}_{1})$, the KL divergence is

$\text{KL}(\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})\|\mathcal{N}(\bm{\mu}_{1},\bm{\Sigma}_{1}))=\frac{1}{2}\left(\operatorname{Tr}\left(\bm{\Sigma}_{1}^{-1}\bm{\Sigma}_{0}\right)-m+\left(\bm{\mu}_{1}-\bm{\mu}_{0}\right)^{\intercal}\bm{\Sigma}_{1}^{-1}\left(\bm{\mu}_{1}-\bm{\mu}_{0}\right)+\ln\left(\frac{\det\bm{\Sigma}_{1}}{\det\bm{\Sigma}_{0}}\right)\right),$

and given that $\operatorname{Tr}(\mathbf{I}_{p})=p$, we can derive the PIF as

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})$
$\displaystyle=\operatorname{KL}\left(p({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})\|p({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})\right)$
$\displaystyle=\frac{1}{2}\Bigg{(}\operatorname{Tr}\left(\bm{\Sigma}_{t}^{-1}\bm{\Sigma}_{t}\right)-m+\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)^{\intercal}\bm{\Sigma}_{t}^{-1}\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)+\ln\left(\frac{\det\bm{\Sigma}_{t}}{\det\bm{\Sigma}_{t}}\right)\Bigg{)}$
$\displaystyle=\operatorname{Tr}\left(\mathbf{I}_{m}\right)-m+\frac{1}{2}\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)^{\intercal}\bm{\Sigma}_{t}^{-1}\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)$
$\displaystyle=\frac{1}{2}\left({\bm{y}}_{t}-{\bm{y}}_{t}^{c}\right)^{\intercal}\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}\left({\bm{y}}_{t}-{\bm{y}}_{t}^{c}\right).$

Here, $\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}$ is
positive definite. To see this, first note that for every invertible matrix
$B$ and conformable positive definite matrix $A$, $B^{\intercal}AB$ is
positive definite (Chapter 7.1 Horn & Johnson, 2012). Since $\mathbf{K}_{t}$
is a product of invertible matrices, it therefore is itself invertible. In
fact, let ${\bm{z}}\in\mathbb{R}^{m}$ be a non-zero vector. We can write
${\bm{z}}^{\intercal}B^{\intercal}AB{\bm{z}}=\tilde{{\bm{z}}}^{\intercal}A\tilde{{\bm{z}}}$,
for $\tilde{{\bm{z}}}=B{\bm{z}}$. We know that $\tilde{{\bm{z}}}=B{\bm{z}}\neq
0$, because $B$ is invertible. Then, since $A$ is positive definite and
$\tilde{{\bm{z}}}\in{}^{m}$, it holds that
$\tilde{{\bm{z}}}^{\intercal}A\tilde{{\bm{z}}}>0$, proving that
$B^{\intercal}AB$ is positive definite. Finally, since $\bm{\Sigma}_{t}^{-1}$
is positive definite, it therefore holds that
$\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}$ is positive
definite, too.

This in turn allows us to lower-bound the PIF by

$\displaystyle\frac{1}{2}\lambda_{\min}(\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t})\|{\bm{y}}_{t}-{\bm{y}}_{t}^{c}\|^{2}_{2}\leq\frac{1}{2}\left({\bm{y}}_{t}-{\bm{y}}_{t}^{c}\right)^{\intercal}\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}\left({\bm{y}}_{t}-{\bm{y}}_{t}^{c}\right)=\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t}),$

where we use the fact that for every positive definite matrix $A$ and vector
${\bm{z}}$, we have
$\lambda_{\min}(A)\|{\bm{z}}\|^{2}\leq{\bm{z}}^{\intercal}A{\bm{z}}\leq\lambda_{\max}(A)\|{\bm{z}}\|^{2}$,
where $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ are the minimum and maximum
eigenvalues of $A$, respectively (Chapter 5.7 Horn & Johnson, 2012). Moreover,
we know that
$\lambda_{\min}\left(\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}\right)>0$,
since $\mathbf{K}_{t}^{\intercal}\bm{\Sigma}_{t}^{-1}\mathbf{K}_{t}$ is
positive definite, and it does not depend on ${\bm{y}}_{t}^{c}$. Therefore, it
indeed holds that
$\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})\rightarrow+\infty$ as
$\|{\bm{y}}_{t}-{\bm{y}}_{t}^{c}\|^{2}_{2}\to+\infty$. ∎

######  Lemma C.2.

Consider the linear Gaussian SSM. The generalised posterior presented in
Proposition 3.1 has bounded PIF and is, therefore, outlier robust for any
weighting function $W$ such that
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{1:t})<\infty$ and
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{1:t})^{2}\,\|{\bm{y}}_{t}\|_{2}<\infty$.

###### Proof.

Let $W$ be a weighting function such that
$\sup_{{\bm{y}}_{t}\in{}^{d}}W({\bm{y}}_{1:t})<\infty$ and
$\sup_{{\bm{y}}_{t}\in{}^{d}}W({\bm{y}}_{1:t})^{2}\|{\bm{y}}_{t}\|_{2}<\infty$.
Lets define $w_{t}:=W^{2}({\bm{y}}_{1:t})$, and
$w_{t}^{c}:=W^{2}({\bm{y}}_{t}^{c},{\bm{y}}_{1:1-t})$. Let
$q({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$,
and
$q({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t}^{c},\bm{\Sigma}_{t}^{c})$,
where

$\displaystyle\bm{\mu}_{t}=\bm{\mu}_{t|t-1}+w_{t}\mathbf{K}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\right),$
$\displaystyle\bm{\mu}_{t}^{c}=\bm{\mu}_{t|t-1}+w_{t}^{c}\mathbf{K}_{t}^{c}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}\right),$
$\displaystyle\bm{\Sigma}_{t}^{-1}=\bm{\Sigma}_{t|t-1}^{-1}+w_{t}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t},$
$\displaystyle(\bm{\Sigma}_{t}^{c})^{-1}=\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t},$

as in Algorithm 1 Then, the PIF has the form

$\displaystyle\operatorname{PIF}({\bm{y}}_{t},{\bm{y}}_{1:t-1})=\frac{1}{2}\Bigg{(}$
$\displaystyle\underbrace{\operatorname{Tr}\left((\bm{\Sigma}_{t}^{c})^{-1}\bm{\Sigma}_{t}\right)-p}_{(1)}+\underbrace{\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)}_{(2)}+\underbrace{\ln\left(\frac{\det\bm{\Sigma}_{t}^{c}}{\det\bm{\Sigma}_{t}}\right)}_{(3)}\Bigg{)}.$

Now, we will get a bound for each term in the PIF. The first term can be
bounded as

$\displaystyle(1)$
$\displaystyle=\operatorname{Tr}\left((\bm{\Sigma}_{t}^{c})^{-1}\bm{\Sigma}_{t}\right)-m\leq\operatorname{Tr}\left((\bm{\Sigma}_{t}^{c})^{-1}\right)\operatorname{Tr}\left(\bm{\Sigma}_{t}\right)-m,$

where we use the fact that for two positive semidefinite matrices $A$, $B$, it
holds that
$\operatorname{Tr}(AB)\leq\operatorname{Tr}(A)\operatorname{Tr}(B)$. Observing
that $\bm{\Sigma}_{t}$ does not depend on $y_{t}^{c}$, we can now write
$C_{1}=\bm{\Sigma}_{t}$, so that by using the arithmetic rules of traces, we
obtain

$\displaystyle(1)$ $\displaystyle\leq
C_{1}\operatorname{Tr}\left(\left(\bm{\Sigma}_{t}^{c}\right)^{-1}\right)-m$
$\displaystyle=C_{1}\operatorname{Tr}\left(\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)-m$
$\displaystyle=C_{1}\operatorname{Tr}\left(\bm{\Sigma}_{t|t-1}^{-1}\right)+C_{1}w_{t}^{c}\operatorname{Tr}\left(\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)-m.$

Here, we use the fact that the trace of a sum is a sum of traces, and the
trace of the constant times matrix is equal to the constant times trace of the
matrix. Finally, since $\sup_{{\bm{y}}_{t}^{c}\in\mathbb{R}^{d}}w_{t}^{c}\leq
C_{2}<\infty$, the entire expression can be bounded by a constant
$C_{3}<\infty$ that does not depend on the contamination ${\bm{y}}_{t}^{c}$ as

$\displaystyle(1)$ $\displaystyle\leq
C_{1}\operatorname{Tr}\left(\bm{\Sigma}_{t|t-1}^{-1}\right)+C_{1}C_{2}\operatorname{Tr}\left(\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)-m=C_{3},$

where we use the fact that both traces are finite since both matrices are
real-valued. Next, we bound the second term by noting that it is the squared
Mahalanobis norm of $\bm{\mu}_{t}-\bm{\mu}_{t}^{c}$ with respect to
$\bm{\Sigma}_{t}^{c}$, and we write it as
$\|\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\|_{\bm{\Sigma}_{t}^{c}}$ (in particular, it
then satisfies all the properties of a norm). Therefore, we apply the triangle
inequality to obtain

$\displaystyle(2)$
$\displaystyle=(\|\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\|_{\bm{\Sigma}_{t}^{c}})^{2}$
$\displaystyle=(\|w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})-w_{t}^{c}\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}})^{2}$
$\displaystyle\leq(\|w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}+\|w_{t}^{c}\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}})^{2}$
$\displaystyle\leq
2(\|w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}^{2}+\|w_{t}^{c}\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}^{2}).$

In the last inequality, we use the fact that, for two real numbers $a$ and
$b$, it holds that $(a+b)^{2}\leq 2(a^{2}+b^{2})$. Then, we bound each term
separately,

$\displaystyle
2\|w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}^{2}$
$\displaystyle=2w_{t}^{2}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})^{\intercal}\mathbf{K}_{t}^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})$
$\displaystyle\leq
2w_{t}^{2}\lambda_{\max}\left(\mathbf{K}_{t}^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\mathbf{K}_{t}\right)\|{\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\|_{2}^{2},$

where
$2\lambda_{\max}\left(\mathbf{K}_{t}^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\mathbf{K}_{t}\right)\leq
C_{4}<\infty$, since
$\mathbf{K}_{t}^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\mathbf{K}_{t}$ is a
positive definite matrix, and thus, all eigenvalues are real-valued. This
property arises from the fact that for every invertible matrix $B$ and a
positive definite matrix $A$, $B^{\intercal}AB$ is positive definite. Finally,
we must verify that $(\bm{\Sigma}_{t}^{c})^{-1}$ is positive definite. Let
${\bm{z}}\in{}^{m}$ be a non-zero vector. Then,

$\displaystyle{\bm{z}}^{\intercal}\big{(}\bm{\Sigma}_{t}^{c}\big{)}^{-1}{\bm{z}}$
$\displaystyle={\bm{z}}^{\intercal}\big{(}\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\big{)}{\bm{z}}={\bm{z}}^{\intercal}\big{(}\bm{\Sigma}_{t|t-1}^{-1}\big{)}{\bm{z}}+w_{t}^{c}\big{(}{\bm{z}}^{\intercal}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}{\bm{z}}\big{)}.$

We know that $\bm{\Sigma}_{t|t-1}^{-1}$ is positive definite, hence
${\bm{z}}^{\intercal}\big{(}\bm{\Sigma}_{t|t-1}^{-1}\big{)}{\bm{z}}>0$.
Moreover,
${\bm{z}}^{\intercal}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}{\bm{z}}>0$
because $\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}$ is
positive definite. Since
$w_{t}^{c}=W^{2}({\bm{y}}_{t}^{c},{\bm{y}}_{1:1-t})\geq 0$ for all
${\bm{y}}_{t}^{c}\in\mathbb{R}^{d}$, then
$w_{t}^{c}\big{(}{\bm{z}}^{\intercal}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}{\bm{z}}\big{)}\geq
0$. Finally, combining these inequalities,
${\bm{z}}^{\intercal}\big{(}\bm{\Sigma}_{t}^{c}\big{)}^{-1}{\bm{z}}>0$.
Therefore, $(\bm{\Sigma}_{t}^{c})^{-1}$ is positive definite. Therefore,

$\displaystyle
2\|w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}^{2}\leq
C_{4}(w_{t}\|{\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\|_{2})^{2}=C_{5}$

since $(w_{t}\|{\bm{y}}_{t}-\hat{{\bm{y}}}_{t}\|_{2})^{2}$ does not depend on
the contamination. Similarly,

$\displaystyle
2\|w_{t}^{c}\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t})\|_{\bm{\Sigma}_{t}^{c}}^{2}\leq
C_{6}(w_{t}^{c}\|{\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}\|_{2})^{2},$

where, using the same argument as before,
$2\lambda_{\max}\left((\mathbf{K}_{t}^{c})^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\mathbf{K}_{t}^{c}\right)\leq
C_{6}<\infty$. Now, since
$\sup_{{\bm{y}}_{t}^{c}\in\mathbf{R}^{d}}w_{t}^{c}\|{\bm{y}}_{t}^{c}\|_{2}\leq
C_{7}<\infty$, we have:

$\displaystyle
C_{6}(w_{t}^{c}\|{\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}\|_{2})^{2}$
$\displaystyle\leq
C_{6}(w_{t}^{c}\|{\bm{y}}_{t}^{c}\|_{2}+w_{t}^{c}\|\hat{{\bm{y}}}_{t}\|_{2})^{2}\leq
C_{6}(C_{7}+C_{2}\|\hat{{\bm{y}}}_{t}\|_{2})^{2}=C_{8}.$

Putting it all together, we find that

$\displaystyle(2)\leq C_{5}+C_{8}.$

Lastly, the third and final term can be rewritten using properties of
determinants as

$\displaystyle(3)$
$\displaystyle=\ln\left(\frac{\det\bm{\Sigma}_{t}^{c}}{\det\bm{\Sigma}_{t}}\right)=\ln\left(\frac{1}{\det\bm{\Sigma}_{t}}\right)+\ln\left(\det\bm{\Sigma}_{t}^{c}\right)=\ln\left(\frac{1}{\det\bm{\Sigma}_{t}}\right)+\ln\left(\frac{1}{\det(\bm{\Sigma}_{t}^{c})^{-1}}\right).$

We define $C_{9}=\ln\left(\frac{1}{\det\bm{\Sigma}_{t}}\right)$ since it does
not depend on the contamination, and write

$\displaystyle(3)$
$\displaystyle=C_{9}+\ln\left(\frac{1}{\det(\bm{\Sigma}_{t}^{c})^{-1}}\right)$
$\displaystyle=C_{9}+\ln\left(\frac{1}{\det\left(\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)}\right)$
$\displaystyle\leq
C_{9}+\ln\left(\frac{1}{\det\left(\bm{\Sigma}_{t|t-1}^{-1}\right)+\det\left(w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)}\right),$

where in the last inequality, we use the fact that for two positive
semidefinite matrices $A$, $B$, it also holds that
$\det(A+B)\geq\det(A)+\det(B)$. Finally,

$\displaystyle(3)$ $\displaystyle\leq
C_{9}+\ln\left(\frac{1}{\det\left(\bm{\Sigma}_{t|t-1}^{-1}\right)+\det\left(w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)}\right)$
$\displaystyle\leq
C_{9}+\ln\left(\frac{1}{\det\left(\bm{\Sigma}_{t|t-1}^{-1}\right)}\right)=C_{10}.$

Here, we use the fact that
$w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}$ is
positive semidefinite, as we showed previously. Therefore
$\det\left(w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}\right)\geq
0$. By putting the bounds for (1), (2), and (3) together, we obtain

$\displaystyle\operatorname{PIF}({\bm{y}}_{t},{\bm{y}}_{1:t-1})$
$\displaystyle=\frac{1}{2}\Bigg{(}\underbrace{\operatorname{Tr}\left((\bm{\Sigma}_{t}^{c})^{-1}\bm{\Sigma}_{t}\right)-d}_{(1)}+\underbrace{\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)^{\intercal}(\bm{\Sigma}_{t}^{c})^{-1}\left(\bm{\mu}_{t}-\bm{\mu}_{t}^{c}\right)}_{(2)}+\underbrace{\ln\left(\frac{\det\bm{\Sigma}_{t}^{c}}{\det\bm{\Sigma}_{t}}\right)}_{(3)}\Bigg{)}.$
$\displaystyle\leq C_{3}+C_{5}+C_{8}+C_{10}<\infty.$

∎

We now extend the result to the linearised approximation of the SSM case.

######  Lemma C.3.

Consider the linearised approximation of the SSM. The standard EKF posterior
has an unbounded PIF and is not outlier robust.

###### Proof.

We can easily replicate the procedure in Lemma C.1 to the EKF since it can be
straightforwardly applied to the approximate posterior presented in Section
2.2. Specifically, for the standard EKF, we compute:

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})=\operatorname{KL}\left(p({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})\|p({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})\right)$
(52)

for
$p({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$,
and
$p({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t}^{c},\bm{\Sigma}_{t}^{c})$,
representing the uncontaminated and contaminated standard EKF posterior. In
particular,

$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}),$
$\displaystyle\bm{\mu}_{t}^{c}$
$\displaystyle=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}),$
$\displaystyle\bm{\Sigma}_{t}^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t},$
$\displaystyle(\bm{\Sigma}_{t}^{c})^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}.$

Here,
$\bm{\mu}_{t|t-1},\mathbf{K}_{t},\mathbf{H}_{t},\bm{\Sigma}^{-1}_{t|t-1},\mathbf{R}_{t}^{-1}$
do not depend on the contamination ${\bm{y}}_{t}^{c}$, and
$\bm{\mu}_{t|t-1}=\mathbb{E}[\bar{\bm{\mu}}_{t|t-1}|\bm{\mu}_{t-1}]=f_{t}(\bm{\mu}_{t-1})$,
$\hat{{\bm{y}}}_{t}=\mathbb{E}[\bar{{\bm{y}}}_{t}]=h_{t}(\bm{\mu}_{t|t-1})$,
where $\bar{\bm{\mu}}_{t|t-1}$ and $\bar{{\bm{y}}}_{t}$ are defined as
Equation 10 and Equation 11 respectively, and $\mathbf{H}_{t}$ is the Jacobian
of $h_{t}$ evaluated at $\bm{\mu}_{t|t-1}$. Since it follows the same
structure as the standard Kalman filter, replicating the procedure in Lemma
C.1, we obtain that the standard EKF is not robust. ∎

######  Lemma C.4.

Consider the linearised approximation of the SSM. The generalised posterior
presented in Proposition 3.1 has bounded PIF and is, therefore, outlier robust
for any weighting function $W$ such that
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})<\infty$
and
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\,\|{\bm{y}}_{t}\|_{2}<\infty$.

###### Proof.

We can easily replicate the procedure in Lemma C.2 since it can be
straightforwardly applied to the approximate posterior presented in Section
3.2. Now consider the weighted EKF. We compute

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})=\operatorname{KL}\left(q({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})\|q({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})\right).$
(53)

Let $W:{}^{d\times d}\to\real$ be a weighting function such that
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})<\infty$
and
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\|{\bm{y}}_{t}\|_{2}<\infty$.
Define $w_{t}:=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$, and
$w_{t}^{c}:=W^{2}({\bm{y}}_{t}^{c},\hat{{\bm{y}}}_{t})$.

Let
$q({\bm{\theta}}_{t}|{\bm{y}}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$,
and
$q({\bm{\theta}}_{t}|{\bm{y}}^{c}_{t},{\bm{y}}_{1:t-1})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t}^{c},\bm{\Sigma}_{t}^{c})$,
where

$\displaystyle\bm{\mu}_{t}$
$\displaystyle=\bm{\mu}_{t|t-1}+w_{t}\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}),$
$\displaystyle\bm{\mu}_{t}^{c}$
$\displaystyle=\bm{\mu}_{t|t-1}+w_{t}^{c}\mathbf{K}_{t}^{c}({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t}),$
$\displaystyle\bm{\Sigma}_{t}^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+w_{t}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t},$
$\displaystyle(\bm{\Sigma}_{t}^{c})^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+w_{t}^{c}\mathbf{H}_{t}^{\intercal}\mathbf{R}_{t}^{-1}\mathbf{H}_{t}.$

for $\bm{\mu}_{t|t-1}$ and $\hat{{\bm{y}}}_{t}$ as before. Then, following the
same procedure as in the proof of Lemma C.2, we show that this PIF is bounded
and the method is robust. ∎

##### Proof of Theorem 3.2

The proof of Theorem 3.2 follows directly from Lemmas C.1, C.2, C.3 and C.4.

###  C.3 Ensemble Kalman Filter

Now, we extend Theorem 3.2 to the ensemble Kalman filter case. There are a
couple of caveats to proving robustness for this case. While we know the
distribution of the particles in the case where the state-space model is
linear and Gaussian, this is not always the case. Therefore, we propose
studying the empirical measures defined by the particles. The problem with
this is that the Kullback-Leibler divergence is not defined for empirical
measures with different supports. This is why we propose using the
2-Wasserstein distance instead since it is unbounded and well-defined for
empirical measures. The 2-Wasserstein distance is defined as follows: If $P$
is an empirical measure with samples $x_{1},...,x_{n}$ and $Q$ is an empirical
measure with samples $y_{1},...,y_{n}$,

$\displaystyle
D_{W_{2}}(P,Q)=\inf_{\pi}\left(\frac{1}{n}\sum_{i=1}^{n}\left\|x_{i}-y_{\pi(i)}\right\|_{2}^{2}\right)^{1/2},$

where the infimum is over all permutations $\pi$ of $n$ elements. Therefore,
the PIF for the ensemble Kalman filter case has the form

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})=D_{W_{2}}\left(\mathbb{P}_{N},\mathbb{P}_{N}^{c}\right),$
(54)

where $\mathbb{P}_{N}$ is the empirical measure of the (non-contaminated)
particles
$\left\\{\hat{\theta}^{(i)}_{t}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\\}_{i=1}^{N}$,
and $\mathbb{P}_{N}^{c}$ is the empirical measure of the contaminated
particles
$\left\\{\hat{\theta}^{(i)}_{t^{c}}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\\}_{i=1}^{N}$.
As in the previous definition, if
$\sup_{{\bm{y}}_{t}^{c}\in\mathbb{R}^{d}}|\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})|<\infty$,
then the posterior is called outlier-robust.

######  Lemma C.5.

The generalised posterior presented in Section D.3 has bounded PIF (as defined
in Equation 54) and is, therefore, outlier robust for any weighting function
$W:{}^{d\times d}\to\real$ such that
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})<\infty$
and
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\,\|{\bm{y}}_{t}\|_{2}<\infty$.

###### Proof.

Let $W:{}^{d\times d}\to\real$ be a weighting function such that
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})<\infty$
and
$\sup_{{\bm{y}}_{t}\in\mathbf{R}^{d}}W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\|{\bm{y}}_{t}\|_{2}<\infty$.
Define $w_{t}:=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$, and
$w_{t}^{c}:=W^{2}({\bm{y}}_{t}^{c},\hat{{\bm{y}}}_{t})$. Then, the PIF

$\displaystyle\operatorname{PIF}({\bm{y}}_{t}^{c},{\bm{y}}_{1:t})=D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})=\inf_{\pi}\left(\frac{1}{N}\sum_{i=1}^{N}\left\|\hat{\theta}^{(i)}_{t}-\hat{\theta}^{(\pi(i))}_{t^{c}}\right\|_{2}^{2}\right)^{1/2},$

where the infimum is over all permutations $\pi$ of $N$ elements,
$\mathbb{P}_{N}$ is the empirical measure of the particles
$\left\\{\hat{\theta}^{(i)}_{t}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\\}_{i=1}^{N}$,
and $\mathbb{P}_{N}^{c}$ is the empirical measure of the contaminated
particles
$\left\\{\hat{\theta}^{(i)}_{t^{c}}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\\}_{i=1}^{N}$.
Therefore, we know that the infimum will be smaller than considering only the
identity as permutation,

$\displaystyle
D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})\leq\left(\frac{1}{N}\sum_{i=1}^{N}\left\|\hat{\theta}^{(i)}_{t}-\hat{\theta}^{(i)}_{t^{c}}\right\|_{2}^{2}\right)^{1/2}.$

Now, using the equation of the particles:

$\displaystyle D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})$
$\displaystyle\leq\left(\frac{1}{N}\sum_{i=1}^{N}\left\|\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)-\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}-w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\right)^{1/2}$
$\displaystyle=\left(\frac{1}{N}\sum_{i=1}^{N}\left\|w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)-w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\right)^{1/2}.$

We use triangle inequality to obtain:

$\displaystyle D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})$
$\displaystyle\leq\left(\frac{1}{N}\sum_{i=1}^{N}\left(\left\|w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}+\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}\right)^{2}\right)^{1/2}$
$\displaystyle\leq\left(\frac{2}{N}\sum_{i=1}^{N}\left\|w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}+\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\right)^{1/2}.$

Lets define
$C_{11}=\sum_{i=1}^{N}\|w_{t}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\|_{2}^{2}$,
a constant that does not depend on the contamination. Therefore:

$\displaystyle D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})$
$\displaystyle\leq\left(\frac{2}{N}\left(C_{11}+\sum_{i=1}^{N}\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\right)\right)^{1/2}.$

First, we observe that
$\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}=\left(w_{t}^{c}\right)^{2}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)^{\intercal}\bar{\mathbf{K}}_{t}^{\intercal}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)$.
Now we bound this expression using the fact that for every positive definite
matrix $A$ and vector ${\bm{z}}$, we have
$\lambda_{\min}(A)\|{\bm{z}}\|^{2}\leq{\bm{z}}^{\intercal}A{\bm{z}}\leq\lambda_{\max}(A)\|{\bm{z}}\|^{2}$,
where $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ are the minimum and maximum
eigenvalues of $A$, respectively (Chapter 5.7 Horn & Johnson, 2012):

$\displaystyle\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\leq(w_{t}^{c})^{2}\lambda_{\max}(\bar{\mathbf{K}}_{t}^{\intercal}\bar{\mathbf{K}}_{t})\left\|{\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right\|_{2}^{2},$

where
$\lambda_{\max}\left(\bar{\mathbf{K}}_{t}^{\intercal}\bar{\mathbf{K}}_{t}\right)\leq
C_{12}<\infty$, since $\bar{\mathbf{K}}_{t}^{\intercal}\bar{\mathbf{K}}_{t}$
is a positive definite matrix, and thus, all eigenvalues are real-valued.
Then,

$\displaystyle\left\|w_{t}^{c}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right)\right\|_{2}^{2}\leq
C_{12}\left(w_{t}^{c}\left\|{\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right\|_{2}\right)^{2}.$

Now, since $\sup_{{\bm{y}}_{t}^{c}\in\mathbb{R}^{d}}w_{t}^{c}\leq
C_{13}<\infty$, and
$\sup_{{\bm{y}}_{t}^{c}\in\mathbf{R}^{d}}w_{t}^{c}\|{\bm{y}}_{t}^{c}\|_{2}\leq
C_{14}<\infty$, we have:

$\displaystyle
C_{12}\left(w_{t}^{c}\left\|{\bm{y}}_{t}^{c}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right\|_{2}\right)^{2}$
$\displaystyle\leq
C_{12}\left(w_{t}^{c}\left\|{\bm{y}}_{t}^{c}\right\|_{2}+w_{t}^{c}\left\|\hat{{\bm{y}}}_{t|t-1}^{(i)}\right\|_{2}\right)^{2}\leq
C_{12}\left(C_{14}+C_{13}\left\|\hat{{\bm{y}}}_{t|t-1}^{(i)}\right\|_{2}\right)^{2}=C_{15}.$

Putting it all together, we find that

$\displaystyle D_{W_{2}}(\mathbb{P}_{N},\mathbb{P}_{N}^{c})$
$\displaystyle\leq\left(\frac{2}{N}\left(C_{11}+\sum_{i=1}^{N}C_{15}\right)\right)^{1/2}<\infty.$

∎

##  Appendix D Extensions of the weighted likelihood filter

Below, we discuss a number of extension to the WoLF methodology including
generalisations to (i) exponential-family members, (ii) multi-output weighting
functions, and (iii) the EnKF.

###  D.1 Exponential family likelihoods

We extend WoLF for measurements modelled using an element of the exponential
family of distributions (as first mentioned in Section 3.2). Classical
examples of exponential families, in addition to the Gaussian distribution,
are the Bernoulli distribution, the Gamma distribution, the Beta distribution,
and the Poisson distribution. These distributions can be considered to tackle
filtering problems when the measurements are generated respectively from a
binary process, a process that only takes values in positive real line, a
process that takes values in the interval $[0,1]$, or a counting process.

We take the mass function of a measurement ${\bm{y}}_{t}\in B\subseteq{}^{d}$
to be of the form

$\displaystyle p({\bm{y}}_{t}|{\bm{\theta}}_{t})$
$\displaystyle=\text{expfam}({\bm{y}}_{t}|\bm{\eta}_{t})=Z^{-1}(\bm{\eta}_{t})\exp\Big{(}\bm{\eta}_{t}^{\intercal}T({\bm{y}}_{t})+b({\bm{y}}_{t})\Big{)}.$
(55)

with $B$ the support of the measurement ${\bm{y}}_{t}$,
$\bm{\eta}_{t}\in{}^{k}$ the natural parameters, $A:B\to{}^{k}$ the sufficient
statistic function, $b:B\to\real$ the base measure, and $Z:{}^{k}\to\real$ the
normalising function. To simplify the notation, let
$\bm{\eta}_{t}=h_{t}({\bm{\theta}}_{t})$ and ${\bm{z}}_{t}=T({\bm{y}}_{t})$.
We also define the dual (moment) parameters

$\bm{\lambda}_{t}=\mathbb{E}\left[{{\bm{z}}_{t}|\bm{\eta}_{t}}\right]=\nabla_{\bm{\eta}_{t}}\log
Z(\bm{\eta}_{t})$ (56)

and the conditional variance

$\mathbf{R}_{t}=\text{Cov}[{\bm{z}}_{t}|\bm{\eta}_{t}]=\nabla_{\bm{\eta}_{t}}\bm{\lambda}_{t}=\nabla_{\bm{\eta}_{t}}^{2}\log
Z(\bm{\eta}_{t}),$ (57)

and we take our predictive model to output the dual parameters:

$\displaystyle\bm{\lambda}_{t}=h_{t}({\bm{\theta}}_{t}).$ (58)

For example, for a Gaussian likelihood with fixed observation noise and a
linear observation model as in (5), we have that
$A({\bm{y}}_{t})={\bm{y}}_{t}$,
$\bm{\lambda}_{t}=\mathbf{H}_{t}{\bm{\theta}}_{t}$, and $\mathbf{R}_{t}$ is
constant.

The exponential family EKF algorithm of Ollivier (2018) approximates the
likelihood in (55) with a moment-matched Gaussian,

$\displaystyle
q({\bm{y}}_{t}|{\bm{\theta}}_{t})=\mathcal{N}({\bm{z}}_{t}|\bar{\bm{\lambda}}_{t},\mathbf{R}_{t}),$
(59)

where
$\bar{\bm{\lambda}}_{t}=\mathbf{H}_{t}({\bm{\theta}}_{t}-\bar{\bm{\mu}}_{t|t-1})+h_{t}(\bar{\bm{\mu}}_{t|t-1})$,
and $\bm{\mu}_{t|t-1}$ is given by (10). We modify this algorithm by including
a weighting term on the log-likelihood to get

$\displaystyle\log q({\bm{y}}_{t}|{\bm{\theta}}_{t})$
$\displaystyle=W_{t}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})^{2}\left(-\frac{1}{2}\bar{\bm{\lambda}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}\bar{\bm{\lambda}}_{t}+\bar{\bm{\lambda}}_{t}^{\top}{\mathbf{R}}_{t}^{-1}T({\bm{y}}_{t})\right)+{\rm
cst},$ (60)

with ${\rm cst}.$ a constant term that does not depend on ${\bm{\theta}}$. We
leave the study of $W({\bm{y}}_{t},\hat{{\bm{y}}}_{t})$ when modelling a non-
Gaussian exponential family for future work.

###  D.2 Dimension-specific weighting

In this section, we provide details about the implementation of a vector-
valued weighting function $W_{t}:{}^{d}\times{}^{d}\to{}^{d}$ introduced in
Section 3.3. See Section 4.3 for an evaluation of this method. We consider
${\bm{y}}_{t}\in{}^{d}$ with $d>1$.

If the likelihood factorises as
$p\left({\bm{y}}_{t}|{\bm{\theta}}_{t}\right)=\prod_{j=1}^{d}p\left(y_{t,j}|{\bm{\theta}}_{t}\right)$,
with $y_{t,j}$ the $j$-th element of the $t$-th measurement, we define the
weighted log-likelihood as

$\log
q\left({\bm{y}}_{t}|{\bm{\theta}}_{t}\right)=\sum_{j=1}^{d}w_{t,j}^{2}\log
p\left(y_{t,j}|{\bm{\theta}}_{t}\right),$ (61)

where $w_{t,j}^{2}=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})_{j}$ is the $j$-th
entry of the vector-valued weight function.

For a Gaussian likelihood
$p({\bm{y}}_{t}|{\bm{\theta}}_{t})=\mathcal{N}\left({\bm{y}}_{t}|h_{t}\left({\bm{\theta}}_{t}\right),\mathbf{R}_{t}\right)$
we define the weighted likelihood as

$\displaystyle q({\bm{y}}_{t}|{\bm{\theta}}_{t})$
$\displaystyle=\mathcal{N}\left({\bm{y}}_{t}|h_{t}\left({\bm{\theta}}_{t}\right),\bar{\mathbf{R}}_{t}\right),$
(62) $\displaystyle\bar{\mathbf{R}}_{t}^{-1}$
$\displaystyle=\mathrm{Diag}\left({\bm{w}}_{t}\right)\mathbf{R}_{t}^{-1}\mathrm{Diag}\left({\bm{w}}_{t}\right),$
(63)

which is a special case of (61) when $\mathbf{R}_{t}$ is
diagonal.888$\mathrm{Diag}({\bm{v}})$ for ${\bm{v}}\in{}^{K}$ is the $K\times
K$ matrix with $\mathrm{Diag}({\bm{v}})_{i,i}={\bm{v}}_{i}$ and
$\mathrm{Diag}({\bm{v}})_{i,j}=0$ for $i\neq j$. This expression scales the
precision of each $y_{t,j}$ by $w_{t,j}^{2}$ while preserving correlations. It
can also be written in terms of dimension-specific weighting on the errors:

$\displaystyle\log q({\bm{y}}_{t}|{\bm{\theta}}_{t})$
$\displaystyle=-\frac{1}{2}\tilde{E}_{t}\left({\bm{\theta}}_{t}\right)^{\intercal}\mathbf{R}_{t}^{-1}\tilde{E}_{t}\left({\bm{\theta}}_{t}\right),$
(64) $\displaystyle\tilde{E}_{t}\left({\bm{\theta}}_{t}\right)$
$\displaystyle=\mathrm{Diag}\left({\bm{w}}_{t}\right)\left({\bm{y}}_{t}-h_{t}\left({\bm{\theta}}_{t}\right)\right).$
(65)

The EKF update from the weighted loglikelihood in (62) mirrors the standard
EKF update with the true observation precision $\mathbf{R}_{t}^{-1}$ replaced
by $\bar{\mathbf{R}}_{t}^{-1}$:

$\displaystyle\bm{\Sigma}_{t}^{-1}$
$\displaystyle=\bm{\Sigma}_{t|t-1}^{-1}+\mathbf{H}_{t}^{\top}\bar{\mathbf{R}}_{t}^{-1}\mathbf{H}_{t},$
(66) $\displaystyle\bm{\mu}_{t}$
$\displaystyle=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}({\bm{y}}_{t}-\hat{{\bm{y}}}_{t}),$
(67) $\displaystyle\mathbf{K}_{t}$
$\displaystyle=\bm{\Sigma}_{t}\mathbf{H}_{t}^{\top}\bar{\mathbf{R}}_{t}^{-1}.$
(68)

This generalises the WoLF update with scalar weight given by Proposition 3.1.

Finally, when $w_{t,j}=0$ for one or more observation dimensions, (61) and
(62) define improper observation distributions because the marginals on these
dimensions are uniform over and the precision $\bar{\mathbf{R}}_{t}^{-1}$ is
singular. Nevertheless, the joint marginal for the dimensions with positive
weights, ${\bm{y}}_{t,J_{t}}$ where $J_{t}=\\{j:w_{t,j}>0\\}$, is a proper
distribution. This is all that is needed because the observations $y_{t,j}$
with $w_{t,j}=0$ are ignored in the update. This can be seen in (66), (67),
and (68), where as a consequence of (63), the error terms and the Jacobian are
both zeroed out on the zero-weighted dimensions.

###  D.3 Weighted ensembles

In this section, we introduce the weighted EnKF first mentioned in Section
3.2. We begin with a discussion of the EnKF originally referenced in Section
2.3.

##### The ensemble Kalman filter:

The ensemble Kalman filter (EnKF) (Roth et al., 2017b) uses an ensemble of
$N\in\mathbb{N}$ particles
$\\{\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}\\}_{i=1}^{N}$. For each
$i=1,\ldots,N$, the update step samples predictions
$\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}$ according to

$\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}\sim\mathcal{N}\left(h_{t}\left(\hat{{\bm{\theta}}}_{t|t-1}^{(i)}\right),\mathbf{R}_{t}\right).$
(69)

and then updates each particle according to

$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right),$
(70)

The EnKF gain $\bar{\mathbf{K}}_{t}$ is

$\bar{\mathbf{K}}_{t}={\rm
cov}_{i}\left[\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)},\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}\right]\left(\mathbb{V}_{{i}}\left[{\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}}\right]\right)^{-1}.$
(71)

We write $\mathbb{E}_{{i}}\left[{\cdot}\right],{\rm
cov}_{i}\left[\cdot,\cdot\right],\mathbb{V}_{{i}}\left[{\cdot}\right]$ to
refer to the distribution over particles as indexed by $i$.

The EnKF converges to standard KF as $N\to\infty$ when the likelihood is
linear-Gaussian and the prior is also Gaussian,

$\displaystyle p\left({\bm{y}}_{t}|{\bm{\theta}}_{t}\right)$
$\displaystyle=\mathcal{N}\left({\bm{y}}_{t}|\mathbf{H}_{t}{\bm{\theta}}_{t},\mathbf{R}_{t}\right),$
(72) $\displaystyle\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}$
$\displaystyle\sim\mathcal{N}\left(\bm{\mu}_{t|t-1},\bm{\Sigma}_{t|t-1}\right).$
(73)

Under these conditions the EnKF gain $\bar{\mathbf{K}}_{t}$ matches the KF
gain $\mathbf{K}_{t}$, because

$\displaystyle{\rm
cov}_{i}\left[\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)},\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}\right]$
$\displaystyle=\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top},$ (74)
$\displaystyle\mathbb{V}_{{i}}\left[{\hat{{\bm{y}}}_{t|t-1}^{\left(i\right)}}\right]$
$\displaystyle=\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}+\mathbf{R}_{t},$
(75) $\displaystyle\bar{\mathbf{K}}_{t}$
$\displaystyle=\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}\left(\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}+\mathbf{R}_{t}\right)^{-1}$
(76) $\displaystyle=\mathbf{K}_{t},$ (77)

and therefore the statistics of the posterior ensemble also match the KF
update

$\displaystyle\bm{\mu}_{t}=\mathbb{E}_{{i}}\left[{\hat{{\bm{\theta}}}_{t}^{\left(i\right)}}\right]$
$\displaystyle=\bm{\mu}_{t|t-1}+\mathbf{K}_{t}\left({\bm{y}}_{t}-\mathbf{H}_{t}\bm{\mu}_{t|t-1}\right),$
(78)
$\displaystyle\bm{\Sigma}_{t}=\mathbb{V}_{{i}}\left[{\hat{{\bm{\theta}}}_{t}^{\left(i\right)}}\right]$
$\displaystyle=\bm{\Sigma}_{t|t-1}-\mathbf{K}_{t}\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}.$
(79)

##### The weighted-likelihood EnKF:

We propose a weighted version of EnKF based on our WoLF with dimension-
specific weights and Gaussian likelihood, where (69) is modified to sample
particle predictions from the weighted likelihood in (62):

$\displaystyle\hat{{\bm{y}}}_{t|t-1}^{(i)}\sim\mathcal{N}\left(h_{t}\left(\hat{{\bm{\theta}}}_{t|t-1}^{(i)}\right),\bar{\mathbf{R}}_{t}\right).$
(80)

The weighted EnKF converges to the WoLF as $N\to\infty$, by the same argument
for the vanilla EnKF and EKF in (74) through (79). Specifically, under this
limit the sampling scheme in (80) and the update in (70) yield

$\displaystyle\bar{\mathbf{K}}_{t}$
$\displaystyle=\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}\left(\mathbf{H}_{t}\bm{\Sigma}_{t|t-1}\mathbf{H}_{t}^{\top}+\bar{\mathbf{R}}_{t}\right)^{-1},$
(81) $\displaystyle\bm{\mu}_{t}$
$\displaystyle=\bm{\mu}_{t|t-1}+\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\mathbf{H}_{t}\bm{\mu}_{t|t-1}\right),$
(82) $\displaystyle\bm{\Sigma}_{t}$
$\displaystyle=\bm{\Sigma}_{t|t-1}-\bar{\mathbf{K}}_{t}\mathbf{H}_{t}\bm{\Sigma}_{t|t-1},$
(83)

which matches the WoLF update in (66), (67), and (68).

For the case of scalar weights with
$w_{t}^{2}=W^{2}({\bm{y}}_{t},\hat{{\bm{y}}}_{t})\in\\{0,1\\}$ as in (19), the
weighted EnKF can be implemented by unweighted sampling as in (69) followed by
a weighted update that replaces (70) with

$\hat{{\bm{\theta}}}^{(i)}_{t}=\hat{{\bm{\theta}}}^{(i)}_{t|t-1}+w_{t}^{2}\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right).$
(84)

That is, the update is simply skipped when $w_{t}^{2}=0$, in agreement with
Proposition 3.1 and Algorithm 1.

For the case of vector weights ${\bm{w}}_{t}$ as in Section D.2, when
${\bm{w}}_{t}\in\\{0,1\\}^{d}$ the weighted EnKF can be implemented by a
generalisation of (84) that ignores only those observation components
${\bm{y}}_{t,j}$ for which ${\bm{w}}_{t,j}=0$. Specifically, we can use the
vanilla EnKF sampling in (69) but sample only the positively weighted
dimensions, i.e. $\hat{{\bm{y}}}^{(i)}_{t|t-1,J_{t}}$ where
$J_{t}=\\{j:{\bm{w}}_{t,j}>1\\}$. Thus $\hat{{\bm{y}}}^{(i)}_{t|t-1}$ has size
$|J_{t}|$ and (71) yields $\bar{\mathbf{K}}_{t}$ with size $m\times|J_{t}|$,
and the particles can be updated according to

$\hat{{\bm{\theta}}}^{(i)}_{t}=\hat{{\bm{\theta}}}^{(i)}_{t|t-1}+\bar{\mathbf{K}}_{t}\left({\bm{y}}_{t,J_{t}}-\hat{{\bm{y}}}_{t|t-1,J_{t}}^{(i)}\right).$
(85)

This agrees with the vanilla EnKF update except that zero-weighted observation
dimensions are ignored. It can be shown to converge to the dimensionally
weighted WoLF described in Section D.2 as $N\to\infty$.999The general case of
${\bm{w}}_{t}\in[0,1]^{d}$ can be implemented by combining this scheme with
weighted sampling of positively weighted dimensions, i.e. using (80) to sample
$\hat{{\bm{y}}}^{(i)}_{t|t-1,J_{t}}$.

For the experiments reported here (Section 4.3), we take the shortcut of
sampling predictions for all dimensions and generalising (84) to vector
weights:

$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+\bar{\mathbf{K}}_{t}\mathrm{Diag}({\bm{w}}_{t})\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right).$
(86)

This approximates the method of sampling only
$\hat{{\bm{y}}}^{(i)}_{t|t-1,J_{t}}$ and updating with (85). This is more
efficient in our Jax (Bradbury et al., 2018) implementation because all arrays
are of constant size. We generalise (19) to define the weight vector as

${\bm{w}}_{t,j}=\begin{cases}1&\frac{1}{N}\sum_{i=1}^{N}\left({\bm{y}}_{t,j}-\hat{{\bm{y}}}_{t|t-1,j}^{(i)}\right)^{2}\leq
c,\\\ 0&\text{otherwise},\end{cases}$ (87)

We call this the average-particle EnKF (AP-EnKF).

Alternatively, we propose the update equation

$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+\bar{\mathbf{K}}_{t}\,\mathrm{Diag}\left({\bm{w}}_{t}^{(i)}\right)\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right),$
(88)

with

${\bm{w}}_{t,j}^{(i)}=\begin{cases}1&\left({\bm{y}}_{t,j}-\hat{{\bm{y}}}_{t|t-1,j}^{(i)}\right)^{2}\leq
c,\\\ 0&\text{otherwise},\end{cases}$ (89)

which we call the per-particle EnKF (PP-EnKF).

##### The Huberised EnKF:

The update equations for the PP-EnKF and AP-EnKF are related to the Huberised
EnKF (H-EnKF) algorithm of Roh et al. (2013) that update particles according
to

$\hat{{\bm{\theta}}}_{t}^{\left(i\right)}=\hat{{\bm{\theta}}}_{t|t-1}^{\left(i\right)}+\bar{\mathbf{K}}_{t}\,\rho\left({\bm{y}}_{t}-\hat{{\bm{y}}}_{t|t-1}^{(i)}\right),$
(90)

with

$\rho(z)=\begin{cases}c&\text{if }z>c,\\\ -c&\text{if }z<-c,\\\
z&\text{otherwise},\end{cases}$ (91)

$c>0$ a given threshold, and the function $\rho$ is applied element-wise.

##  Appendix E Additional numerical experriments

This section provides additional numerical experiments. In particular, in
Appendix E.1, we discuss Figure 2 in more detail; in Appendix E.2, we conduct
robustness checks for the numerical experiments presented in Section 4; and in
Appendix E.3, we provide an additional experiment for online learning in a
non-stationary environment with outlier measurements.

###  E.1 Description of the PIF

In this section, we discuss the 2d tracking problem of Section 4.1 in more
detail. We generate $t=20$ steps from (20) and the last measurement
${\bm{y}}_{t}$ is replaced with ${\bm{y}}_{t}^{c}={\bm{y}}_{t}+\bm{\epsilon}$,
where $\bm{\epsilon}\in[-5,5]\times[-5,5]$ is the outlier noise. Figure 9
(left) shows the weighting term for the WoLF methods as a function of the
outlier noise $\bm{\epsilon}$. Alternatively, Figure 9 (right) shows the
Mahalanobis distance between the prior predictive $\hat{{\bm{y}}}_{t}$ and the
outlier measurement ${\bm{y}}_{t}^{c}$.

Figure 9:  Weighting term $W_{t}({\bm{y}}_{1:t})$ given an the outlier
measurement ${\bm{y}}_{t}^{c}$. The $x$-axis and $y$-axis describe the errors
$\bm{\epsilon}\in[-5,5]\times[-5,5]$.

We observe that the midpoint of the contours for all methods is not centred at
$\bm{\epsilon}=(0,0)$. This is because the weighting term is a function of the
prior predictive and the measurement at time $t$. Hence, for the IMQ, the
weighting is $1$ only if the prior predictive equals the measurement, i.e.,
$\hat{{\bm{y}}}_{t}=\mathbf{H}_{t}\bm{\mu}_{t|t-1}$. This result explains the
distored PIF observed in Figure 2.

Figure 10 shows the PIF for the WoLF-IMQ and the WoLF-TMD first shown in
Figure 2.

Figure 10:  PIF of Figure 2 for extended domain of the errors
$\bm{\epsilon}\in[-15,10]\times[-15,10]$.

We observe that the PIF for the WoLF-TMD has a hard cutoff that strongly
bounds the PIF values at the expense of higher PIF around a region near the
cutoff boundary. In contrast, the WoLF-IMQ bounds the PIF values more softly
and has lower maximum PIF than the WoLF-TMD.

###  E.2 Robustness

####  E.2.1 2d tracking

In this section, we present further results from experiment 4.1. Figure 11
shows the distribution over errors over all 4 state components for both
outlier variants of the 2d tracking problem. For the mixture variant, we see
elongated tails in the error distribution for the  WoLF-TMD algorithm.

Figure 11:  Distribution of $J_{T,i}$ for all state components. The left panel
is for the student variant (21) and the right panel is for the mixture variant
(22).

To see why, consider the trace of weights over time in Figure 12. We observe
that both  WoLF-IMQ and  WoLF-TMD set the weighting term close to (or equal)
to zero in outlier events. However, the  WoLF-TMD is more prone to have false
positives (at a rate of about 7% in this example), in which it fails to update
the posterior state. This explains the elongated tails.

Figure 12: Trace of the weighting term for WoLF-IMQ and WoLF-MD.

##### Right panel of Figure 3

Figure 13: Sample of the top-left panel in Figure 3

##### 2d tracking: hyperparameter choice

Here we show the performance of each method for the experiment in Section 4.1
as we vary the hyperparameters. We do this for both, the Mixture case and the
Student-t case.

Figure 14 shows the RMSE for estimating the first state component as we vary
the $c$ hyperparameter for the IMQ and the TDM weighting functions; for  KF-
IW, we use two inner iterations and we vary the $h$ hyperparameter; for  KF-B,
we fix $\beta=1.0$, use four inner iterations and vary the value of $\alpha$.
The definition for the RMSE $J_{T,0}$ is as in Section 4.1.

Figure 14: $J_{t,0}$ over 100 runs for the methods as a function of their
hyperparameter.

##### Empirical comparison to (Boustati et al., 2020)

In this section, we compare the particle-filter-based method of (Boustati et
al., 2020), which we denote the RBPF. Figure 15 extends Figure 4 to include
the RBPF.

Figure 15:  Distribution (across 500 2d tracking trials) of RMSE for first
component of the state vector, $J_{T,0}$. Left panel: Student observation
model. Right panel: Mixture observation model. We extend Figure 4 to include
the RBPF.

The RBPF performs comparably to our proposed method; however, it has much
higher computational cost and does not have a closed-form solution.

####  E.2.2 UCI datasets

In this section, we give more details on the MLP regression experiments in
Section 4.2. The size of the datasets and models which we use are summarised
in Table 3.

| #Examples $T$ | #Features $d$ | #Parameters $m$
---|---|---|---
Dataset |  |  |
Boston | $506$ | $14$ | $321$
Concrete | $1,030$ | $9$ | $221$
Energy | $768$ | $9$ | $221$
Kin8nm | $8,192$ | $9$ | $221$
Naval | $11,934$ | $18$ | $401$
Power | $9,568$ | $5$ | $141$
Protein | $45,730$ | $10$ | $241$
Wine | $1,599$ | $12$ | $281$
Yacht | $308$ | $7$ | $181$
Table 3: Description of UCI datasets. Number of parameters refers to the size
of the one-layer MLP.

##### Speed vs accuracy

Figure 16 (left) shows the time-step and RMedSe for the Kin8nm dataset over
multiple trials. Our methods, the  WoLF-IMQ and the  WoLF-TMD, strike the best
balance between RMedSE and running time among the competing methods. The  EKF-
IW and the  EKF-B have comparable error rates to the WoLF methods, but their
running time is 6x and 12x slower than the  WoLF-IMQ, respectively. The  EKF
has comparable running time, but 3x higher average RMedSE than the  WoLF-IMQ.
Finally, the  OGD is 40% faster, but has 2x higher average RMedSE than the
WoLF-IMQ.

##### Sensitivity to outlier rate

We evaluate the sensitivity of the methods to the choice of corrosion rate
$p_{\epsilon}$ for the Kin8nm dataset. Figure 16 (right) shows the RMedSE
after 100 trials of each method as a function of $p_{\epsilon}$. We set the
hyperparameters using BO with $p_{\epsilon}=0.1$ and evaluate the RMedSE with
$p_{\epsilon}\in\\{0.00,0.05,\ldots,0.45\\}$. We observe that all robust
methods — the  WoLF-IMQ, the  WoLF-TMD, the  EKF-B, and the  EKF-IW— have
similar RMEdSE for any choice of $p_{\epsilon}$ and have similar RMedSE rate
of increase. Conversely, the  EKF and the  OGD are much less stable and their
RMedSE error rate increases at a faster pace.

Figure 16:  Left panel shows the RMedSE (root median squared error) vs running
time (per observation) to fit a neural network to the Kin8nm UCI regression
dataset. Each point corresponds to a different trial. Right panel shows the
RMedSE (root median squared error) vs percentage probability of outlier,
$p_{\epsilon}$ when fitting a neural network to the Kin8nm UCI regression
dataset. (We show the bootstrp average over 100 trials and 95% confidence
interval) Crucially, the same hyper-parameters are used for all experiments
(and are estimated under $p_{\epsilon}=0.1$).

####  E.2.3 Lorenz96 model

In this subsection, we evaluate the EnKF methods in the Lorenz96 model
presented in subsection 4.3 when we only have $N=20$ particles but $d=100$
states. In this scenario, the EnKF is usually modified to incorporate the
covariance inflation proposed in Anderson (1999). We present the results for a
single run of the modified methods in the left panel of Figure 17. In the