# Robust fixed-lag smoothing under model perturbations

Shenglun Yi<EMAIL_ADDRESS>Mattia Zorzi<EMAIL_ADDRESS>Department of Information Engineering, University of Padova, Via Gradenigo
6/B, 35131 Padova, Italy

###### Abstract

A robust fixed-lag smoothing approach is proposed in the case there is a
mismatch between the nominal model and the actual model. The resulting robust
smoother is characterized by a dynamic game between two players: one player
selects the least favorable model in a prescribed ambiguity set, while the
other player selects the fixed-lag smoother minimizing the smoothing error
with respect to least favorable model. We propose an efficient implementation
of the proposed smoother. Moreover, we characterize the corresponding least
favorable model over a finite time horizon. Finally, we test the robust fixed-
lag smoother in two examples. The first one regards a target tracking problem,
while the second one regards a parameter estimation problem.

###### keywords:

Robust fixed-lag smoothing; minimax problem; reduced order smoothing; least
favorable model.

##  1 Introduction

Fixed-lag smoothing aims to estimate the state of a dynamical system at time
$t$ using the observations in the interval $[0,t+L-1]$ with $L>1$. This
algorithm is fundamental in various applications, e.g. tracking and navigation
because it can handle online requirements, see [1, 2, 3]. Moreover, it can be
used in the expectation maximization (EM) algorithm to compute the maximum
likelihood (ML) estimator of the unknown parameters characterizing the
matrices of a state space model, see e.g. [4, 5, 6]. Indeed, although the EM
algorithm is based on the Rauch-Tung-Striebel (RTS) smoother, [7], its
estimate can be approximated by the one given by the fixed-lag smoother
provided that $L$ is taken sufficiently large. However, in all the
aforementioned applications, the actual model is typically known only
imprecisely that is only the nominal model is known. In this situation, these
smoothers could perform poorly.

Model uncertainty is traditionally addressed by risk sensitive filtering, see
e.g. [8, 9, 10, 11, 12, 13] . Here, the robust estimator minimizes an
exponential loss function which severely penalizes large errors. These filters
can be also interpreted as the solution of a dynamic minimax game [14], see
also [15, 16, 17, 18]. Then, [19] proposed a robust estimator which solves an
incremental minimax game. At time $t$ the actual model belongs to the
ambiguity set which is a ball, in the Kullback-Leibler (KL) topology, about
the nominal model. In this way, the uncertainty is “spread” along the time and
not concentrated in specific time steps. Then, there are two players which
operate against. One player, say nature, selects the least favorable model in
this prescribed “ball”, and the other player designs the optimum estimator for
the least favorable model. It is worth noting that many extensions of this
paradigm have been proposed such as: the case with different ambiguity sets
[20, 21, 22, 23]; the distributed case [24, 25]; the case with external input
[26]; the case of degenerate densities [27, 28].

In the literature, the robust smoothing problems mainly consider two
situations. In the first case, the noise distribution is known but it is not
necessarily Gaussian [29, 30, 31, 32, 33, 34, 35, 36, 37, 38], for instance
the noise process is assumed to have a non-Gaussian distribution in order to
model outliers, temporary model uncertainties, missing observations or sensor
delays. Some of these robust paradigms are adaptive because the parameters of
the noises characterizing the state space model are inferred from the
collected data. In the second situation, the noise distribution is not known
but this process takes values in a bounded set, e.g. an ellipsoidal set [39,
40]. However, there are relatively few studies on robust smoothing problems
which use the risk sensitive philosophy, see [41, 42].

In this paper, we propose a new robust fixed-lag smoothing problem where the
model uncertainty is expressed incrementally as in [19, 27]. Thus, at each
time step we have to solve a dynamic game between two players: the nature
which selects the least favorable model in the ambiguity set and the other
player which designs the optimal fixed-lag smoother according to the least
favorable model. The resulting smoother is characterized by matrices whose
dimension is proportional to the lag $L$. On the other hand, the typical value
of the lag is large. Accordingly, numerical instabilities and high
computational burden characterize the algorithm. In order to avoid those
issues, we propose an efficient implementation drawing inspiration from the
reduced order fixed-lag smoother proposed in [43]. Then, the corresponding
least favorable model over a finite simulation horizon is derived to evaluate
the performance of the smoother. Finally, we consider a target tracking
problem and a parameter estimation problem to test the performance of the
robust fixed-lag smoother.

The outline of the paper is as follows. In Section 2, we introduce the problem
formulation about robust fixed-lag smoothing. In Section 3 we derive the
robust fixed-lag smoother. The corresponding algorithm is then reformulated in
an efficient way to reduce the computational burden in Section 4. In Section
5, we derive the least favorable model corresponding to the robust fixed-lag
smoother. The numerical examples are provided in Section 6, which is devoted
to the target tracking problem, and in Section 7, which is devoted to the
parameter estimation problem. Finally, in Section 8 we draw the conclusions.

##  2 Problem Formulation

We consider the nominal state space model:

$\begin{array}[]{cc}x_{t+1}&=Ax_{t}+Bv_{t}\\\ y_{t}&=Cx_{t}+Dv_{t}\end{array}$
(1)

where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times(m+n)}$,
$C\in\mathbb{R}^{m\times n}$ and $D\in\mathbb{R}^{m\times(m+n)}$, $x_{t}$ is
the state vector, $y_{t}$ is the observation vector, and
$v_{t}\in\mathbb{R}^{m+n}$ is normalized white Gaussian noise. Moreover,
$x_{0}\sim\mathcal{N}\left(\hat{x}_{0},V_{0}\right)$, with $V_{0}>0$, which is
independent from $v_{t}$. We also assume that $BD^{\top}=0$,
$\mathrm{rank}(B)=n$ and $\mathrm{rank}(D)=m$. In this way, the process noise,
say $\mathrm{w}_{t}:=Bv_{t}$, and the measurement noise, say
$\mathrm{v}_{t}:=Dv_{t}$, are independent and their covariance matrices are
$BB^{\top}$ and $DD^{\top}$, respectively. Fixed-lag smoothing aims to find an
estimate of $\hat{x}_{t-L+1|t}$ of $x_{t-L+1}$ given
$Y_{t}=\left\\{{y_{0}}\cdots y_{t}\right\\}$ and $L$ denotes the lag. It is
well-known that such a problem can be interpreted as a Kalman prediction
problem corresponding to the following nominal augmented state space model
[43]:

$\begin{array}[]{rl}\xi_{t+1}&=\tilde{A}\xi_{t}+\tilde{B}v_{t}\\\
y_{t}&=\tilde{C}\xi_{t}+\tilde{D}v_{t}\\\
x_{t-L+1}&=\tilde{H}\xi_{t+1}\end{array}$ (2)

where $\tilde{A}\in\mathbb{R}^{(L+1)n\times(L+1)n}$,
$\tilde{B}\in\mathbb{R}^{(L+1)n\times(m+n)}$,
$\tilde{C}\in\mathbb{R}^{m\times(L+1)n}$,
$\tilde{D}\in\mathbb{R}^{m\times(m+n)}$ and
$\tilde{H}\in\mathbb{R}^{n\times(L+1)n}$ are such that

$\xi_{t}=\left[\begin{array}[]{c}{x_{t}}\\\ {x_{t-1}}\\\ {\vdots}\\\
{x_{t-L}}\end{array}\right],~{}\tilde{A}=\left[\begin{array}[]{cccc}{A}&{0}&{\cdots}&{0}\\\
{I}&{\cdots}&{0}&{0}\\\ {\vdots}&{\ddots}&{\vdots}&{\vdots}\\\
{0}&{\cdots}&{I}&{0}\end{array}\right],~{}\tilde{B}=\left[\begin{array}[]{c}{B}\\\
{0}\\\ {\vdots}\\\ {0}\end{array}\right]$
$\tilde{C}=\left[\begin{array}[]{llll}{C}&{0}&{\cdots}&{0}\end{array}\right],\tilde{D}=D,\tilde{H}=\left[\begin{array}[]{llll}{0}&{0}&{\cdots}&{I}\end{array}\right].$

Then,

${\hat{x}_{t-L+1|t}}=\tilde{H}\hat{\xi}_{t+1}$

where $\hat{\xi}_{t+1}$ is the one step-ahead predictor of $\xi_{t+1}$ given
$Y_{t}$. We define $z_{t}:=[\,\xi_{t+1}^{\top}\;y_{t}^{\top}\,]^{\top}.$ Let
$\phi(z_{t}|\xi_{t})$ be the transition probability density function of
$z_{t}$ given $\xi_{t}$ corresponding to the nominal model (2). Then,
$\phi_{t}(z_{t}|\xi_{t})\sim\mathcal{N}\left(m_{z_{t}},K_{z_{t}}\right)$ with

$m_{z_{t}}=\left[\begin{array}[]{c}{\tilde{A}}\\\
{\tilde{C}}\end{array}\right]\xi_{t},~{}~{}~{}K_{z_{t}}=\left[\begin{array}[]{cc}{\tilde{B}\tilde{B}^{\top}}&{0}\\\
{0}&{\tilde{D}\tilde{D}^{\top}}\end{array}\right].$

Notice that $\phi_{t}(z_{t}|\xi_{t})$ is a degenerate probability density
function because $\tilde{B}\tilde{B}^{\top}$ is singular and thus $K_{z_{t}}$
as well. More precisely, $\mathrm{rank}(K_{z_{t}})=n+m$. Accordingly, the
support of $\phi_{t}(z_{t}|\xi_{t})$ is the $n+m$-dimensional affine subspace

$\mathcal{A}_{t}=\left\\{m_{z_{t}}+v,\quad
v\in\mathrm{Im}\left(K_{z_{t}}\right)\right\\}$

which depends on $\xi_{t}$. Then,

$\displaystyle\phi_{t}(z_{t}|\xi_{t})=\left[(2\pi)^{n+m}\operatorname{det}^{+}\left(K_{z_{t}}\right)\right]^{-1/2}\times\exp\left[-\frac{1}{2}\left(z_{t}-m_{z_{t}}\right)^{\top}K_{z_{t}}^{+}\left(z_{t}-m_{z_{t}}\right)\right]$

where $K^{+}_{z_{t}}$ is the pseudo-inverse of $K_{z_{t}}$ and
$\operatorname{det}^{+}(K_{z_{t}})$ is the pseudo-determinant of $K_{z_{t}}$.

The nominal model in (2) in the time horizon $[0,N]$ is described by the joint
probability density

$f\left(\Xi_{N+1},Y_{N}\right)={\tilde{f}_{0}\left(\xi_{0}\right)}\prod_{t=0}^{N}\phi_{t}\left(z_{t}|\xi_{t}\right)$
(3)

where

$\displaystyle\Xi_{N+1}^{\top}$
$\displaystyle=\left[\begin{array}[]{ccccc}{\xi_{0}^{\top}}&{\ldots}&{\xi_{t}^{\top}}&{\ldots}&{\xi_{N+1}^{\top}}\end{array}\right]$
$\displaystyle Y_{N}^{\top}$
$\displaystyle=\left[\begin{array}[]{ccccc}{y_{0}^{\top}}&{\ldots}&{y_{t}^{\top}}&{\ldots}&{y_{N}^{\top}}\end{array}\right],$

${\tilde{f}_{0}(\xi_{0})}\sim\mathcal{N}(\hat{\xi}_{0},\tilde{V}_{0})$ with

$\displaystyle\hat{\xi}_{0}=\left[\begin{array}[]{c}\hat{x}_{0}\\\
\star\end{array}\right],~{}\tilde{V}_{0}=\left[\begin{array}[]{cc}V_{0}&0\\\
0&\star\end{array}\right].$ (8)

In the above equations the star symbol means that it is an arbitrary vector or
matrix. Indeed, since $x_{0}$ is independent from $x_{t}$ with $-L\leq t<0$,
the smoother does not depend on those parameters for $t\geq 0$. Accordingly,
without loss of generality, we assume that $\tilde{V}_{0}>0$.

In this paper, we consider the situation in which the actual model does not
coincide with the nominal one in (2). In particular, we assume that the
probability density of the actual model has a structure similar to the one in
(3):

$\tilde{f}\left(\Xi_{N+1},Y_{N}\right)=\tilde{f}_{0}\left(\xi_{0}\right)\prod_{t=0}^{N}\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)$

where we assume that $\tilde{\phi}_{t}$ has the same support of $\phi_{t}$. In
this way, we can measure the discrepancy between $f$ and $\tilde{f}$ through
the KL-divergence:

$\displaystyle
D(\tilde{f}_{Z},f_{Z}):=\int_{\mathcal{A}_{\Xi}}\int_{\mathcal{A}_{Y}}\tilde{f}\left(\Xi_{N+1},Y_{N}\right)\ln\frac{\tilde{f}\left(\Xi_{N+1},Y_{N}\right)}{f\left(\Xi_{N+1},Y_{N}\right)}dY_{N}d\Xi_{N+1}$

where $\tilde{E}[\cdot]$ is the expected value operator with respect to the
actual probability density $\tilde{f}$ and
$\mathcal{A}_{\Xi}\times\mathcal{A}_{Y}$ is the support of $f$ and
$\tilde{f}$. It is worth noting that the KL divergence is the natural metric
to measure such a mismatch in the case the nominal model is inferred from
data, see [44]. It is not difficult to see that

$\displaystyle
D(\tilde{f},f)=\sum_{t=0}^{N}D\left(\tilde{\phi}_{t},\phi_{t}\right)$ (9)

where

$\displaystyle D($
$\displaystyle\tilde{\phi}_{t},\phi_{t})=\tilde{\mathbb{E}}\left[\ln\left(\frac{\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)}{\phi_{t}\left(z_{t}|\xi_{t}\right)}\right)\right]$
$\displaystyle:=\int_{\bar{\mathcal{A}}_{t}}\int_{\mathcal{A}_{t}}\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}\right)\ln\left(\frac{\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)}{\phi_{t}\left(z_{t}|\xi_{t}\right)}\right)dz_{t}d\xi_{t}$

where $\bar{\mathcal{A}}_{t}$ is the support of $\tilde{f}_{t}(\xi_{t})$ which
denotes the actual marginal density of $\xi_{t}$.

Since the actual model is not known, we assume the latter belongs to the
ambiguity set which is a ball about $f$ formed by placing an upper bound on
$D(\tilde{f},f)$. However, this ambiguity set contains models which
concentrate the uncertainty in a unique time step, i.e. a situation which is
unrealistic in practice. On the other hand, in view of Equation (9), we can
express such a mismatch incrementally through $\phi_{t}$ and
$\tilde{\phi}_{t}$. Accordingly, we assume that $\tilde{\phi}_{t}$ given
$Y_{t-1}$ belongs to the following ambiguity set:

$\mathcal{B}_{t}:=\left\\{\,\tilde{\phi}_{t}\hbox{ s.t.
}\tilde{\mathbb{E}}\left[\ln\left(\frac{\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)}{\phi_{t}\left(z_{t}|\xi_{t}\right)}\right)\bigg{|}Y_{t-1}\right]\leq
c_{t}\right\\}$

where

$\displaystyle\tilde{\mathbb{E}}\left[\ln\left(\frac{\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)}{\phi_{t}\left(z_{t}|\xi_{t}\right)}\right)\bigg{|}Y_{t-1}\right]$
(10)
$\displaystyle:=\int_{\check{\mathcal{A}}_{t}}\int_{\mathcal{A}_{t}}\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)\ln\left(\frac{\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)}{\phi_{t}\left(z_{t}|\xi_{t}\right)}\right)dz_{t}d\xi_{t},$

and $\check{A}_{t}$ is the support of $\tilde{f}_{t}(\xi_{t}|Y_{t-1})$. It is
worth noting that $c_{t}>0$, hereafter called tolerance, is the mismodeling
budget allowed at time step $t$. Our aim is to address the following problem.

###### Problem 1

Design a fixed-lag smoother with respect to the ambiguity set
$\mathcal{B}_{t}$ for $t=0\ldots N$.

##  3  Robust smoothing

We propose a robust fixed-lag smoother of $x_{t-L+1}$ given $Y_{t}$ with
respect to $\mathcal{B}_{t}$ solving the following minimax problem:

$\displaystyle\hat{x}_{t-L+1|t}=\tilde{H}\hat{\xi}_{t+1}$
$\displaystyle\hat{\xi}_{t}=\underset{g_{t}\in\mathcal{G}_{t}}{\mathrm{argmin}}\max_{\tilde{\phi}_{t}\in\mathcal{B}_{t}}J_{t}(\tilde{\phi}_{t},g_{t})$
(11)

where

$\displaystyle J_{t}(\tilde{\phi}_{t},g_{t})=$
$\displaystyle\frac{1}{2}\tilde{\mathbb{E}}\left[\|\tilde{H}\left(\xi_{t+1}-g_{t}\left(y_{t}\right)\right)\|^{2}|Y_{t-1}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{\check{A}_{t}}\int_{\mathcal{A}_{t}}\|\tilde{H}\left(\xi_{t+1}-g_{t}\left(y_{t}\right)\right)\|^{2}\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\times\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)dz_{t}d\xi_{t},$

$\mathcal{G}_{t}$ denotes the class of estimators with finite second-order
moments with respect to all the densities
$\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}(\xi_{t}|Y_{t-1})$
such that $\tilde{\phi}_{t}\in\mathcal{B}_{t}$. Notice that $\tilde{\phi}_{t}$
must satisfy the constraint:

$I_{t}(\tilde{\phi}_{t})\triangleq\int_{\check{A}_{t}}\int_{\mathcal{A}_{t}}\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)dz_{t}d\xi_{t}=1.$
(12)

It is worth noting that problems like (11) can be written as a risk-sensitive
problem, i.e. as a minimization problem where the standard quadratic cost
function is replaced by an exponential cost function, see [45] for more
details.

###### Lemma 1

For a fixed estimator $g_{t}\in\mathcal{G}_{t}$, the density
$\tilde{\phi}_{t}\left(z_{t}|\xi_{t}\right)\in\mathcal{B}_{t}$ that maximizes
the objective function

$J_{t}(\tilde{\phi}_{t},g_{t})=\tilde{\mathbb{E}}\left[\|\tilde{H}\left(\xi_{t+1}-g_{t}\left(y_{t}\right)\right)\|^{2}|Y_{t-1}\right]$

under constraint $D_{t}(\tilde{\phi}_{t},\phi_{t})\leq c_{t}$ is given by

$\tilde{\phi}_{t}^{0}=\frac{1}{M_{t}\left(\lambda_{t}\right)}\exp\left(\frac{1}{2\lambda_{t}}\left\|\tilde{H}(x_{t+1}-g_{t}\left(y_{t})\right)\right\|^{2}\right)\phi_{t}$
(13)

where $M_{t}(\lambda_{t})$ is the normalizing constant defined as follows:

$\displaystyle M_{t}\left(\lambda_{t}\right)=$
$\displaystyle\int_{\check{\mathcal{A}}_{t}}\int_{\mathcal{A}_{t}}\exp\left(\frac{1}{2\lambda_{t}}\|\tilde{H}(\xi_{t+1}-g_{t}\left(y_{t})\right)\|^{2}\right)\phi_{t}\times\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)dz_{t}d\xi_{t}.$

Moreover, for $c_{t}>0$ sufficiently small, there exists a unique
$\lambda_{t}>0$ such that $D(\tilde{\phi}^{0}_{t},\phi_{t})=c_{t}$.

###### Proof 1

The proof is similar to the one of [27, Lemma2].

Once we get the function $\tilde{\phi}^{0}_{t}$, the estimator
$g_{t}\in\mathcal{G}_{t}$ minimizing the objective function
$J_{t}(\tilde{\phi}_{t}^{0},g_{t})$ is given by

$\displaystyle\hat{\xi}_{t+1}=g_{t}^{0}\left(y_{t}\right)=\tilde{\mathbb{E}}\left[\xi_{t+1}|Y_{t}\right]=\int_{\check{\mathcal{A}}_{t+1}}\xi_{t+1}\tilde{f}_{t+1}\left(\xi_{t+1}|Y_{t}\right)d\xi_{t+1}$

where

$\tilde{f}_{t+1}\left(\xi_{t+1}|Y_{t}\right)=\frac{\int_{\check{\mathcal{A}_{t}}}\tilde{\phi}_{t}^{0}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)d\xi_{t}}{\int_{\check{\mathcal{A}}_{t}}\int_{\mathcal{A}_{t}^{\star}}\tilde{\phi}_{t}^{0}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)d\xi_{t+1}d\xi_{t}}$
(14)

where $\mathcal{A}^{\star}_{t}$ is defined as follows:
$\xi_{t+1}\in\mathcal{A}^{\star}_{t}$ if and only if there exists at least one
$y_{t}$ for which
$[\,\xi_{t+1}^{\top}\;y_{t}^{\top}\,]^{\top}\in\mathcal{A}_{t}$. The optimal
estimator $g_{t}^{0}$ solution to (11) relies on the least-favorable density
$\tilde{\phi}_{t}^{0}\left(z_{t}|\xi_{t}\right)$. On the other hand, the
latter depends on the estimator $g^{0}_{t}$. In order to break this deadlock
problem, an additional assumption is needed. More precisely, we assume that
the $a~{}priori$ conditional density $\tilde{f}_{t}(\xi_{t}|Y_{t-1})$ is
Gaussian
$\tilde{f}_{t}(\xi_{t}|Y_{t-1})\sim\mathcal{N}(\hat{\xi}_{t},\tilde{V}_{t})$.
In view of (2), the marginal density

$\bar{f}_{t}\left(z_{t}|Y_{t-1}\right):=\int_{\mathcal{A}_{t}}\phi_{t}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)d\xi_{t}$

is Gaussian so that

$\bar{f}_{t}\left(z_{t}|Y_{t-1}\right)\sim\mathcal{N}\left(m_{z_{t}|Y_{t-1}},K_{z_{t}|Y_{t-1}}\right)$
(15)

where

$\displaystyle m_{z_{t}|Y_{t-1}}$
$\displaystyle=\left[\begin{array}[]{c}{\tilde{A}}\\\
{\tilde{C}}\end{array}\right]\hat{\xi}_{t},$ $\displaystyle K_{z_{t}|Y_{t-1}}$
$\displaystyle=\left[\begin{array}[]{c}{\tilde{A}}\\\
{\tilde{C}}\end{array}\right]\tilde{V}_{t}\left[\begin{array}[]{cc}{\tilde{A}^{\top}}&{\tilde{C}^{\top}}\end{array}\right]+\left[\begin{array}[]{c}{\tilde{B}}\\\
{\tilde{D}}\end{array}\right]\left[\begin{array}[]{cc}{\tilde{B}^{\top}}&{\tilde{D}^{\top}}\end{array}\right].$

Then, on the basis of Lemma 1, we have that the least favorable density of
$z_{t}$ given $Y_{t-1}$ is

$\displaystyle\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right):=\int_{\mathcal{A}_{t}}\tilde{\phi}_{t}^{0}\left(z_{t}|\xi_{t}\right)\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)d\xi_{t}$
(16)
$\displaystyle=\frac{1}{M_{t}\left(\lambda_{t}\right)}\exp\left(\frac{1}{2\lambda_{t}}\|\tilde{H}\left(\xi_{t+1}-g_{t}\left(y_{t}\right)\right)\|^{2}\right)\bar{f}_{t}\left(z_{t}|Y_{t-1}\right).$

Accordingly, $\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right)$ is a Gaussian
probability density.

###### Lemma 2

Consider the state space model (1) with $\mathrm{rank}(B)=n$ and
$\mathrm{rank}(D)=m$. If $\tilde{V}_{t}>0$ then $K_{z_{t}|Y_{t-1}}>0$, i.e.
$\bar{f}_{t}\left(z_{t}|Y_{t-1}\right)$ is a non-degenerate density. Moreover,

$\displaystyle\tilde{P}_{t+1}:=\tilde{A}\tilde{V}_{t}\tilde{A}^{\top}-\tilde{A}\tilde{V}_{t}\tilde{C}^{\top}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}\tilde{V}_{t}\tilde{A}^{\top}+\tilde{B}\tilde{B}^{\top}$

is positive definite.

###### Proof 2

First, notice that the block in position (2,2) of $K_{z_{t}|Y_{t-1}}$ is
$\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top}$. The latter
is positive definite because
$\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top}\geq\tilde{D}\tilde{D}^{\top}>0$.
Accordingly, in order to prove that $K_{z_{t}|Y_{t-1}}>0$ it is sufficient to
prove that $\tilde{P}_{t+1}$, which is the Schur complement of the block (2,2)
of $K_{z_{t}|Y_{t-1}}$, is positive definite. Since $\tilde{V}_{t}$ is
invertible we can rewrite $\tilde{P}_{t+1}$ by using the Woodbury formula:

$\displaystyle\tilde{P}_{t+1}=\tilde{A}[\tilde{V}_{t}^{-1}+\tilde{C}^{\top}(\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}]^{-1}\tilde{A}^{\top}+\tilde{B}\tilde{B}^{\top}.$

Let $v=[\,v_{1}\;v_{2}\,]^{\top}\in\mathbb{R}^{(n+1)L}$ be such that
$v_{1}\in\mathbb{R}^{n}$ and $v_{2}\in\mathbb{R}^{nL}$. Notice that

$\displaystyle
v^{\top}\tilde{P}_{t+1}v=v^{\top}\tilde{A}[\tilde{V}_{t}^{-1}+\tilde{C}^{\top}(\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}]^{-1}\tilde{A}^{\top}v+v_{1}^{\top}BB^{\top}v_{1}$

and thus

$\displaystyle v^{\top}\tilde{P}_{t+1}v$ $\displaystyle\geq
v_{1}^{\top}BB^{\top}v_{1}$ (17) $\displaystyle v^{\top}\tilde{P}_{t+1}v$
$\displaystyle\geq
v^{\top}\tilde{A}[\tilde{V}_{t}^{-1}+\tilde{C}^{\top}(\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}]^{-1}\tilde{A}^{\top}v.$
(18)

Assume that $v^{\top}\tilde{P}_{t+1}v=0$. Since $BB^{\top}>0$, by (17) we have
that $v_{1}=0$. Accordingly, the inequality in (18) for $v_{1}=0$ becomes

$\displaystyle v^{\top}\tilde{P}_{t+1}v$
$\displaystyle\geq\left[\begin{array}[]{cc}0&v_{2}^{\top}\end{array}\right]\tilde{A}[\tilde{V}_{t}^{-1}+\tilde{C}^{\top}(\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}]^{-1}\tilde{A}^{\top}\left[\begin{array}[]{c}0\\\
v_{2}\end{array}\right].$ (22)

Moreover, in view of the particular structure of $\tilde{A}$, we have that
$[\,0\;v_{2}^{\top}\,]\tilde{A}=[\,v_{2}^{\top}\;0\,]=0$ if and only if
$v_{2}=0$. Since
$[\tilde{V}_{t}^{-1}+\tilde{C}^{\top}(\tilde{D}\tilde{D}^{\top})^{-1}\tilde{C}]^{-1}>0$,
because $\tilde{V}_{t}>0$, it follows that $v_{2}=0$ and thus $v=0$. We proved
that if $v^{\top}\tilde{P}_{t+1}v=0$ then $v=0$, i.e. $\tilde{P}_{t+1}$ is
positive definite.

Finally, if $\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)$ is Gaussian then, in
view of (14), also $\tilde{f}_{t+1}\left(\xi_{t+1}|Y_{t}\right)$ is Gaussian.
Accordingly, the assumption that $\tilde{f}_{0}\left(\xi_{0}\right)$ is
Gaussian, implies that $\tilde{f}_{t}\left(\xi_{t}|Y_{t-1}\right)$ is Gaussian
for any $t$.

###### Theorem 3

Consider the state space model (1) where we recall that $\mathrm{rank}(B)=n$
and $\mathrm{rank}(D)=m$. Let
$\tilde{f}_{t}(\xi_{t}|Y_{t-1})\sim\mathcal{N}(\hat{\xi}_{t},\tilde{V}_{t})$
with $\tilde{V}_{t}>0$. Then, the estimator

${g}_{t}^{0}\left(y_{t}\right)=\tilde{A}\hat{\xi}_{t}+\tilde{G}_{t}(y_{t}-\tilde{C}\hat{\xi}_{t})$
(23)

with

$\displaystyle\tilde{G}_{t}$
$\displaystyle=\tilde{A}\tilde{V}_{t}\tilde{C}^{\top}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top})^{-1}$

solves Problem (11). The nominal error covariance of $\xi_{t+1}$ given $y_{t}$
is

$\displaystyle\tilde{P}_{t+1}$
$\displaystyle=\tilde{A}\tilde{V}_{t}\tilde{A}^{\top}-\tilde{G}_{t}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top}){\tilde{G}_{t}}^{\top}+\tilde{B}\tilde{B}^{\top}$

and the perturbed error covariance of $\xi_{t+1}$ given $Y_{t}$ is

$\tilde{V}_{t+1}=(\tilde{P}_{t+1}^{-1}-\lambda^{-1}\tilde{H}^{\top}\tilde{H})^{-1}.$

If we denote $r(P)$ as the largest eigenvalue of $P$, then the Lagrange
multiplier ${\lambda_{t}}>r(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top})$ is
unique and such that

$\displaystyle\gamma\left(\lambda_{t}\right)=\frac{1}{2}\left[\operatorname{tr}\left((I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}-I\right)+\ln\operatorname{det}(I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})\right]=c_{t}.$
(24)

Finally, the least favorable density $\tilde{f}_{t}^{0}(z_{t}|Y_{t-1})$
corresponding to the solution of (11) is a non-degenerate Gaussian density.

###### Proof 3

As we already noticed, $\bar{f}_{t}(z_{t}|Y_{t-1})$ is Gaussian, and in view
of Lemma 2, non-degenerate. Accordingly, in view of (16),
$\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right)$ is Gaussian and non-degenerate. Let

$\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right)\sim\mathcal{N}(\tilde{m}_{z_{t}|Y_{t-1}},\tilde{K}_{z_{t}|Y_{t-1}})$

with

$\tilde{m}_{z_{t}|Y_{t-1}}=\left[\begin{array}[]{c}{\tilde{m}_{\xi_{t+1}|Y_{t-1}}}\\\
{\tilde{m}_{y_{t}|Y_{t-1}}}\end{array}\right],$

and

$\tilde{K}_{z_{t}|Y_{t-1}}=\left[\begin{array}[]{cc}{\tilde{K}_{\xi_{t+1}|Y_{t-1}}}&{\tilde{K}_{\xi_{t+1}y_{t}|Y_{t-1}}}\\\
{\tilde{K}_{y_{t}\xi_{t+1}|Y_{t-1}}}&{\tilde{K}_{y_{t}|Y_{t-1}}}\end{array}\right].$

In view of (15) and (16), the conditional KL-divergence in (10) admits the
closed-form expression

$\displaystyle\tilde{\mathbb{E}}$
$\displaystyle\left[\ln({\tilde{\phi}_{t}^{\circ}}/\phi_{t})|Y_{t-1}\right]=\tilde{\mathbb{E}}\left[\ln(\tilde{f}_{t}/\bar{f}_{t})|Y_{t-1}\right]$
(25) $\displaystyle=$ $\displaystyle\frac{1}{2}\left[\left\|\Delta
m\right\|_{K_{z_{t}|Y_{t-1}}^{-1}}^{2}+\operatorname{tr}(K_{z_{t}|Y_{t-1}}^{-1}\tilde{K}_{z_{t}|Y_{t-1}}-I)-\ln\operatorname{det}(K_{z_{t}|Y_{t-1}}^{-1}\tilde{K}_{z_{t}|Y_{t-1}})\right]$

where

$\Delta m=\tilde{m}_{z_{t}|Y_{t-1}}-m_{z_{t}|Y_{t-1}}.$

Therefore, we can rewrite the minimax game (11) with respect to
$\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right)$ and $\bar{f}_{t}(z_{t}|Y_{t-1})$:

$(g_{t}^{0},\tilde{f}^{0}_{t})=\arg\min_{g_{t}\in\mathcal{\bar{G}}_{t}}\max_{\tilde{f}_{t}\in\mathcal{\bar{B}}_{t}}J_{t}(\tilde{f}_{t},g_{t})$

where

$\mathcal{\bar{B}}_{t}=\left\\{\tilde{f}_{t}(z_{t}|Y_{t-1})~{}s.t.~{}\tilde{\mathbb{E}}\left[\ln(\tilde{f}_{t}/\bar{f}_{t})|Y_{t-1}\right]\leq
c_{t}\right\\},$

$\mathcal{\bar{G}}_{t}$ is the set of estimators with finite second order
moments with respect to all the densities
$\tilde{f}_{t}\in\mathcal{\bar{B}}_{t}$ and

$\displaystyle
J_{t}(\tilde{f}_{t},g_{t}):=\frac{1}{2}\int_{\mathbb{R}^{(L+1)n+m}}\|\tilde{H}(\xi_{t+1}$
$\displaystyle-
g_{t}(y_{t}))|Y_{t-1}\|^{2}\times\tilde{f}_{t}\left(z_{t}|Y_{t-1}\right)dz_{t}.$

Next, we prove that $\tilde{f}^{0}_{t}$ and $g^{0}_{t}$ are such that

$J_{t}(\tilde{f}_{t},g_{t}^{0})\leq J_{t}(\tilde{f}_{t}^{0},g_{t}^{0})\leq
J_{t}(\tilde{f}_{t}^{0},g_{t}),$ (26)

where
$\tilde{f}^{0}_{t}\left(z_{t}|Y_{t-1}\right)\sim\mathcal{N}(\tilde{m}^{0}_{z_{t}|Y_{t-1}},\tilde{K}^{0}_{z_{t}|Y_{t-1}})$
with

$\displaystyle\tilde{m}^{0}_{z_{t}|Y_{t-1}}=m_{z_{t}|Y_{t-1}},$ (27)
$\displaystyle\tilde{K}^{0}_{z_{t}|Y_{t-1}}=\left[\begin{array}[]{cc}{\tilde{K}_{\xi_{t+1}|Y_{t-1}}}&{K_{\xi_{t+1}y_{t}|Y_{t-1}}}\\\
{K_{y_{t}\xi_{t+1}|Y_{t-1}}}&{K_{y_{t}|Y_{t-1}}}\end{array}\right].$

Since $\tilde{f}^{0}_{t}$ is Gaussian, the optimal estimator satisfying the
second inequality in (26) is (23). Then, it remains to prove that the least
favorable density $\tilde{f}^{0}_{t}$ is such that (27) holds.

It is not difficult to see that

$\displaystyle J_{t}(\tilde{f}_{t},g_{t}^{0})=$
$\displaystyle\frac{1}{2}\operatorname{tr}\left\\{\left[\begin{array}[]{c}{I}\\\
{-\tilde{G}_{t}^{\top}}\end{array}\right]\tilde{H}^{\top}\tilde{H}\left[I\quad-\tilde{G}_{t}\right]\times\left(\tilde{K}_{z_{t}|Y_{t-1}}+\Delta
m\Delta m^{\top}\right)\right\\}.$ (28)

Then, based on the parametric structure of the KL divergence in (25) and the
objective function in (28), we consider the corresponding Lagrangian as a
function of $\tilde{m}_{z_{t}|Y_{t-1}}$ and $\tilde{K}_{z_{t}|Y_{t-1}}$ as
follows:

$\displaystyle\mathcal{L}($
$\displaystyle\tilde{m}_{z_{t}|Y_{t-1}},\tilde{K}_{z_{1}|Y_{t-1}},\lambda_{t})$
$\displaystyle=$ $\displaystyle
J(\tilde{m}_{z|Y_{t-1}},\tilde{K}_{z_{t}|Y_{t-1}})+\lambda_{t}\left(c_{t}-\tilde{\mathbb{E}}[\ln(\tilde{f}_{t}/\bar{f}_{t})|Y_{t-1}]\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\operatorname{tr}\left\\{\left[\begin{array}[]{c}{I}\\\
{-{\tilde{G}_{t}}^{\top}}\end{array}\right]\tilde{H}^{\top}\tilde{H}\left[I\quad-{\tilde{G}_{t}}\right](\tilde{K}_{z_{t}|Y_{t-1}}+\Delta
m\Delta m^{\top})\right\\}$ $\displaystyle-\frac{\lambda_{t}}{2}\Delta
m^{\top}K_{z_{t}|Y_{t-1}}^{-1}\Delta
m-\frac{\lambda_{t}}{2}\operatorname{tr}(K_{z_{t}|Y_{t-1}}^{-1}\tilde{K}_{z_{t}|Y_{t-1}})+\frac{\lambda_{t}}{2}\operatorname{tr}(I)$
$\displaystyle+\frac{\lambda_{t}}{2}\operatorname{In}\operatorname{det}(K_{z_{t}|Y_{t-1}}^{-1}\tilde{K}_{z_{t}|Y_{t-1}})+\lambda_{t}c_{t}$
$\displaystyle=$
$\displaystyle\lambda_{t}c_{t}+\frac{\lambda_{t}}{2}\operatorname{tr}(I)+\frac{\lambda_{t}}{2}\ln\operatorname{det}(K_{z_{t}|Y_{t-1}}^{-1}\tilde{K}_{z_{t}|Y_{t-1}})$
$\displaystyle+\frac{1}{2}\operatorname{tr}(W(\lambda_{t})\tilde{K}_{z_{t}|Y_{t-1}})+\frac{1}{2}\Delta
m^{\top}W(\lambda_{t})\Delta m$

where

$W(\lambda_{t})\triangleq\left[\begin{array}[]{c}{I}\\\
{-{\tilde{G}_{t}}^{\top}}\end{array}\right]\tilde{H}^{\top}\tilde{H}\left[\begin{array}[]{cc}{I}&{-{\tilde{G}_{t}}}\end{array}\right]-\lambda_{t}K_{z_{t}|Y_{t-1}}^{-1}.$

The first variation and the second variation of $\mathcal{L}$ with respect to
$\tilde{m}_{z_{t}|Y_{t-1}}$ are, respectively,

$\displaystyle\delta\mathcal{L}$
$\displaystyle(\tilde{m}_{z_{t}|Y_{t-1}},\tilde{K}_{z_{t}|Y_{t-1}},\lambda_{t};\delta\tilde{m}_{z_{t}|Y_{t-1}})$
$\displaystyle~{}~{}~{}~{}~{}=\frac{1}{2}\delta\tilde{m}^{\top}_{z_{t}|Y_{t-1}}W(\lambda_{t})\Delta
m+\Delta m^{\top}W(\lambda_{t})\delta\tilde{m}_{z_{t}|Y_{t-1}},$
$\displaystyle\delta^{2}{\mathcal{L}}(\tilde{m}_{z_{t}|Y_{t-1}},\tilde{K}_{z_{t}|Y_{t-1}},\lambda_{t};\delta\tilde{m}_{z_{t}|Y_{t-1}},\delta\tilde{m}_{z_{t}|Y_{t-1}})=\delta\tilde{m}_{z_{t}|Y_{t-1}}W(\lambda_{t})\delta\tilde{m}_{z_{t}|Y_{t-1}}$

so that $\delta^{2}{\mathcal{L}}<0$, for any $\tilde{m}_{z_{t}|Y_{t-1}}\neq 0$
if and only if $W(\lambda_{t})<0$, which means ${\mathcal{L}}$ is strictly
concave if and only if $W(\lambda_{t})$ is negative definite. Next, we find
the condition on $\lambda_{t}$ for which $W(\lambda_{t})<0$. We denote
$M=\tilde{H}[\,I\;{-\tilde{G}_{t}}\,]O_{z_{t}}$, where
$K_{z_{t}|Y_{t-1}}=O_{z_{t}}O_{z_{t}}^{\top}$, so that

$\displaystyle
MM^{\top}=\tilde{H}\left[\begin{array}[]{cc}{I}&{-\tilde{G}_{t}}\end{array}\right]O_{z_{t}}O_{z_{t}}^{\top}\left[\begin{array}[]{c}{I}\\\
{-\tilde{G}^{\top}_{t}}\end{array}\right]\tilde{H}^{\top}$
$\displaystyle=\tilde{H}\left(K_{\xi_{t+1}|Y_{t-1}}-K_{\xi_{t+1}y_{t}|Y_{t-1}}K_{y_{t}|Y_{t-1}}^{-1}K_{y_{t}\xi_{t+1}|Y_{t-1}}\right)\tilde{H}^{\top}=\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}.$

Therefore,

$\displaystyle W(\lambda_{t})$
$\displaystyle=O^{-\top}_{z_{t}}\left(O_{z_{t}}^{\top}\left[\begin{array}[]{c}{I}\\\
{-\tilde{G}_{t}^{\top}}\end{array}\right]\tilde{H}^{\top}\tilde{H}\left[\begin{array}[]{cc}{I}&-\tilde{G}_{t}\end{array}\right]O_{z_{t}}-\lambda_{t}I\right)O_{z_{t}}^{-1}$
$\displaystyle=O_{z_{t}}^{-\top}\left(M^{\top}M-\lambda
I\right)O_{z_{t}}^{-1}$

and it is congruent to $\bar{W}(\lambda_{t})=M^{\top}M-\lambda I$, which means
$\bar{W}(\lambda_{t})<0$ if and only if $W(\lambda_{t})<0$, meanwhile, it is
not difficult to see that
$r(M^{\top}M)=r(MM^{\top})=r(\tilde{H}P_{t+1}\tilde{H}^{\top})$. Hence,
$W(\lambda_{t})<0$ as long as the Lagrange multiplier
$\lambda_{t}>r(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top})$. In such a situation
the minimum is such that $\Delta m=0$, which implies

$\tilde{m}^{0}_{z_{t}|Y_{t-1}}={m}_{z_{t}|Y_{t-1}}.$ (29)

The first variation and the second variation of $\mathcal{L}$ with respect to
$\tilde{K}_{z_{1}|Y_{t-1}}$ are, respectively,

$\displaystyle\delta L$
$\displaystyle(\tilde{m}_{z_{t}|{Y_{t-1}}},\tilde{K}_{z_{t}|Y_{t-1}},\lambda_{t};\delta\tilde{K}_{z_{t}|Y_{t-1}})$
$\displaystyle~{}~{}~{}~{}~{}=\frac{1}{2}\operatorname{tr}\left\\{{\lambda_{t}}\tilde{K}_{z_{t}|Y_{t-1}}^{-1}\delta\tilde{K}_{z_{t}|Y_{t-1}}+W(\lambda_{t})\delta\tilde{K}_{z_{t}|Y_{t-1}}\right\\}$

and

$\displaystyle\delta^{2}L(\tilde{m}$
$\displaystyle{}_{z_{t}|Y_{t-1}},\tilde{K}_{z_{t}|Y_{t-1}},\lambda_{t};\delta\tilde{K}_{z_{t}|Y_{t-1}},{\delta\tilde{K}_{z_{t}|Y_{t-1}}})$
$\displaystyle~{}~{}~{}~{}~{}=-\frac{1}{2}\lambda_{t}\operatorname{tr}\left\\{(\tilde{K}_{z_{t}|Y_{t-1}}^{-1}\delta\tilde{K}_{z_{t}|Y_{t-1}})^{2}\right\\}<0.$

Accordingly, $\mathcal{L}$ is strictly concave in $\tilde{K}_{z_{t}|Y_{t-1}}$.
Thus, the minimum point $\tilde{K}^{0}_{z_{t}|Y_{t-1}}$ is given by imposing
the stationarity condition

$\delta
L(\tilde{m}_{z_{t}|{Y_{t-1}}},\tilde{K}^{0}_{z_{t}|Y_{t-1}},\lambda_{t};\delta\tilde{K}_{z_{t}|Y_{t-1}})=0.$

The latter implies that such a point is

$\tilde{K}^{0,-1}_{z_{t}|Y_{t-1}}=K^{-1}_{z_{t}|Y_{t-1}}-\frac{1}{\lambda_{t}}\left[\begin{array}[]{c}{I}\\\
{-{\tilde{G}_{t}}^{\top}}\end{array}\right]\tilde{H}^{\top}\tilde{H}\left[\begin{array}[]{cc}{I}&{-\tilde{G}_{t}}\end{array}\right].$
(30)

Notice that the block upper diagonal lower (UDL) form of ${K}_{z_{t}|Y_{t-1}}$
is

$K_{z_{t}|Y_{t-1}}=\left[\begin{array}[]{cc}{I}&{\tilde{G}_{t}}\\\
{0}&{I}\end{array}\right]\left[\begin{array}[]{cc}{\tilde{P}_{t+1}}&{0}\\\
{0}&{K_{y_{t}|Y_{t-1}}}\end{array}\right]\left[\begin{array}[]{cc}{I}&{0}\\\
{{\tilde{G}_{t}}^{\top}}&{I}\end{array}\right]$ (31)

and its inverse admits the following UDL decomposition

$K_{z_{t}|Y_{t-1}}^{-1}=\left[\begin{array}[]{cc}{I}&{0}\\\
{-{\tilde{G}_{t}}^{\top}}&{I}\end{array}\right]\left[\begin{array}[]{cc}{\tilde{P}^{-1}_{t+1}}&{0}\\\
{0}&{K_{y_{t}|Y_{t-1}}^{-1}}\end{array}\right]\left[\begin{array}[]{cc}{I}&{-\tilde{G}_{t}}\\\
{0}&{I}\end{array}\right].$ (32)

Therefore, substituting Equation (32) in Equation (30), we have

$\displaystyle{\tilde{K}_{z_{t}|Y_{t-1}}^{0,-1}}=\left[\begin{array}[]{cc}{I}&{0}\\\
{-{\tilde{G}_{t}}^{\top}}&{I}\end{array}\right]\left[\begin{array}[]{cc}{\tilde{P}_{t+1}^{-1}-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}}&{0}\\\
{0}&{K_{y_{t}|Y_{t-1}}^{-1}}\end{array}\right]\left[\begin{array}[]{cc}{I}&{-{\tilde{G}_{t}}}\\\
{0}&{I}\end{array}\right],$

so that

$\displaystyle{\tilde{K}_{z_{t}|Y_{t-1}}^{0}}=\left[\begin{array}[]{cc}{I}&{\tilde{G}_{t}}\\\
{0}&{I}\end{array}\right]\left[\begin{array}[]{cc}{\tilde{V}_{t+1}}&{0}\\\
{0}&{K_{y_{t}|Y_{t-1}}}\end{array}\right]\left[\begin{array}[]{cc}{I}&{0}\\\
{{\tilde{G}_{t}}^{\top}}&{I}\end{array}\right]$ (33)

where

$\tilde{V}_{t+1}=(\tilde{P}_{t+1}^{-1}-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H})^{-1}.$

Then, let
$\gamma(\lambda_{t}):=\tilde{\mathbb{E}}[\ln(\tilde{f}_{t}^{0}/\bar{f}_{t})|Y_{t-1}].$
By taking into account (25) and using Equations (29), (31) and (33), we obtain
(24). The first derivative of $\gamma(\lambda_{t})$ is

$\displaystyle\frac{\partial\gamma(\lambda_{t};\delta\lambda_{t})}{\partial\lambda_{t}}=$
$\displaystyle\frac{\lambda_{t}^{-2}}{2}\operatorname{tr}\left[(I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}^{\top}\tilde{P}_{t+1})^{-1}\right.$
$\displaystyle\left.\times\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1}(I-(I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1})\right]$
$\displaystyle=$
$\displaystyle\frac{\lambda_{t}^{-2}}{2}\operatorname{tr}\left[(I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1}\right.$
$\displaystyle\left.\times(I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}((I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}^{\top}\tilde{P}_{t+1})-I)\right]$
$\displaystyle=$
$\displaystyle-\frac{\lambda_{t}^{-3}}{2}\operatorname{tr}\left[((I-\lambda_{t}^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{2}\right]<0.$

Therefore, $\gamma(\lambda_{t})$ is strictly monotone decreasing. Moreover, it
is not difficult to see that

$\lim_{\lambda_{t}\rightarrow\infty}\gamma(\lambda_{t})=0,\quad\lim_{\lambda_{t}\rightarrow
r(\tilde{Q})}\gamma(\lambda_{t})=+\infty$ (34)

where $r(\tilde{Q})=r(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top})$. As a
consequence, there exists a unique Lagrangian multiplier
$\lambda_{t}>r(\tilde{Q})>0$ such that $D(\tilde{f},f)=c_{t}$. The fact that
$\tilde{f}_{t}^{0}(z_{t}|Y_{t-1})$ is non-degenerate follows from the fact
that ${\tilde{K}_{z_{t}|Y_{t-1}}^{0}}\geq K_{z_{t}|Y_{t-1}}>0$.

###### Corollary 1

Consider the state space model (1) where we recall: $\mathrm{rank}(B)=n$,
$\mathrm{rank}(D)=m$ and
$\xi_{0}\sim{\mathcal{N}}(\hat{\xi}_{0},\tilde{V}_{0})$ with
$\tilde{V}_{0}>0$. Then, $\tilde{f}_{t}^{0}(z_{t}|Y_{t-1})$ is non-degenerate
for any $t\geq 0$.

###### Proof 4

We prove the claim by induction. Let
$\tilde{f}_{t}(\xi_{t}|Y_{t-1})\sim\mathcal{N}(\hat{\xi}_{t},\tilde{V}_{t})$
with $\tilde{V}_{t}>0$. By Theorem 3 we have that
$\tilde{f}_{t}^{0}(z_{t}|Y_{t-1})$ is Gaussian non-degenerate and $V_{t+1}>0$.
From, $\tilde{f}_{t}(z_{t}|Y_{t-1})$ we have that
$\tilde{f}_{t+1}(\xi_{t+1}|Y_{t})\sim{\mathcal{N}}(\hat{\xi}_{t+1},\tilde{V}_{t+1})$
which is non-degenerate. Finally, at the initial time $t=0$, we have
$\tilde{f}_{0}(\xi_{0}|Y_{-1}):=\tilde{f}_{0}(\xi_{0})\sim\mathcal{N}(\hat{\xi}_{0},\tilde{V}_{0})$
and $\tilde{V}_{0}$ is positive definite by assumption.

The resulting robust fixed-lag smoother is outlined in Algorithm 1 where
$\theta_{t}:=\lambda_{t}^{-1}$ is the risk sensitivity parameter and

$\displaystyle\begin{aligned}
\gamma(\tilde{P}_{t+1},\theta_{t}):=\frac{1}{2}\left[\operatorname{tr}\left((I-\theta_{t}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}-I\right)+\ln\operatorname{det}(I-\theta_{t}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})\right].\end{aligned}$
(35)

Algorithm 1 Robust fixed-lag smoother with lag L

1:$y_{0}\ldots y_{N}$, $\hat{\xi}_{0}$, $\tilde{V}_{0}$,
$\color[rgb]{0,0,0}c_{t}$

2:$\hat{x}_{t-L+1|t}$, $t=L-1\ldots N$

3:for $t=0:N$ do

4:
$\tilde{G}_{t}=\tilde{A}\tilde{V}_{t}\tilde{C}^{\top}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top})^{-1}$

5:
$\hat{\xi}_{t+1}=\tilde{A}\hat{\xi}_{t}+\tilde{G}_{t}(y_{t}-\tilde{C}\hat{\xi}_{t})$

6: $\hat{x}_{t-L+1|t}=\tilde{H}\hat{\xi}_{t+1}$

7:
$\tilde{P}_{t+1}=\tilde{A}\tilde{V}_{t}\tilde{A}^{\top}-\tilde{G}_{t}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top}){\tilde{G}_{t}^{\top}}+\tilde{B}\tilde{B}^{\top}$

8: Find $\theta_{t}$ s.t. $\gamma(\tilde{P}_{t+1},\theta_{t})=c_{t}$

9:
$\tilde{V}_{t+1}=(\tilde{P}^{-1}_{t+1}-\theta_{t}\tilde{H}^{\top}\tilde{H})^{-1}$

10:end for

Finally, in the case that $c_{t}=0$, i.e. the nominal model coincides with the
actual one, it is not difficult to see that $\theta_{t}=0$ that is we obtain
the standard fixed-lag smoother.

###### Remark 1

In the presence of a deterministic input $u_{t}$, then it is possible to
derive the corresponding robust fixed-lag smoother by using arguments similar
to the ones in [26]. For instance, if the input acts only in the state
equation, i.e. we have $x_{t+1}=Ax_{t}+Bv_{t}+u_{t}$, then Step 5 in Algorithm
1 is substituted with
$\hat{\xi}_{t+1}=\tilde{A}\hat{\xi}_{t}+\tilde{G}_{t}(y_{t}-\tilde{C}\hat{\xi}_{t})+w_{t}$,
where $w_{t}:=[\,u_{t}^{\top}\,u_{t-1}^{\top}\ldots
u_{t-L+1}^{\top}\,]^{\top}$.

##  4 Efficient implementation

Algorithm 1 is not numerically robust and efficient in terms of computational
burden. Since the dimension of $\tilde{P}_{t}$ and $\tilde{V}_{t}$ is
proportional to $L$, which is typically large, their inversion is time
consuming and not accurate. Accordingly, there is the need to develop an
efficient strategy which avoids those matrix inversions as it has been done in
[43] for the standard fixed-lag smoother. The efficient procedure for our
robust smoother is outlined in Algorithm 2. Next, we explain how to derive the
salient steps. In what follows we always refer to the steps of Algorithm 2 if
not specified.

First, we rewrite the risk-sensitive Riccati iteration in Step 7 of Algorithm
1 as:

$\displaystyle\tilde{L}_{t}$
$\displaystyle=\tilde{V}_{t}\tilde{C}^{\top}(\tilde{C}\tilde{V}_{t}\tilde{C}^{\top}+\tilde{D}\tilde{D}^{\top})^{-1}$
(36) $\displaystyle\tilde{P}_{t|t}$
$\displaystyle=(I-\tilde{L}_{t}\tilde{C})\tilde{V}_{t}$ (37)
$\displaystyle\tilde{P}_{t+1}$
$\displaystyle=\tilde{A}\tilde{P}_{t|t}\tilde{A}^{\top}+\tilde{B}\tilde{B}^{\top}$
(38)

where
$\tilde{P}_{t|t}:=\tilde{\mathbb{E}}[(\xi_{t}-\hat{\xi}_{t|t})(\xi_{t}-\hat{\xi}_{t|t})^{\top}]$
and $\hat{\xi}_{t|t}$ is the estimator of $\xi_{t}$ given $Y_{t}$. Then, we
parameterize $\tilde{V}_{t}$ and $\tilde{L}_{t}$ in blocks of $n\times n$
matrices as follows:

$\displaystyle\tilde{V}_{t}$
$\displaystyle=\left[\begin{array}[]{cccccc}V_{t}&(V_{t}^{1})^{\top}&\cdots&(V_{t}^{j})^{\top}&\cdots&(V_{t}^{L})^{\top}\\\
V_{t}^{1}&V_{t}^{1,1}&\cdots&(V_{t}^{j,1})^{\top}&\cdots&(V_{t}^{L,1})^{\top}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
V_{t}^{j}&V_{t}^{j,1}&\cdots&V_{t}^{j,k}&\cdots&(V_{t}^{L,k})^{\top}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
V_{t}^{L}&V_{t}^{L,1}&\cdots&V_{t}^{L,k}&\cdots&V_{t}^{L,L}\end{array}\right]$
$\displaystyle{\color[rgb]{0,0,0}\tilde{L}_{t}}$
$\displaystyle{\color[rgb]{0,0,0}=\left[\begin{array}[]{lllllll}(L_{t})^{\top}&(L^{1}_{t})^{\top}&\cdots&(L^{j}_{t})^{\top}&\cdots&(L^{L}_{t})^{\top}\end{array}\right]^{\top}.}$

With some abuse of notation: $V_{t}$ is also denoted by $V_{t}^{0,0}$;
$V_{t}^{j}$ is also denoted by $V_{t}^{j,0}$. Substituting the above
parametrizations in (36), we obtain Steps 5, 9. In regard to the initial
conditions: we only need $\hat{x}_{0}$ and $V_{0}$, we set $V_{0}^{j}=0$ and
$V_{0}^{j,k}>0$, with $j,k>1$, are set arbitrary such that $\tilde{V}_{0}>0$,
see (8). Using a parametrization for $\tilde{P}_{t|t}$ and $\tilde{P}_{t+1}$
as the one for $\tilde{V}_{t}$: from (37) and (38) we obtain Steps 6, 8, 11
and Steps 15, 17, 19, respectively. Regarding Step 14, recall that from Step 5
in Algorithm 1 we have

$\displaystyle\hat{\xi}_{t+1}=\tilde{A}\hat{\xi}_{t}+\tilde{G}_{t}(y_{t}-\tilde{C}\hat{\xi}_{t}),$
(39)

where $\tilde{G}_{t}=\tilde{A}\tilde{L}_{t}.$ Notice that $\hat{\xi}_{t}$ is
the predictor of $\xi_{t}$ given $Y_{t-1}$, and it can be partitioned as
$\hat{\xi}_{t}:=[\,(\hat{x}_{t|t-1}^{0})^{\top}\;(\hat{x}^{1}_{t|t-1})^{\top}\;\cdots\;(\hat{x}^{j}_{t|t-1})^{\top}\;\cdots\;(\hat{x}^{L}_{t|t-1})^{\top}\,]^{\top}$
where $\hat{x}^{j}_{t|t-1}=\hat{x}_{t-j|t-1},j\geq 0$. Substituting the
definitions of $\tilde{A}$, $\tilde{C}$, $\tilde{L}_{t}$ in (39) we obtain
Step 14. It remains to find an efficient way to find the risk sensitivity
parameter $\theta_{t}$. Indeed, in order to evaluate
$\gamma(\tilde{P}_{t+1},\cdot)$ for a specific $\theta$ in Step 8 of Algorithm
1, we have to perform the inversion and the eigenvalue decomposition of a
matrix whose dimension is proportional to $L$. The next result shows that it
is possible to find $\theta_{t}$ by considering a function which requires to
perform the inversion and the eigenvalue decomposition of matrices of
dimension $n\times n$, see Step 22. In this way, $\theta_{t}$ can be computed
in a numerically robust way in the case that $L$ is large.

###### Proposition 1

Consider $\gamma$ defined in (35). Then, $\theta_{t}$ is the unique solution
to $\gamma(\tilde{P}_{t+1},\theta_{t})=c_{t}$ if and only if $\theta_{t}$ is
the unique solution to $\bar{\gamma}(P^{L,L}_{t+1},\theta_{t})=c_{t}$ where

$\displaystyle\bar{\gamma}(P^{L,L}_{t+1},\theta_{t})=$
$\displaystyle-\frac{1}{2}\left\\{\operatorname{tr}\left[P_{t+1}^{L,L}(P^{L,L}_{t+1}-\theta_{t}I)^{-1}\right]\right.$
$\displaystyle\left.+\ln\operatorname{det}\left[I-P^{L,L}_{t+1}{(P^{L,L}_{t+1}-\theta_{t}I)^{-1}}\right]\right\\}=c_{t}.$

###### Proof 5

First, notice that $\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}=P_{t+1}^{L,L}$.
Thus, condition

$\theta_{t}<r(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top})^{-1}$

is equivalent to $\theta_{t}<r(P_{t+1}^{L,L})^{-1}$. By (35), we know

$\gamma(\tilde{P}_{t+1},\theta_{t})=\frac{1}{2}\left[\operatorname{tr}(K_{1})+\ln\operatorname{det}(K_{2})\right]$
(40)

where

$K_{1}=(I-\theta_{t}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1})^{-1}-I,\;\;K_{2}=I-\theta_{t}\tilde{H}^{\top}\tilde{H}\tilde{P}_{t+1}.$

Then,

$\displaystyle\operatorname{tr}(K_{1})=\operatorname{tr}\left\\{\left[(\tilde{P}^{-1}_{t+1}-\theta_{t}\tilde{H}^{\top}\tilde{H})\tilde{P}_{t+1}\right]^{-1}\right\\}-(L+1)n$
$\displaystyle=\operatorname{tr}\left[\tilde{P}^{-1}_{t+1}(\tilde{P}^{-1}_{t+1}-\theta_{t}\tilde{H}^{\top}\tilde{H})^{-1}\right]-(L+1)n$
$\displaystyle=\operatorname{tr}\left\\{\tilde{P}^{-1}_{t+1}\left[\tilde{P}_{t+1}-\tilde{P}_{t+1}\tilde{H}^{\top}\times(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}-\theta_{t}^{-1}I)^{-1}\tilde{H}\tilde{P}_{t+1}\right]\right\\}-(L+1)n$
$\displaystyle=\operatorname{tr}\left\\{\left[I-(\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}-\theta_{t}^{-1}I)^{-1}\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}\right]\right\\}-n$
$\displaystyle=-\operatorname{tr}(\Gamma_{t+1}^{-1}P^{L,L}_{t+1})$ (41)

where $\Gamma_{t+1}:=P_{t+1|t}^{L,L}-\theta_{t}^{-1}I_{n}$ and we exploited
the Woodbury matrix identity. Then,

$\displaystyle\ln\operatorname{det}(K_{2})$
$\displaystyle=\ln\operatorname{det}\left[(\tilde{P}^{-1}_{t+1}-\theta_{t}\tilde{H}^{\top}\tilde{H})\tilde{P}_{t+1}\right]$
$\displaystyle=-\ln\operatorname{det}\left[\tilde{P}^{-1}_{t+1}(\tilde{P}^{-1}_{t+1}-\theta_{t}\tilde{H}^{\top}\tilde{H})^{-1}\right]$
$\displaystyle=-\ln\operatorname{det}\left[\tilde{P}^{-1}_{t+1}(\tilde{P}_{t+1}-\tilde{P}_{t+1}\tilde{H}^{\top}\Gamma_{t+1}^{-1}\tilde{H}\tilde{P}_{t+1})\right]$
$\displaystyle=-\ln\operatorname{det}(I-\tilde{H}^{\top}\Gamma_{t+1}^{-1}\tilde{H}\tilde{P}_{t+1})$
$\displaystyle=-\ln\operatorname{det}(I-\tilde{H}\tilde{P}_{t+1}\tilde{H}^{\top}\Gamma_{t+1}^{-1})$
$\displaystyle=-\ln\operatorname{det}(I-P^{L,L}_{t+1}\Gamma_{t+1}^{-1})$ (42)

where we exploited the fact that
$\operatorname{det}(I_{m}+AB)=\operatorname{det}(I_{n}+BA)$ and
$A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{n\times m}$. Substituting (5)
and (5) in (40), we get the claim.

Using the Woodbury formula, we have

$\displaystyle\tilde{V}_{t+1}$
$\displaystyle=(\tilde{P}_{t+1}^{-1}-\theta_{t}\tilde{H}^{\top}\tilde{H})^{-1}$
$\displaystyle=\tilde{P}_{t+1}-\tilde{P}_{t+1}\tilde{H}^{\top}\Gamma_{t+1}^{-1}\tilde{H}^{\top}\tilde{P}_{t+1}$

and using the parametrization of $\tilde{V}_{t+1}$ and $\tilde{P}_{t+1}$ in
terms $V_{t+1}^{j,k}$ and $P_{t+1}^{j,k}$, we obtain Steps 23, 24, 26, 28.
Finally, Step 31 is derived as follows. The estimate of $\xi_{t+1}$ given
$Y_{t+1}$ is obtained by the update step of the standard Kalman filter:

$\displaystyle\hat{\xi}_{t+1|t+1}=\tilde{A}\hat{\xi}_{t|t}+\tilde{L}_{t+1}(y_{t+1}-\tilde{C}\tilde{A}\hat{\xi}_{t|t}).$
(43)

Notice that $\hat{\xi}_{t|t}$ can be partitioned as

$\hat{\xi}_{t|t}:=[\,(\hat{x}_{t|t}^{0})^{\top}\;(\hat{x}^{1}_{t|t})^{\top}\;\cdots\;(\hat{x}^{j}_{t|t})^{\top}\;\cdots\;(\hat{x}^{L}_{t|t})^{\top}\,]^{\top}$

where $\hat{x}^{j}_{t|t}=\hat{x}_{t-j|t},j\geq 0$. Substituting the
definitions of $\tilde{A}$, $\tilde{C}$, $\tilde{L}_{t}$ in (43), we obtain

$\displaystyle\hat{x}^{j}_{t+1|t+1}=\hat{x}^{j-1}_{t|t}+L^{j}_{t+1}(y_{t+1}-CA\hat{x}_{t|t}),\quad{j\geq
0}.$ (44)

Note that when $j=0$,
$\hat{x}^{j-1}_{t|t}=\hat{x}^{-1}_{t|t}=\hat{x}_{t+1|t}=A\hat{x}_{t|t}.$ In
addition, from Equation (29), it is not difficult to see (the derivation is
the same of the one in [43])

$L^{j}_{t+1}=V^{j}_{t+1}C^{\top}(CV_{t}C^{\top}+DD^{\top})^{-1}={V^{j}_{t+1}V^{-1}_{t+1}L_{t+1}},\quad
j\geq 0$

where $L^{0}_{t+1}=L_{t+1}$ and while $V^{0}_{t+1}=V_{t+1}$. Substituting the
latter in (44) we have

$\displaystyle\hat{x}^{j}_{t+1|t+1}=\hat{x}^{j-1}_{t|t}+V^{j}_{t+1}V_{t+1}^{-1}L_{t+1}(y_{t+1}-CA\hat{x}_{t|t}).$
(45)

from which Step 31 can be established.

Algorithm 2 Robust efficient version

1:$y_{0}\ldots y_{N}$, $\hat{x}_{0}$, $V_{0}$, $\color[rgb]{0,0,0}c_{t}$.

2:$\hat{x}_{t-L+1|t},t=L-1\ldots N$

3:Set $V_{0}^{j,k}=0$ for $j\neq k$, $V_{0}^{j,j}$ positive definite for $j>0$

4:for $t=0:N$ do

5: $L_{t}=V_{t}C^{\top}(CV_{t}C^{\top}+DD^{\top})^{-1}$

6: $P_{t|t}=V_{t}(I-C^{\top}L_{t}^{\top})$;

7: for $j=1:L$ do

8: $P^{j}_{t|t}=V^{j}_{t}(I-C^{\top}L_{t}^{\top})$;

9: $L^{j}_{t}=V^{j}_{t}C^{\top}(CV_{t}C^{\top}+R)^{-1}$;

10: for $k=1:L$ do

11: $P^{j,k}_{t|t}=V^{j,k}_{t}-V^{j}_{t}C^{\top}(L^{k}_{t})^{\top}$;

12: end for

13: end for

14: $\hat{x}_{t+1}=A\hat{x}_{t}+AL_{t}(y_{t}-C\hat{x}_{t})$;

15: $P_{t+1}=AP_{t|t}A^{\top}+BB^{\top}$;

16: for $j=1:L$ do

17: $P^{j}_{t+1}=P^{j-1}_{t|t}A^{\top}$;

18: for $k=1:L$ do

19: $P^{j,k}_{t+1}=P^{j-1,k-1}_{t|t}$;

20: end for

21: end for

22: Find $\theta_{t}$ s.t. ${\bar{\gamma}}(P^{L,L}_{t+1},\theta_{t})=c_{t}$;

23: $\Gamma_{t+1}=P^{L,L}_{t+1}-\theta^{-1}_{t}I_{n}$;

24: $V_{t+1}=P_{t+1}-(P^{L}_{t+1})^{\top}\Gamma_{t+1}^{-1}P^{L}_{t+1}$;

25: for $j=1:L$ do

26:
$V^{j}_{t+1}=P^{j}_{t+1}-(P^{L,j}_{t+1})^{\top}\Gamma_{t+1}^{-1}P^{L}_{t+1}$;

27: for $k=1:L$ do

28:
$V^{j,k}_{t+1}=P^{j,k}_{t+1}-(P^{L,j}_{t+1})^{\top}\Gamma_{t+1}^{-1}P^{L,k}_{t+1}$;

29: end for

30: end for

31:end for

32:for $t=L-1:N$ do

33:
$\displaystyle\hat{x}_{t-L+1|t}=\hat{x}_{t-L+1}+\sum_{j=t-L+1}^{t}V_{j}^{j+L-t-1}V_{j}^{-1}L_{j}\left(y_{j}-C\hat{x}_{j}\right)$

34:end for

Figure 1: Computational time of ARFLS (corresponding to the red part) and RFLS
(corresponding to the green part).

Computational complexity. We perform the asymptotic analysis of the
computational complexity, understood as the number of floating point
operations (flops) by using big O notation, of the efficient version of the
robust fixed-lag smoother (RFLS), i.e. Algorithm 2, versus the one of the
augmented robust fixed-lag smoother (ARFLS), i.e. Algorithm 1.

First, referring to ARFLS, Steps 2-5 have the same complexity of the standard
Kalman predictor, that is:

$O(n^{3}(L+1)^{3})+O(mn^{2}(L+1)^{2})+Q(m^{2}n(L+1))+O(m^{3}).$

Here, it is worth noting that $\tilde{Q}=\tilde{B}\tilde{B}^{\top}$ and
$\tilde{R}=\tilde{D}\tilde{D}^{\top}$ are computed offline. Then, in regard to
Step 6, the complexity to evaluate $\gamma(\tilde{P}_{t+1},\theta_{t})$ is
$O(n^{3}(L+1)^{3})$, see [46, Section 13.1 and 13.4]. Then, the computation of
$\theta_{t}\in(0,r(\tilde{P}_{t+1})^{-1})$ is accomplished by a bisection
method, see Algorithm 2 in [47]. Since at each step we spend constant time to
reduce the problem to an instance half its size [46, Section 4.10.2], the
complexity of Step 6 is

$O(n^{3}(L+1)^{3}\textrm{log}_{2}(r(\tilde{P}_{t+1})^{-1}/\varepsilon))$

where $\varepsilon>0$ is the selected accuracy, i.e. the solution found
satisfies the condition
$|\gamma(\tilde{P}_{t+1},\theta_{t})-c_{t}|\leq\varepsilon$. Step 7 has
complexity $O(n^{3}(L+1)^{3})$. Thus, the computational complexity of
Algorithm 1 is:

$\displaystyle O(n^{3}$
$\displaystyle(L+1)^{3})+O(mn^{2}(L+1)^{2})+Q(m^{2}n(L+1))+O(m^{3})$
$\displaystyle+O(n^{3}(L+1)^{3}\textrm{log}_{2}(r(\tilde{P}_{t+1})^{-1}/\varepsilon)).$

Accordingly, the complexity of Algorithm 1 with respect to the instance $L$ is
$O(L^{3})$.

Referring to RFLS, the total complexity of Steps 3-4, 12-13 and 21-22 is
$O(n^{3})+O(mn^{2})+Q(m^{2}n)+O(m^{3})$. Then, the complexity of Steps 5-11 is
$O(n^{3}L)+O(mn^{2}L^{2})+Q(m^{2}nL^{2})+O(Lm^{3})$; the complexity of Steps
14-19 is $O(n^{3}L)$; the complexity of Steps 23-28 is $O(n^{3}L^{2})$. Next,
Step 20 has the complexity of
$O(n^{3}\textrm{log}_{2}(r(P^{L,L}_{t+1})^{-1}/\varepsilon))$; it is worth
noting thay in Step 20 the computation of $\theta_{t}$ is done by using the
same bisection method of Step 6 in Algorithm 1: the difference is the
dimension of matrices $\tilde{P}_{t+1}$, of dimension $n(L+1)$, and
$P^{L,L}_{t+1}$, of dimension $n$. Finally, Step 31 has the complexity of
$O(mn^{2}L)+O(n^{3}L).$ Hence, the computational complexity of Algorithm 2 is:

$O(n^{3}L^{2})+O(mn^{2}L^{2})+Q(m^{2}nL^{2})+O(m^{3}L)+O(n^{3}\textrm{log}_{2}(r(P^{L,L}_{t+1})^{-1}/\varepsilon)).$

Thus, the complexity of Algorithm 2 with respect to the instance $L$ is
$O(L^{2})$. We conclude that Algorithm 2 is computationally more efficient
than Algorithm 1 and this advantage will become more pronounced as $L$ grows.

Finally, we also analyze the computational time with respect to the lag $L$
through a Monte Carlo study. In the latter, the lag ranges from $L=3$ up to
$L=50$. Each case is composed by 100 trials. In each trial, the matrices
$A,B,C,D$ of Model (1) with $n=2$ and $m=1$ are randomly generated as follows.
Each entry is drawn according to a uniform distribution in the interval
$[0,1]$. Then, matrix $A$ is rescaled in such a way that its maximum
eigenvalue (in modulus) is equal to 0.95. Then, an output sequence $Y_{N}$
with $N=500$ is generated. Fig. 1 shows the average value of the computational
time over 100 trials required by RFLS and ARFLS with $c_{t}=10^{-3}$ to
estimate the state trajectory from $Y_{N}$. The results were obtained using a
Huawei MateBook X Pro Laptop with Intel Core I5-8250U CPU and 8GB RAM. The
dashed lines defines the corresponding confidence intervals (with level 0.95).
It is possible to note that the computational time of these two algorithms
grows polynomially. As expected, the growth rate of RFLS is much smaller than
the one of ARFLS, i.e. RFLS drastically reduces the computational time.

##  5 Least-Favorable Model

In order to evaluate the performance of the robust fixed-lag smoother, we need
to construct its least favorable model solution to (11). The latter can be
characterized over a finite time interval $[0,T]$ by using arguments similar
to the ones in [27, Section V]. More precisely, the least favorable model
takes the form

$\displaystyle\eta_{t+1}$
$\displaystyle=\bar{A}_{t}\eta_{t}+\bar{B}_{t}\epsilon_{t}$ (46)
$\displaystyle y_{t}$
$\displaystyle=\bar{C}_{t}\eta_{t}+\bar{D}_{t}\epsilon_{t}$

where $\eta_{t}\triangleq[\,\xi_{t}^{\top}\;\tilde{e}_{t}^{\top}\,]^{\top}$.
Moreover,

$\displaystyle\begin{aligned}
&\bar{A}_{t}:=\left[\begin{array}[]{cc}\tilde{A}&\tilde{B}\tilde{S}_{t}\\\
0&\tilde{A}-\tilde{G}_{t}\tilde{C}+(\tilde{B}-\tilde{G}_{t}\tilde{D})\tilde{S}_{t}\end{array}\right]\\\
&\bar{B}_{t}:=\left[\begin{array}[]{c}\tilde{B}\\\
\tilde{B}-\tilde{G}_{t}\tilde{D}\end{array}\right]L_{t}\\\
&\bar{C}_{t}:=\left[\begin{array}[]{lll}\tilde{C}&\tilde{D}\tilde{S}_{t}\end{array}\right],\quad\bar{D}_{t}:=\tilde{D}L_{t}.\\\
\end{aligned}$ (47)

The matrices above are computed through the backward recursion illustrated in
Algorithm 3.

Algorithm 3 Backward recursion

1:$\tilde{G}_{0}\ldots\tilde{G}_{N}$, $\theta_{0}\ldots\theta_{N}$,
$\tilde{\Omega}_{N+1}^{-1}$

2:$\bar{A}_{t},\bar{B}_{t},\bar{C}_{t},\bar{D}_{t},t=0\ldots N$

3:$\tilde{\Omega}^{-1}_{N+1}=0$

4:for $t=N:0$ do

5:
$W^{-1}_{t+1}=\tilde{\Omega}_{t+1}^{-1}+\theta_{t}\tilde{H}^{\top}\tilde{H}$

6:
$\tilde{K}_{v_{t}}=(I-(\tilde{B}-\tilde{G}_{t}\tilde{D})^{\top}W_{t+1}^{-1}(\tilde{B}-\tilde{G}_{t}\tilde{D}))^{-1}$

7:
$\tilde{S}_{t}=\tilde{K}_{v_{t}}(\tilde{B}-\tilde{G}_{t}\tilde{D})^{\top}W_{t+1}^{-1}(\tilde{A}-\tilde{G}_{t}\tilde{C})$

8: Compute $L_{t}$ such that $\tilde{K}_{v_{t}}=L_{t}L_{t}^{\top}$

9: Compute $\bar{A}_{t}$, $\bar{B}_{t}$, $\bar{C}_{t}$, $\bar{D}_{t}$ as in
(47)

10:
$\tilde{\Omega}_{t}^{-1}=(\tilde{A}-\tilde{G}_{t}\tilde{C})^{\top}W_{t+1}^{-1}(\tilde{A}-\tilde{G}_{t}\tilde{C})+\tilde{S}_{t}^{\top}\tilde{K}_{v_{t}}^{-1}\tilde{S}_{t}$

11:end for

It is worth noting that (46) is the least favorable model corresponding to the
augmented state $\xi_{t}$. It is then natural to wonder whether such a least
favorable model reduces to a least favorable model corresponding to the state
$x_{t}$. The answer is affirmative. This justifies why in the minimax problem
(11) we did not need to impose that $\tilde{\phi}_{t}(z_{t}|x_{t})$ preserve
the same structure of the augmented state space model in (2). Substituting
$\bar{A}_{t}$, $\bar{B}_{t}$, $\bar{C}_{t}$, $\bar{D}_{t}$ in (46), we obtain

$\displaystyle\xi_{t+1}$
$\displaystyle=\tilde{A}\xi_{t}+\tilde{B}(\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t})$
$\displaystyle y_{t}$
$\displaystyle=\tilde{C}\xi_{t}+\tilde{D}(\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t})$

which is the augmented state space of the least favorable model

$\displaystyle x_{t+1}$
$\displaystyle={A}x_{t}+{B}(\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t})$
(48) $\displaystyle y_{t}$
$\displaystyle={C}x_{t}+{{D}}(\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t}).$

Finally, consider a fixed-lag smoother of the form

$\displaystyle\hat{x}_{t-L+1|t}^{\prime}=\hat{x}_{t-L+1}^{\prime}+\sum_{j=t-L+1}^{t}G_{j}^{\prime\,L-t+j}(y_{j}-C\hat{x}_{j}^{\prime})$

where $G_{t}^{\prime\,j}$, with $0\leq j\leq L$, are the arbitrary gains of a
fixed-lag smoothing algorithm. It is not difficult to see that we can rewrite
such a smoother as
$\hat{x}_{t-L+1|t}^{\prime}=\tilde{H}\hat{\xi}_{t+1}^{\prime}$ where

$\displaystyle\hat{\xi}_{t+1}^{\prime}$
$\displaystyle=\tilde{A}\hat{\xi}_{t}^{\prime}+\tilde{G}_{t}^{\prime}(y_{t}-\tilde{C}\hat{\xi}_{t}^{\prime})$
$\displaystyle\tilde{G}_{t}^{\prime}$
$\displaystyle=\left[\begin{array}[]{ccc}(G_{t}^{\prime\,0})^{\top}&\ldots&(G_{t}^{\prime\,L})^{\top}\end{array}\right]^{\top}.$

To evaluate its performance under the least favorable model in (48), we define
the corresponding smoothing error
$e_{t-L+1}^{\prime}=x_{t-L+1}-\hat{x}_{t-L+1|t}^{\prime}$. Then it is not
difficult to see that $e_{t-L+1}^{\prime}$ is zero mean and with covariance
matrix

$\bar{\Pi}_{t+1}=\tilde{H}\tilde{\Pi}_{t+1}\tilde{H}^{\top}$

where: $\tilde{\Pi}_{t+1}$ is the $(L+1)n\times(L+1)n$ submatrix of
$\Pi_{t+1}$ in position $(1,1)$; $\Pi_{t+1}$ is the solution to the Lyapunov
equation

$\displaystyle\Pi_{t+1}$
$\displaystyle=\left(\bar{A}_{t}-\left[\begin{array}[]{c}\tilde{G}_{t}^{\prime}\\\
0\end{array}\right]\bar{C}_{t}\right)\Pi_{t}\left(\bar{A}_{t}-\left[\begin{array}[]{c}\tilde{G}_{t}^{\prime}\\\
0\end{array}\right]\bar{C}_{t}\right)^{\top}$ $\displaystyle\hskip
14.22636pt+\left(\bar{B}_{t}-\left[\begin{array}[]{c}\tilde{G}_{t}^{\prime}\\\
0\end{array}\right]\bar{D}_{t}\right)\left(\bar{B}_{t}-\left[\begin{array}[]{c}\tilde{G}_{t}^{\prime}\\\
0\end{array}\right]\bar{D}_{t}\right)^{\top}$

with initial condition $\Pi_{0}=\mathbbm{1}_{2}\otimes\tilde{V}_{0}$ and
$\mathbbm{1}_{2}$ is the $2\times 2$ matrix whose entries are equal to one.

##  6 Robust Target Tracking

We compare the performance of the robust and standard fixed-lag smoothers in a
maneuvering target tracking problem where model uncertainty is present. More
precisely, we consider as nominal model the second-order Singer model with an
exponentially autocorrelated noise, see [31, 48, 49]. The state vector is
defined as $x:=\left[p^{lat},v^{lat},p^{lon},v^{lon}\right]^{\top}$ where
$p^{lat}$ and $v^{lat}$ denote the target position and velocity along the
latitudinal direction, respectively; $p^{lon}$ and $v^{lon}$ denote the target
position and velocity along the longitudinal direction, respectively. This
model can be written as (1) with

$A=\left[\begin{array}[]{cccc}1&T&0&0\\\ 0&1&0&0\\\ 0&0&1&T\\\
0&0&0&1\end{array}\right]$

where $T=0.01$ denotes the sampling period; $B$ is such that $BB^{\top}=Q$ and

$Q=2\alpha\sigma_{m}^{2}\left[\begin{array}[]{cccc}T^{3}/3&T^{2}/2&0&0\\\
T^{2}/2&T&0&0\\\ 0&0&T^{3}/3&T^{2}/2\\\ 0&0&T^{2}/2&T\end{array}\right].$

where we assume that $\sigma_{m}^{2}=5$ and $1/\alpha=0.5$ which are the
instantaneous variance of the velocity and the time constant of the target
velocity autocorrelation, respectively. Moreover,

$C=\left[\begin{array}[]{cccc}1&0&0&0\\\ 0&0&1&0\end{array}\right],\quad
DD^{\top}=I_{2}$

and the output $y_{t}$ denotes the noisy position measurements along the two
directions. Finally, $x_{0}$ is Gaussian distributed with zero mean and
covariance matrix $V_{0}=\operatorname{diag}(50,5,50,5)$.

In practice, the nominal model above does not coincide with the actual one
(e.g. the nominal parameters $\alpha$ and $\sigma_{m}^{2}$ are typically
imprecise). In what follows, we assume that the actual model belongs to the
ambiguity set $\mathcal{B}_{t}$. More precisely, we consider two scenarios:
the first one considers the ambiguity set with $c_{t}=10^{-3}$, while the
second one with $c_{t}=5\cdot 10^{-3}$, i.e. the latter is larger than the
former. We compare the proposed robust fixed-lag smoother, denoted by RFLS,
and the standard fixed-lag smoother, denoted by FLS, both with $L=20$.

Fig. 2(a) shows the variances of the smoothing error under the least favorable
model in the first ambiguity set with $c_{t}=10^{-3}$, while Fig. 2(b) shows
the variances of the smoothing error under the least favorable model in the
second ambiguity set with $c_{t}=5\cdot 10^{-3}$. It is possible to see that
RFLS outperforms FLS. It is also worth noting that the higher $c_{t}$ is, the
more RFLS outperforms FLS.

(a) Ambiguity set with $c=10^{-3}$

(b) Ambiguity set with $c=5\cdot 10^{-3}$

Figure 2:  Scalar variances of the smoothing error for RFLS and FLS under the
least favorable models in the different ambiguity sets.

(a) Ambiguity set with $c=10^{-3}$

(b) Ambiguity set with $c=5\cdot 10^{-3}$

Figure 3:  Trajectories generated by the least favorable models in the
different ambiguity sets (black line) and the estimated trajectories with RFLS
(blue line) and FLS (red line). The cross denotes the endpoint of the target
trajectory.

In what follows, we consider a sample $Y_{N}=\\{y_{0},y_{1}\ldots y_{N}\\}$
with $N=500$ generated by the least favorable model (48) in the first
ambiguity set with $c_{t}=10^{-3}$ and $\xi_{0}=0$.

| RFLS | FLS | RFLS | FLS
---|---|---|---|---
| $c_{t}=10^{-3}$ | $c_{t}=10^{-3}$ | $c_{t}=5\cdot 10^{-3}$ | $c_{t}=5\cdot 10^{-3}$
$\operatorname{RMSE}_{lat}$ | 0.2918 | 0.4093 | 0.4962 | 0.6112
$\operatorname{RMSE}_{lon}$ | 0.2891 | 0.3280 | 0.4804 | 0.6115
Table 1: RMSE along the latitudinal direction and longitudinal direction.

Fig. 3(a) shows the target trajectory (black line) and the ones estimated with
RFLS (blue line) and FLS (red line). We also perform the same experiment using
the least favorable model in the second ambiguity set with $c_{t}=5\times
10^{-3}$. Moreover, we quantify the performance of the smoothers in these two
experiments through the root mean-square error (RMSE) along the latitudinal
direction and longitudinal direction

$\operatorname{RMSE}_{lat}=\sqrt{\frac{1}{N}\sum_{t=1}^{N}\left(p_{t}^{lat}-\hat{p}_{t}^{lat}\right)^{2}},$
$\operatorname{RMSE}_{lon}=\sqrt{\frac{1}{N}\sum_{t=1}^{N}\left(p_{t}^{lon}-\hat{p}_{t}^{lon}\right)^{2}}$

whose values are displayed in Table 1. As expected, RFLS better reduces the
influence of the modeling error on the estimation accuracy than others.

##  7 Robust parameter estimation

Consider the following state space model

$\displaystyle\tilde{\mathcal{M}}(\alpha):x_{t+1}$
$\displaystyle={A}(\alpha)x_{t}+{B}(\alpha)\tilde{\upsilon}_{t}$ (49)
$\displaystyle y_{t}$
$\displaystyle={C}(\alpha)x_{t}+{{D}(\alpha)}\tilde{\upsilon}_{t}$

where $x_{0}\sim\mathcal{N}(\hat{x}_{0},V_{0})$,
$\tilde{\upsilon}_{t}\sim\mathcal{N}(0,\tilde{R}_{t})$, i.e. $\tilde{v}_{t}$
is a nonstationary process, and the matrices ${A}(\alpha)$, ${B}(\alpha)$,
${C}(\alpha)$ and ${D}(\alpha)$ are parameterized by $\alpha$. In many
practical applications, $\alpha$ is not known and needs to be estimated from
the observed data $Y_{N}=\\{y_{0},y_{1}\ldots y_{N}\\}$. In plain words, the
latter is a system identification problem where the model class is
$\tilde{\mathcal{M}}=\\{\tilde{\mathcal{M}}(\alpha),\alpha\in\Theta\\}$ and
$\Theta$ is the parameter space. A well established paradigm to find $\alpha$
is the maximum likelihood (ML) principle. However, it is usually difficult to
find an explicit expression of the ML function under Model (49). Such a
problem is typically addressed by using the expectation-maximization (EM)
algorithm, see [6, Algorithm 12.3], which computes a lower bound of the ML
function through the iterative scheme:

  * 1.

Set up an initial guess $\alpha^{0}$;

  * 2.

For $n=0,1,...$:

    1. $\bullet$

E-step: compute

$\displaystyle\mathcal{Q}(\alpha,\alpha^{n}):=\int\tilde{f}_{\alpha^{n}}(X_{N+1}|Y_{N})\log\tilde{f}_{\alpha}(Y_{N},X_{N+1})dX_{N+1}$

    2. $\bullet$

M-step: compute
$\alpha^{(n+1)}=\arg\max_{\alpha}\mathcal{Q}(\alpha,\alpha^{n})$

where $\tilde{f}_{\alpha}(Y_{N},X_{N+1})$ is the joint density of
$Y_{N}=\\{y_{0},y_{1}\ldots,y_{N}\\}$ and
$X_{N+1}=\\{x_{0},x_{1},\ldots,x_{N+1}\\}$ under
$\tilde{\mathcal{M}}(\alpha)$. In most cases, however, the covariance matrix
$\tilde{R}_{t}$ of the noise process $\tilde{v}_{t}$ is not known and it is
also time-varying. Such a matrix is typically designed empirically. However,
this would require to have the possibility to make more experiments on the
system, i.e. a requirement that is not always met. Alternatively, we can
select a nominal covariance matrix for the noise process using some a priori
knowledge. However, this causes the nominal model to be inconsistent with the
actual one and thus the reliability of the estimate of $\alpha$ will be
compromised.

A possible way to address this model uncertainty is to understand (49) as the
least favorable model (48) where $\tilde{\upsilon}_{t}$ is equal to
$\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t}$. More precisely, assume that we
want to estimate $\alpha$ only knowing the nominal state space model
$\mathcal{M}$, i.e.

$\displaystyle\mathcal{M}(\alpha):x_{t+1}$
$\displaystyle={A}(\alpha)x_{t}+{B}(\alpha)\upsilon_{t}$ (50) $\displaystyle
y_{t}$ $\displaystyle={C}(\alpha)x_{t}+{{D}(\alpha)}\upsilon_{t}$

where $v_{t}$ is normalized WGN, in particular it is a stationary process. It
is worth noting that the least favorable model in (49) does not belong to
$\mathcal{M}$. Notice that, it is not restrictive to assume that the
covariance matrix of $v_{t}$ is equal to the identity. Indeed, in the case its
covariance matrix is $R$, then we can always take
$\tilde{B}(\alpha)=B(\alpha)R^{1/2}$ and $\tilde{D}(\alpha)=D(\alpha)R^{1/2}$.
Then, the least favorable model, solution to (11), is (49) with
$\tilde{v}_{t}=\tilde{S}_{t}\tilde{e}_{t}+L_{t}\epsilon_{t}$. Hence,

$\tilde{R}_{t}=\tilde{S}_{t}\left[\begin{array}[]{cc}0&I_{n}\end{array}\right]\Pi_{t}\left[\begin{array}[]{c}0\\\
I_{n}\end{array}\right]\tilde{S}_{t}^{T}+L_{t}L_{t}^{T}.$

At this point we can use the density, say
$\tilde{f}^{0}_{\alpha^{n}}(X_{N+1},Y_{N})$, of the least favorable model in
order to compute the lower bound $\mathcal{Q}(\alpha,\alpha^{n})$. On the
other hand, an approximation of the moments of $\tilde{f}_{\alpha^{n}}^{0}$
required in $\mathcal{Q}(\alpha,\alpha^{n})$ can be constructed by the robust
fixed-lag smoother of Algorithms 2 and 3 leading to

$\displaystyle\mathcal{Q}$
$\displaystyle(\alpha,\alpha^{n})\propto-\frac{1}{2}\sum_{t=0}^{N-1}\log|2\pi
B(\alpha)\tilde{R}_{t}B^{T}(\alpha)|$
$\displaystyle-\frac{1}{2}\operatorname{tr}\left\\{\sum_{t=0}^{N-1}\left[\left(B(\alpha)\tilde{R}_{t}B^{T}(\alpha)\right)^{-1}\left(\Phi_{1,t}\right.\right.\right.$
$\displaystyle\left.\left.\left.-\Phi_{2,t}A^{T}(\alpha)-A(\alpha)\Phi_{2,t}^{T}+A(\alpha)\Phi_{3,t}A^{T}(\alpha)\right)\right]\right\\}$
$\displaystyle-\frac{1}{2}\sum_{t=0}^{N-1}\log|2\pi
D(\alpha)\tilde{R}_{t+1}D^{T}(\alpha)|$
$\displaystyle-\frac{1}{2}\operatorname{tr}\left\\{\sum_{t=0}^{N-1}\left[\left(D(\alpha)\tilde{R}_{t+1}D^{T}(\alpha)\right)^{-1}\left(\Phi_{4,t}\right.\right.\right.$
$\displaystyle\left.\left.\left.-C(\alpha)\Phi_{5,t}^{T}-\Phi_{5,t}C^{T}(\alpha)+C(\alpha)\Phi_{1,t}C^{T}(\alpha)\right)\right]\right\\}$

where

$\displaystyle\Phi_{1,t}$
$\displaystyle={\bar{\Pi}^{L,L}}_{t+L-1}+\hat{x}_{t+1|t+L-1}\hat{x}^{T}_{t+1|t+L-1},$
$\displaystyle\Phi_{2,t}$
$\displaystyle=({\bar{\Pi}}^{L,L-1}_{t+L-1})^{T}+\hat{x}_{t+1|t+L-1}\hat{x}^{T}_{t|t+L-1},$
$\displaystyle\Phi_{3,t}$
$\displaystyle={\bar{\Pi}}^{L,L}_{t+L-1}+\hat{x}_{t|t+L-1}\hat{x}^{T}_{t|t+L-1},$
$\displaystyle\Phi_{4,t}$
$\displaystyle=y_{t+1}y^{T}_{t+1},~{}~{}\Phi_{5,t}=y_{t+1}\hat{x}^{T}_{t+1|t+L-1},$

and ${\bar{\Pi}^{j,k}}_{t}$ is the $n\times n$ block in position $(j,k)$ of
$\tilde{\Pi}_{t}$. Then, $\hat{x}_{t|t+L-1}$, $\hat{x}_{t-1|t+L-1}$ and
$y_{t}$ are given by Algorithm 2. Clearly, such an approximation is legitimate
if the lag $L$ is chosen big enough.

Figure 4: Estimation error for EM, EM-FL and REM in the Monte Carlo
experiment.

Next, we show a numerical example. We consider the problem to estimate the
parameter $\alpha=[\,\alpha_{1}\;\alpha_{2}\,]$ using the nominal model class
(50) with

$\displaystyle A(\alpha)$
$\displaystyle=\left[\begin{array}[]{cc}\alpha_{1}&1\\\
0&\alpha_{2}\end{array}\right],~{}~{}~{}~{}~{}B=\left[\begin{array}[]{ccc}0.01&0&0\\\
0&0.45&0\end{array}\right]$ $\displaystyle C$
$\displaystyle=\left[\begin{array}[]{ll}1&-1\end{array}\right],~{}~{}~{}~{}~{}~{}D=\left[\begin{array}[]{lll}0&0&0.01\end{array}\right],$

and the collected data $Y_{N}=\\{y_{0},y_{1}\ldots y_{N}\\}$. We assume
$x_{0}\sim\mathcal{N}\left(0,V_{0}\right)$ with $V_{0}=0.0001I_{2}$. We assume
that the actual model has the same structure of (50) with
$\alpha^{\star}=[\,0.1\;0.9\,]$, but the actual noise process, say
$\tilde{v}_{t}$, is not stationary and not known. We solve the aforementioned
system identification problem by means of the REM method introduced in Section
7 with $L=50$. Moreover, we compare it with: the standard EM method where the
state estimation task is performed by the RTS smoother; the “fixed-lag” EM
(EM-FL) method where the state estimation task is performed by the standard
fixed-lag smoother. To estimate the effectiveness of the REM method we assume
that the actual model is the least favorable one belonging to the ambiguity
set with $c=2\cdot 10^{-2}$. Moreover, we consider a Monte Carlo experiment
with $100$ trials. More precisely, in each trial, we generate the data set
$Y_{N}$ according to (48) with $\alpha=\alpha^{\star}$ and $N=1000$. The
initial parameter estimates $\alpha_{1}^{0}$ and $\alpha_{2}^{0}$, with
$\alpha^{0}:=[\,\alpha_{1}^{0}\;\alpha_{2}^{0}\,]$, are drawn from a uniform
distribution with interval $[0.4,0.9]$ and $[0.07,0.13]$, respectively. Then,
the termination condition is
$\left\|\alpha^{n+1}-\alpha^{n}\right\|\leq\epsilon$ where $\epsilon=10^{-3}$.

Then, in order to compare the performance of these algorithms, we consider the
estimation error $\|\hat{\alpha}-\alpha^{\star}\|$ where $\hat{\alpha}$ is the
parameter estimate obtained at the last stage by EM/EM-FL and REM. Fig. 4
shows the boxplot of the estimation error for EM, EM-FL and REM. We see that
REM outperforms EM and EM-FL. It is worth noting that EM and EM-FL perform in
the same way. This means that the value of the lag $L$ has been chosen large
enough in EM-FL and thus the fixed-lag smoother represents a good
approximation of the RTS smoother.

(a) First state component

(b) Second state component

Figure 5:  Comparison among the actual state components and the estimated ones
obtained by EM and REM over the interval $[600,800]$.

Finally, Fig. 5 shows the estimated state trajectory in the last stage by EM
and REM, with $t\in[600,800]$, in a trial of the Monte Carlo experiment. As we
can see, the one obtained with REM, and thus using RFLS, is slightly better
than the one with EM, and thus using RTS. Although this advantage is not
prominent, it made a dramatic improvement in the performance of the parameter
estimator.

##  8 Conclusion

In this paper, we have proposed a robust fixed-lag smoother in the case that
the actual model is different from the nominal one. More precisely, this
paradigm solves a minimax game with two players: one selects the least
favorable model in a prescribed ambiguity set, the other designs the optimal
estimator based on the least favorable model. We also proposed an efficient
implementation of the robust fixed-lag smoother in order to reduce the
computational burden and avoid numerical instabilities. Then, we have
characterized the least favorable model for the robust smoother over a finite
time horizon. Finally, we have presented some numerical examples showing the
effectiveness of the proposed robust fixed-lag smoother.

## References

  * [1] F. S. Cattivelli, A. H. Sayed, Diffusion strategies for distributed Kalman filtering and smoothing, IEEE Transactions on Automatic Control 55 (9) (2010) 2069–2084.
  * [2] T. C. Dong-Si, A. I. Mourikis, Motion tracking with fixed-lag smoothing: Algorithm and consistency analysis, in: IEEE International Conference on Robotics and Automation, 2011, pp. 5655–5662.
  * [3] F. Papi, M. Bocquel, M. Podt, Y. Boers, Fixed-lag smoothing for Bayes optimal knowledge exploitation in target tracking, IEEE Transactions on Signal Processing 62 (12) (2014) 3143–3152.
  * [4] R. H. Shumway, D. S. Stoffer, An approach to time series smoothing and forecasting using the EM algorithm, Journal of time series analysis 3 (4) (1982) 253–264.
  * [5] A. Wills, T. B. Schön, B. Ninness, Estimating state-space models in innovations form using the expectation maximisation algorithm, in: IEEE Conference on Decision and Control (CDC), 2010, pp. 5524–5529.
  * [6] S. Särkkä, Bayesian filtering and smoothing, Cambridge University Press, 2013.
  * [7] T. Kailath, A. H. Sayed, B. Hassibi, Linear estimation, Prentice Hall, 2000.
  * [8] P. Whittle, Risk-sensitive Optimal Control, J. Wiley, Chichester, England, 1980\.
  * [9] R. Banavar, J. Speyer, Properties of risk-sensitive filters/estimators, IEE Proceedings-Control Theory and Applications 145 (1) (1998) 106–112.
  * [10] B. Hassibi, A. Sayed, T. Kailath, Indefinite-Quadratic Estimation and Control- A Unified Approach to $H^{2}$ and $H^{\infty}$ Theories, Society for Industrial and Applied Mathematics, Philadelphia, 1999.
  * [11] J. Huang, D. Shi, T. Chen, Distributed robust state estimation for sensor networks: A risk-sensitive approach, in: IEEE Conference on Decision and Control (CDC), 2018, pp. 6378–6383.
  * [12] B. C. Levy, M. Zorzi, A contraction analysis of the convergence of risk-sensitive filters, SIAM Journal on Control and Optimization 54 (4) (2016) 2154–2173.
  * [13] J. L. Speyer, C. H. Fan, R. N. Banavar, Optimal stochastic estimation with exponential cost criteria, in: IEEE Conference on Decision and Control (CDC), 1992, pp. 2293–2298.
  * [14] M.-G. Yoon, V. A. Ugrinovskii, I. R. Petersen, Robust finite horizon minimax filtering for discrete-time stochastic uncertain systems, Systems & control letters 52 (2) (2004) 99–112.
  * [15] R. K. Boel, M. R. James, I. R. Petersen, Robustness and risk-sensitive filtering, IEEE Transactions on Automatic Control 47 (3) (2002) 451–461.
  * [16] L. P. Hansen, T. J. Sargent, Robust estimation and control under commitment, Journal of Economic Theory 124 (2) (2005) 258–301.
  * [17] L. Hansen, T. J. Sargent, Recursive robust estimation and control without commitment, Journal of Economic Theory 136 (1) (2007) 1–27.
  * [18] M. Zorzi, B. C. Levy, Robust kalman filtering: Asymptotic analysis of the least favorable model, in: 2018 IEEE Conference on Decision and Control (CDC), IEEE, 2018, pp. 7124–7129.
  * [19] B. C. Levy, R. Nikoukhah, Robust state-space filtering under incremental model perturbations subject to a relative entropy tolerance, IEEE Transactions on Automatic Control 58 (2013) 682–695.
  * [20] S. Abadeh, V. Nguyen, D. Kuhn, P. Esfahani, Wasserstein distributionally robust Kalman filtering, in: Advances in Neural Information Processing Systems, 2018, pp. 8474–8483.
  * [21] M. Zorzi, Robust Kalman filtering under model perturbations, IEEE Transactions on Automatic Control 62 (6) (2016) 2902–2907.
  * [22] M. Zorzi, Convergence analysis of a family of robust Kalman filters based on the contraction principle, SIAM Journal on Control and Optimization 55 (5) (2017) 3116–3131.
  * [23] M. Zorzi, On the robustness of the bayes and wiener estimators under model uncertainty, Automatica 83 (2017) 133–140.
  * [24] M. Zorzi, Distributed Kalman filtering under model uncertainty, IEEE Transactions on Control of Network Systems 7 (2) (2019) 990–1001.
  * [25] A. Emanuele, F. Gasparotto, G. Guerra, M. Zorzi, Robust distributed Kalman filtering: On the choice of the local tolerance, Sensors 20 (11) (2020) 3244.
  * [26] A. Zenere, M. Zorzi, On the coupling of model predictive control and robust Kalman filtering, IET Control Theory and Applications 12 (13) (2018) 1873–1881.
  * [27] S. Yi, M. Zorzi, Robust Kalman filtering under model uncertainty: The case of degenerate densities, IEEE Transactions on Automatic Control 67 (7) (2022) 3458–3471.
  * [28] S. Yi, M. Zorzi, Low-rank Kalman filtering under model uncertainty, in: IEEE Conference on Decision and Control (CDC), 2020, pp. 2930–2935.
  * [29] R. Piché, S. Särkkä, J. Hartikainen, 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, 2012, pp. 1–6.
  * [30] A. Aravkin, B. Bell, J. Burke, G. Pillonetto, An $\ell_{1}$-laplace robust Kalman smoother, IEEE Transactions on Automatic Control 56 (12) (2011) 2898–2911.
  * [31] S. Farahmand, G. B. Giannakis, D. Angelosante, Doubly robust smoothing of dynamical processes via outlier sparsity constraints, IEEE Transactions on Signal Processing 59 (10) (2011) 4529–4543.
  * [32] Y. Huang, Y. Zhang, N. Li, J. Chambers, A robust gaussian approximate fixed-interval smoother for nonlinear systems with heavy-tailed process and measurement noises, IEEE Signal Processing Letters 23 (4) (2016) 468–472.
  * [33] J. Navarro-Moreno, R. M. Fernández-Alcalá, J. D. Jiménez-López, J. C. Ruiz-Molina, Widely linear estimation for multisensor quaternion systems with mixed uncertainties in the observations, Journal of the Franklin Institute 356 (5) (2019) 3115–3138.
  * [34] R. M. Fernández-Alcalá, J. Navarro-Moreno, J. D. Jiménez-López, J. C. Ruiz-Molina, Semi-widely linear estimation algorithms of quaternion signals with missing observations and correlated noises, Journal of the Franklin Institute 357 (5) (2020) 3075–3096.
  * [35] S. Zhao, J. Wang, Y. S. Shmaliy, F. Liu, Discrete time $q$-lag maximum likelihood FIR smoothing and iterative recursive algorithm, IEEE Transactions on Signal Processing 69 (2021) 6342–6354.
  * [36] Y. Huang, Y. Zhang, Y. Zhao, L. Mihaylova, J. A. Chambers, Robust Rauch–Tung–Striebel smoothing framework for heavy-tailed and/or skew noises, IEEE Transactions on Aerospace and Electronic Systems 56 (1) (2019) 415–441.
  * [37] S. Zhao, B. Huang, Y. S. Shmaliy, Bayesian state estimation on finite horizons: The case of linear state–space model, Automatica 85 (2017) 91–99.
  * [38] Y. Huang, F. Zhu, G. Jia, Y. Zhang, A slide window variational adaptive Kalman filter, IEEE Transactions on Circuits and Systems II: Express Briefs 67 (12) (2020) 3552–3556.
  * [39] B. Kwon, Z. Quan, S. Han, A robust fixed-lag receding horizon smoother for uncertain state space models, International Journal of Adaptive Control and Signal Processing 29 (11) (2015) 1354–1366.
  * [40] O. U. Rehman, I. R. Petersen, A robust continuous-time fixed-lag smoother for nonlinear uncertain systems, International Journal of Robust and Nonlinear Control 26 (2) (2016) 345–364.
  * [41] U. Theodor, U. Shaked, C. E. de Souza, A game theory approach to robust discrete-time H∞-estimation, IEEE Transactions on Signal Processing 42 (6) (1994) 1486–1495.
  * [42] P. Bolzern, P. Colaneri, G. De Nicolao, On discrete-time H∞ fixed-lag smoothing, IEEE transactions on signal processing 52 (1) (2004) 132–141.
  * [43] J. Moore, Discrete-time fixed-lag smoothing algorithms, Automatica 9 (2) (1973) 163–173.
  * [44] L. Hansen, T. Sargent, Robustness, Princeton University Press, Princeton, NJ, 2008\.
  * [45] B. C. Levy, R. Nikoukhah, Robust least-squares estimation with a relative entropy constraint, IEEE Transactions on Information Theory 50 (1) (2004) 89–104.
  * [46] S. S. Skiena, The algorithm design manual, Vol. 2, Springer, 1998.
  * [47] A. Zenere, M. Zorzi, On the coupling of model predictive control and robust Kalman filtering, IET Control Theory & Applications 12 (13) (2018) 1873–1881.
  * [48] R. A. Singer, Estimating optimal tracking filter performance for manned maneuvering targets, IEEE Transactions on Aerospace and Electronic Systems AES-6 (4) (1970) 473–483.  doi:10.1109/TAES.1970.310128.
  * [49] S. Yi, X. Jin, T. Su, Z. Tang, F. Wang, N. Xiang, J. Kong, Online denoising based on the second-order adaptive statistics model, Sensors 17 (7) (2017) 1668\.