m #1 θ μ μ_p μ_q μ^KF ¯x

# Variational Kalman Filtering with H∞-Based Correction for Robust Bayesian
Learning in High Dimensions

Niladri Das, Jed A. Duersch, and Thomas A. Catanach  Niladri Das
(corresponding author, ndas@sandia.gov), Jed Duersch (jaduers@sandia.gov), and
Thomas A. Catanach<EMAIL_ADDRESS>are with Sandia National Laboratories,
Livermore, CA 94550, USA. Sandia National Laboratories is a multimission
laboratory managed and operated by National Technology & Engineering Solutions
of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for
the U.S. Department of Energy’s National Nuclear Security Administration under
contract DE-NA0003525. This paper describes objective technical results and
analysis. Any subjective views or opinions that might be expressed in the
paper do not necessarily represent the views of the U.S. Department of Energy
or the United States Government.

###### Abstract

In this paper, we address the problem of convergence of sequential variational
inference filter (VIF) through the application of a robust variational
objective and H∞-norm based correction for a linear Gaussian system. As the
dimension of state or parameter space grows, performing the full Kalman update
with the dense covariance matrix for a large scale system requires increased
storage and computational complexity, making it impractical. The VIF approach,
based on mean-field Gaussian variational inference, reduces this burden
through the variational approximation to the covariance usually in the form of
a diagonal covariance approximation. The challenge is to retain convergence
and correct for biases introduced by the sequential VIF steps. We desire a
framework that improves feasibility while still maintaining reasonable
proximity to the optimal Kalman filter as data is assimilated. To accomplish
this goal, a H∞-norm based optimization perturbs the VIF covariance matrix to
improve robustness. This yields a novel VIF-H∞ recursion that employs
consecutive variational inference and H∞ based optimization steps. We explore
the development of this method and investigate a numerical example to
illustrate the effectiveness of the proposed filter.

##  I INTRODUCTION

The sequential estimation of a system’s state is essential for many problems
in science and engineering. Filtering allow sequential observations to be
integrated into a time-varying estimate of a system’s state. A large number of
such systems can be modelled as a linear Gaussian system. The Kalman filter
(KF), derived from the optimal Bayesian filter, is used for state/parameter
estimation of such problems. As the dimension of the system becomes very
large, the standard formulations of KF become computationally intractable due
to matrix storage and computation requirements. We anticipate this task
becoming increasingly problematic as traditional data assimilation methods,
like Kalman filters, are used to sequentially update very high-dimensional
machine learning models, as is the case in learning for deep neural networks
[8]. Therefore, the main challenge is to develop an approximating algorithm
that is robust enough to asymptotically converge to the data generating
process.

There are variants of KF that are computationally efficient for high
dimensional problems. The reduced rank Kalman Filter [5, 6, 11] searches for
lower dimensional subspace to operate in. However, their performance assumes
there exists a low enough dimensional subspace capturing the bulk of the
system’s structure while still computationally tractable. This often does not
hold in very high dimensions. Further, these methods suffer from the lack of
structured mechanism to find the effective lower dimensional projection
operator. Added to that, such operators are typically fixed in time, unable to
switch to a different optimal lower dimensional space. Direct approximation
techniques [1, 2, 3] have properties that are favourable in high dimensions.
However, the previously proposed methods provide no guarantee that the
approximated covariance matrix remains non-negative definite.

In this paper we propose a variational formulation aided by the H∞ filter to
avoid the complexity of storage and computation with large dense state error
covariance matrices. Variational methods, which optimize an approximating
distribution to a Bayesian posterior, have become critical in large scale
machine learning because they can handle both high dimensional models and
large amounts of data [12, 4]. Variational methods often approximate the
posterior distribution with a multivariate normal distribution with diagonal
covariance, making it tractable to store and manipulate in high dimensions
[10]. We explore combining the robustness of H∞ and the tractability of
variational inference to provide an algorithm that is efficient and can handle
uncertainty due to approximations.

The outline is as follows: First in Section II we introduce the model along
with traditional KF, then in Sections III and IV we discuss the VI and H∞
filter respectively. In Section V we derive our proposed augmented filter.
Next, in the Section VI we demonstrate the effectiveness of our technique and
finally in Section VII, we provide a discussion of the limitations and
possible future directions.

##  II The Finite-Dimensional Linear-Gaussian Model and Kalman Filter

Many large scale systems, where the problem entails state or parameter
estimation, can be modeled or approximated by a linear system with additive
Gaussian noise. We enumerate observations and update by the index $t$, so that
$\bm{x}_{t}\in\mathbb{R}^{n}$ represents input variables that predicts outputs
$\bm{y}_{t}\in\mathbb{R}^{m}$. We are particularly interested in the case
where the dimension of inputs is much larger than dimension of outputs, $m\ll
n$. We are going to focus on the simplest case where $m=1$, so our measurement
and model update equations becomes,

$\displaystyle\bm{y}_{t}$ $\displaystyle=\bm{x}_{t}^{T}\vTh_{t}+\bm{\eta}_{t}$
(1) $\displaystyle\vTh_{t+1}$ $\displaystyle=\bm{A}_{t}\vTh_{t}+\bm{w}_{t}$
(2)

where $\vTh_{t}\in\mathbb{R}^{n}$ is the model parameter and/or states, that
is to be estimated. The variables $\bm{w}_{t}$ and $\bm{\eta}_{t}$ are assumed
to be i.i.d. zero mean Gaussian noise with covariance matrix $\bm{Q}$ and
$\bm{R}$ respectively. The matrix $\bm{A}_{t}$ is identity and $\bm{w}_{t}=0$
for all $i$ when we perform parameter estimation. We assume a prior on
$\vTh\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\vTh_{0}})$.

Under the model (II), the optimal update scheme based on Bayesian inference is
to use the Kalman Filter. Given a current estimate
$\vTh_{t-1}\sim\mathcal{N}\left(\vMu_{t-1},\bm{P}_{t-1}\right)$, the KF
defines how to update our estimate given $\bm{A}_{t}$, $\bm{x}_{t}$, and
$\bm{y}_{i}$ to arrive at
$\vTh_{i}\sim\mathcal{N}\left(\vMuKF_{t},\bm{P}^{\text{KF}}_{t}\right)$. If we
have a set of efficient covariance matrix approximations, e.g. diagonal
matrices $\bm{\mathcal{P}}$, then the approximate KF problem relies on finding
a $\bm{P}_{t}\in\bm{\mathcal{P}}$ that captures information about
$\bm{P}^{\text{KF}}_{t}$ without introducing bias (difference between expected
estimate and the true value) in the filter. A motivating approach is to choose
$\bm{P}_{i}$ such that
$\bm{P}^{\text{KF}}_{t}\preccurlyeq\bm{P}_{t}\preccurlyeq\bm{P}_{t-1}$.
$\bm{P}_{t}\preccurlyeq\bm{P}_{t-1}$ implies information has been learned
while $\bm{P}^{\text{KF}}_{t}\preccurlyeq\bm{P}_{t}$ implies that the
inference is conservative. However, this approach is too conservative to
enable efficient learning in the case of diagonal matrices. Finding a
satisfactory $\bm{P}_{t}$ is near impossible because the principle axes of
$\bm{P}^{\text{KF}}_{t}$ almost never align with the standard basis so
$\bm{P}_{t-1}=\bm{P}_{t}$ (Theorem 1). This motivates the need for a less
conservative yet robust approach.

###### Theorem 1

For parameter estimation (e.g. $A_{t}=\mathbb{I}$ and $w_{t}=0$) with a random
input vector $\bm{x}_{t}\sim\mathcal{N}\left(\bar{x},\Sigma_{x}\right)$, the
constraint
$\bm{P}^{\text{KF}}_{t}\preccurlyeq\bm{P}_{t}\preccurlyeq\bm{P}_{t-1}\xRightarrow{a.s.}\bm{P}_{t}=\bm{P}_{t-1}$,
when $\bm{P}_{t}$, $\bm{P}_{t-1}$ are diagonal.

###### Proof:

Let $\bm{d}^{t-1}$, $\bm{d}^{t}$ be the diagonals of $\bm{P}_{t-1}$ and
$\bm{P}_{t}$. Then $\bm{P}_{t}\preccurlyeq\bm{P}_{t-1}\implies d^{t}_{j}\leq
d^{t-1}_{j}\forall j$. For the KF,
$\bm{P}^{\text{KF}}_{t}=\bm{P}_{t-1}-\frac{\bm{P}_{t-1}\bm{x}^{T}_{t}\bm{x}_{t}\bm{P}_{t-1}}{\bm{x}_{t}\bm{P}_{t-1}\bm{x}^{T}_{t}+R_{t}}$.
$\bm{P}^{\text{KF}}_{t}\preccurlyeq\bm{P}_{t}\Leftrightarrow\bm{u}\left(\bm{P}_{t-1}-\bm{P}_{t}\right)\bm{u}^{T}\leq\bm{u}\frac{\bm{P}_{t-1}x^{T}_{t}x_{t}\bm{P}_{t-1}}{x_{t}\bm{P}_{t-1}x^{T}_{t}+R_{t}}\bm{u}^{T}$
$\forall\bm{u}\in\mathbb{R}^{n}$. Since $\bm{P}_{t-1}$ and $\bm{P}_{t}$ are
diagonal, the LHS is $\sum_{j}u_{j}^{2}(d^{t-1}_{j}-d^{t}_{j})$. Because
$\bm{P}_{t-1}\bm{x}^{T}_{t}\bm{x}_{t}\bm{P}_{t-1}$ is a rank-1 matrix,
$\exists\bm{u}_{\perp}\in
Null\left(\bm{P}_{t-1}\bm{x}^{T}_{t}\bm{x}_{t}\bm{P}_{t-1}\right)$ s.t. the
RHS is 0. When this holds then,
$\sum_{j}u_{j}^{2}({d}^{t-1}_{j}-{d}^{t}_{j})\leq 0$. Since both $u_{j}^{2}$
and $({d}^{t-1}_{j}-{d}^{t}_{j})$ are always non-negative, then either
$u_{j}^{2}$ or $({d}^{t-1}_{j}-{d}^{t}_{j})$ must be $0$ $\forall j$. Since
$\bm{x}_{t}$ are randomly distributed then the resulting $u_{\perp}$ will
almost surely be all non-zero. Therefore
a.s.,$({d}^{t-1}_{j}={d}^{t}_{j})\implies\bm{P}_{t-1}=\bm{P}_{t}$. ∎

##  III Discrete-time VI Filter

Variational Inference(VI) is an approximate inference method in Bayesian
statistics. Given a model we infer its posterior density by updating our prior
belief about the model with observations. Solving for the exact posterior
$p(\vTh\mid\bm{y})$ is often impossible, necessitating the use of different
types of approximations such as approximating it with samples or a more
tractable distribution. In VI this approximation is done with a distribution
$q(\vTh\mid\bm{\phi})$ which comes from a set of possible distributions
parameterized by $\bm{\phi}$.

For linear Gaussian models, the analytic posterior maybe constructed; however,
the cost of storing $\bm{P}_{t}$ and potentially computing its inverse
$\bm{P}_{t}^{-1}$ for subsequent analysis has motivated us to simplify the
structure of the Gaussian distribution. We assume a variational structure such
that the individual $\vTh_{t}$ elements are independent. This method is called
the mean-field approximation [4], which imposes a diagonal structure on
$\bm{P}_{t}=\text{diag}(\bm{d}_{t})$. This is much more tractable in high
dimensional problems. In this case the variational distribution is
parameterized by the mean and diagonal vector $\bm{d}_{t}$.

Traditional VI seeks to minimize the Kullback–Leibler (KL) divergence between
the approximating variational distribution and the target posterior. For two
distribution $p\left(\vTh\right)$ and $q\left(\vTh\right)$ the KL divergence
is given by,

$\text{D}_{\text{KL}}(p(\vTh)||q(\vTh))=\int
p(\vTh)\log\frac{p\left(\vTh\right)}{q\left(\vTh\right)}d\vTh.$ (3)

The KL divergence stems from information theory and is similar to a distance
between the distributions; it is nonnegative and only zero if the
distributions are the same almost everywhere. Unlike a distance metric,
however, it is not symmetric nor does it obey the triangle inequality. If both
$p$ and $q$ are n-dimensional multivariate Gaussian distributions, with means
and covariances $(\vMuP,\bm{\Sigma}_{p})$ and $(\vMuQ,\bm{\Sigma}_{q})$
respectively, then the KL-divergence becomes

$\displaystyle\text{D}_{\text{KL}}(p(\vTh)||q(\vTh))$ $\displaystyle=$
$\displaystyle\frac{1}{2}\biggl{(}\mathbf{tr}\left(\bm{\Sigma}_{q}^{-1}\bm{\Sigma}_{p}\right)$
$\displaystyle-n+||\vMuP-\vMuQ||^{2}_{\bm{\Sigma}_{q}^{-1}}+\log\frac{|\bm{\Sigma}_{q}|}{|\bm{\Sigma}_{p}|}\biggr{)}.$
(4)

When using the KL divergence for VI there is a choice to minimize either
direction of the KL divergence e.g.
$\text{D}_{\text{KL}}(p(\vTh\mid\bm{y})||q(\vTh\mid\bm{\phi}))$ or
$\text{D}_{\text{KL}}(q(\vTh\mid\bm{\phi})||p(\vTh\mid\bm{y}))$. Choosing
$\text{D}_{\text{KL}}(q(\vTh\mid\bm{\phi})||p(\vTh\mid\bm{y}))$ is common
because it does not require sampling or otherwise manipulating the posterior
$p(\vTh\mid\bm{y})$. This objective is equivalent to maximizing the evidence
lower bound (ELBO) so is often referred to as ELBO optimization. For Gaussian
VI with mean-field approximation we optimize this KL divergence over the
positive space of $\bm{d}_{t}$.

$\displaystyle\mathcal{L}(\bm{d}_{t},\bm{y}_{t},\bm{x}_{t})=\text{D}_{\text{KL}}(q(\vTh\mid\bm{d}_{t})||p(\vTh_{t}\mid\bm{y}_{t}))$
(5)

where $q(\vTh\mid\bm{d}_{t})$ is the Gaussian VI distribution parameterized by
$\bm{d}_{t}$. The $p(\vTh_{t}\mid\bm{y}_{t})$ is the posterior distribution of
the state given the observation. We assume that the $p,q$ have the same mean
so it is not a parameter of the optimization because it can be easily shown
that the optimal mean is the posterior mean. Therefore, we would use the
standard Kalman update to update the mean as data is assimilated.

Similarly, we could also choose the other direction of the KL divergence,
$\text{D}_{\text{KL}}(p(\vTh\mid\bm{y})||q(\vTh\mid\bm{\phi}))$. Within VI
this is know as Expectation Propagation (EP). While EP is more challenging, it
may be preferred when tractable because it is well known that the ELBO
formulation of VI generally underestimates uncertainty [4] compared to EP.
Therefore sequential ELBO-VI is more likely to lead to bias compared to EP-VI.
Fortunately, when the dimension of the streaming data is small, EP can be
computed analytically because the posterior is just a low-rank update to the
mean and covariance. Therefore, for the remainder of this work we will focus
on the EP formulation.

######  Theorem 2 (Optimal $P_{\text{VI-EP}}$)

Suppose the target posterior
$p\left(\vTh\mid\bm{y}\right)\sim\mathcal{N}\left(\vMuKF,\bm{P}_{\text{KF}}\right)$
and the variational distribution is
$q(\vTh\mid\bm{\phi})\sim\mathcal{N}\left(\vMu,\bm{P}(\bm{d})\right)$, where
$\bm{P}(\bm{d})$ is constrained to be a diagonal matrix with diagonal elements
$d_{i}$. Then the EP optimal $\bm{\mu}_{\text{VI-EP}}$ and $\bm{d}_{\text{VI-
EP}}$ are given by

$\displaystyle\left(\bm{\mu}_{\text{VI-EP}},\bm{d}_{\text{VI-EP}}\right)$
$\displaystyle=\operatorname*{\arg\\!\min}_{(\bm{\mu},\bm{d})}\text{D}_{\text{KL}}(p(\vTh\mid\bm{y})||q(\vTh\mid\bm{\mu},\bm{d}))$
$\displaystyle=\left(\vMuKF,\text{diag}(\bm{P}_{\text{KF}})\right)$ (6)

###### Proof:

The optimal $\bm{\mu}$ is found by setting the $\bm{\mu}$’s to be equal and
eliminating the norm term in (4). The optimal $\bm{d}$ is found by setting the
derivative of $\eqref{eq:guass_kl}$ wrt $d_{i}$ to zero and solving for
$d_{i}$. ∎

Figure 1: Plot showing the scaled estimation error response of the VI-EP
filter with dimension of $\vTh$, $n=50$ and dimension of $\bm{y}$, $m=1$. The
scaling is wrt to the standard deviation from the covariance matrix. We see
some state are not converging because their error relative to the standard
deviation is growing and approaches 60$\sigma$ by the end of the simulation.
This indicates that not only are they being poorly estimated but that their
uncertainty is very poorly quantified.

Even EP-VI filter suffers from convergence issue (having bias) when the $\vTh$
dimension is considerably higher than the observation dimension (Fig. 1),
unlike in [10] where both dimensions are comparable. This is because learning
sequentially from a limited number of observations requires capturing
parameter correlation which is cannot be done with a diagonal matrix.
Therefore, the covariance update must be more conservative than the EP update.

With the same diagonal structure of $\bm{P}_{t}$ we considered the general
family of f-divergence to measure the difference between $p$ and $q$. In [7],
the notions of information density $\mathbb{D}(q|p)$ and $\mathbb{L}^{r}$
information psuedometrics,

$\displaystyle\mathbb{L}^{r}=\Big{(}\int
p(\vTh)\Big{|}\log\frac{p\left(\vTh\right)}{q\left(\vTh\right)}\Big{|}^{r}d\vTh\Big{)}^{1/r};\quad
r\geq 1$ (7)

enables us to investigate the performance of this information theory based
f-divergence measure to identify the VI distribution. In fig. 2 we see that
optimizing the information pseudometric for $r=2$, results in an approximation
that captures more of the high posterior probability domain, which dominates
the objective. As a consequence, it also includes regions with lower posterior
probability than traditional VI EP or the ELBO optimization.

Figure 2: Comparing 90% confidence level of true posterior (in black),
traditional VI update based on KL divergence with EP (in red) and with the
ELBO (in blue). VI with the information pseudometric update using $r=2$ (in
green). The pseudometric optimum avoids suppressing probability in domains
that have significant posterior density.

Motivated by this observation, we propose a filter, whose covariance update
step is given by the following,

$\mathbb{L}^{r}$ Filter Update Step: The diagonal posterior covariance
$\bm{P}_{t}(\bm{d})$, is the outcome of the following information pseudometric
optimization,

$\displaystyle\operatorname*{\arg\\!\min}_{\bm{d}_{i}}\int
p(\vTh)\Big{|}\log\frac{p\left(\vTh\right)}{q\left(\vTh\mid\bm{d}\right)}\Big{|}^{r}d\vTh=\bm{d}^{*}\rightarrow\bm{P}_{t}(\bm{d}^{*})$
(8)

For high-dimensional models, the diagonal structure of the $\bm{P}_{t}$ can
severely constrain the movement of estimates towards the true values. Due to
the lack of cross-correlation, information gain only occurs along standard
basis dimensions, causing some estimates to become much more confident than
others. To mitigate this issue, and enable each of the elements in the
estimate of $\vTh_{t}$ to be confident, we draw motivation from the H∞ filter,
where we increase robustness is exchanged for slower information gain. For the
remainder of this work we will focus on the $\mathbb{L}^{2}$ formulation. For
if $p,q$ are both n-dimensional multivariate Gaussian distributions with
common mean $\vMu$ and but different covariances $\bm{\Sigma}_{p}$ and
$\bm{\Sigma}_{q}=\text{diag}(\bm{d})$ respectively, then equation. (8) for
$r=2$ simplifies to,

$\displaystyle\int p(\vTh)\Big{|}\log$
$\displaystyle\frac{p\left(\vTh\right)}{q\left(\vTh\mid\bm{d}\right)}\Big{|}^{2}d\vTh=\frac{n}{2}-\mathbf{tr}\left(\bm{M}\right)+\frac{\mathbf{tr}\left(\bm{M}^{2}\right)}{2}$
$\displaystyle+\frac{1}{4}\left(\mathbf{tr}\left(\bm{M}\right)-\log|\bm{M}|-n\right)^{2}$
(9a)

where $\bm{M}=\bm{\Sigma}_{q}^{-1}(\bm{d})\bm{\Sigma}_{p}$. With the $\bm{d}$
exposed in equations. (9), we optimize the $\mathbb{L}^{2}$ psuedometric with
respect to $\bm{d}$ to calculate $\bm{d}^{*}$. This objective and its
gradients can be efficiently computed when $\bm{\Sigma}_{p}$ is diagonal plus
low rank.

##  IV Discrete-time $H_{\infty}$ filter

We consider the same set of equations (2) and (1), and add that $\vTh_{t}$
itself is the signal of interest rather than its linear projection. The goal
of the H∞ filter[9] is to estimate $\vTh_{t}$, such that the cost function,

$\displaystyle
J=\frac{\sum\limits_{t=0}^{N-1}||\vTh_{t}-\hat{\vTh}_{t}||^{2}_{\bm{S}_{t}}}{||\vTh_{0}-\hat{\vTh}_{0}||^{2}_{\bm{P}_{0}^{-1}}+\sum\limits_{t=0}^{N-1}\Big{(}||\bm{w}_{t}||^{2}_{\bm{Q}^{-1}}+||\bm{\eta}_{t}||^{2}_{\bm{R}^{-1}}\Big{)}}$
(10)

can be made less than $1/\gamma$, where $\bm{S}_{t}$ is a symmetric, positive
definite matrix, chosen to be identity. The posterior covariance matrix of the
state/parameter estimate is $\bm{P}_{t}$. The estimator that achieves this is,

$\displaystyle\bm{P}_{t}$
$\displaystyle=\tilde{\bm{P}}_{t}[\bm{I}-\gamma\bm{S}_{t}\tilde{\bm{P}}_{t}+\bm{x}_{t}\bm{R}^{-1}\bm{x}_{t}^{T}\tilde{\bm{P}}_{t}]^{-1}$
(11a) $\displaystyle\bm{K}_{t}$
$\displaystyle=\bm{P}_{t}\bm{x}_{t}\bm{R}^{-1}$ (11b)
$\displaystyle\hat{\vTh}_{t}$
$\displaystyle=\tilde{\vTh}_{t}+\bm{K}_{t}(\bm{y}_{t}-\bm{x}_{t}^{T}\tilde{\vTh}_{t})$
(11c)

where $\tilde{\bm{P}}_{t}$ is the prior covariance and $\hat{\vTh}_{t}$ is the
posterior estimate from the prior $\tilde{\vTh}_{t}$. The condition,
$\tilde{\bm{P}}_{t}^{-1}-\gamma\bm{I}+\bm{x}_{t}\bm{R}^{-1}\bm{x}_{t}^{T}\succ
0$ must hold for the estimator to provide solution to the $H_{\infty}$
problem.

##  V Augmented H∞ Filter Update Step

The $\gamma$ parameter in the H∞ Filter determines the upper bound on the H∞
cost function in (10). We seek to calculate $\gamma$ to minimize the
discrepancy between the Kalman gains resulting from the regular Kalman update
and the H∞ Kalman gain calculated after the VI or the $\mathbb{L}^{2}$ stage.
Doing so we are ensuring the diagonal structure over the $\bm{P}_{t}$ while
improving the robustness of the estimates. The proposed $\mathbb{L}^{2}$-VI H∞
filter algorithm is presented in algorithm. 1. The algorithm drastically
simplifies for static parameter estimation for which $\bm{A}_{t}$ is identity
and $\bm{w}_{t}=0$.

Figure 3: Illustration of a robust Variational Filter. The Filter is composed
of three primary steps. First the Kalman update assimilates data using
Bayesian Inference to compute the true update
$(\mu_{t+1},\bm{P}^{\text{KF}}_{t+1})$. Then Variational Inference is used to
find the best low memory approximation $\bm{P}^{VI}_{t+1}$ (using either EP or
$\mathbb{L}^{r}$). Finally, comparing the true update and variational
approximation a final update $\bm{P}^{H_{\infty}}_{t+1}$is optimized to
robustly assimilate the next observation at $\bm{x}_{t+1}$

$\bm{P}^{H_{\infty}}_{t-1}$: Diagonal posterior matrix at step (t-1)
$\bm{P}^{H_{\infty}}_{t}$: Diagonal posterior matrix at step t Calculate:

$\displaystyle\bm{K}_{t}$
$\displaystyle=(\bm{A}_{t-1}\bm{P}^{H_{\infty}}_{t-1}\bm{A}_{t-1}^{T}+\bm{Q})\bm{x}_{t}$
$\displaystyle\times\Big{[}\bm{x}_{t}^{T}(\bm{A}_{t-1}\bm{P}^{H_{\infty}}_{t-1}\bm{A}_{t-1}^{T}+\bm{Q})\bm{x}_{t}+\bm{R}\Big{]}^{-1}$
$\displaystyle\bm{P}_{t}^{\text{KF}}$
$\displaystyle=(\bm{I}-\bm{K}_{t}\bm{x}_{t}^{T})(\bm{A}_{t-1}\bm{P}^{H_{\infty}}_{t-1}\bm{A}_{t-1}^{T}+\bm{Q})$
(12)

Optimization 1:

$\displaystyle\bm{d}^{\mathbb{L}^{r}}=\operatorname*{\arg\\!\min}_{\bm{d}_{t}}\int
p(\vTh\mid\bm{y}_{t})\Big{|}\text{log}\frac{p(\vTh\mid\bm{y})}{q(\theta\mid\bm{d}_{t})}\Big{|}^{r}d\vTh\rightarrow\bm{P}^{\mathbb{L}^{r}}_{t}(\bm{d})$
(13)

Propagate $\bm{P}^{\mathbb{L}^{r}}_{t}(\bm{d})$:
$\tilde{\bm{P}}^{\mathbb{L}^{r}}_{t+1}=(\bm{A}_{t+1}\bm{P}^{\mathbb{L}^{r}}_{t}\bm{A}_{t+1}^{T}+\bm{Q})$

Compute the Kalman update for the next step using $\bm{P}_{t}^{\text{KF}}$:

$\displaystyle\bm{K}^{\text{KF}}_{t+1}$
$\displaystyle=(\bm{A}_{t}\bm{P}^{\text{KF}}_{t}\bm{A}_{t}^{T}+\bm{Q})\bm{x}_{t+1}$
$\displaystyle\times\Big{[}\bm{x}_{t+1}^{T}(\bm{A}_{t}\bm{P}^{\text{KF}}_{t}\bm{A}_{t}^{T}+\bm{Q})\bm{x}_{t+1}+\bm{R}\Big{]}^{-1}$
(14)

Compute the expression for $\bm{K}^{\text{H}_{\infty}}_{t+1}$:

$\displaystyle\bm{K}^{\text{H}_{\infty}}_{t+1}=$
$\displaystyle\tilde{\bm{P}}^{\mathbb{L}^{r}}_{t+1}[\bm{I}-\gamma\tilde{\bm{P}}^{\mathbb{L}^{r}}_{t+1}+\bm{x}_{t}\bm{R}^{-1}\bm{x}_{t}^{T}\tilde{\bm{P}}^{\mathbb{L}^{r}}_{t+1}]^{-1}$
$\displaystyle\times\bm{x}_{t}\bm{R}^{-1}$ (15)

Optimization 2:

$\displaystyle\gamma^{*}_{t}$
$\displaystyle=\operatorname*{\arg\\!\min}_{\gamma}||\bm{K}^{\text{H}_{\infty}}_{t+1}-\bm{K}^{\text{KF}}_{t+1}||_{2}$
$\displaystyle\bm{P}^{H_{\infty}}_{t}$
$\displaystyle=\bm{P}^{\mathbb{L}^{r}}_{t}(\bm{d})[\bm{I}-\gamma_{t}^{*}\bm{P}^{\mathbb{L}^{r}}_{t}(\bm{d})+\bm{x}_{t}\bm{R}^{-1}\bm{x}_{t}^{T}\bm{P}^{\mathbb{L}^{r}}_{t}(\bm{d})]^{-1}$
(16)

Algorithm 1 Algorithm of Augmented H∞ Filter Update

When the Optimization step 1 is replaced by the EP objective, we recover the
VI-H∞ filtering. The fig. 3 illustrates how this augmented filtering is
performed. It is important to note that all these steps can be done in a
computationally efficient manner, assuming $\bm{A}_{t}$ is sparse, without
having to manipulate dense matrices. When we assume that $\bm{A}_{t}$ is the
identity, as is the case for parameter estimation, computing is highly
efficient, even when working with $\bm{P}_{t}^{\text{KF}}$ because it is a
diagonal plus low-rank matrix.

##  VI Results

Figure 4: Mean Square Error for the five different filtering algorithms with
varying problem dimensions after 1000 sequentially observed data points. The
error bars correspond to 93% confidence intervals. Two main observations are
that 1) the $\mathbb{L}^{2}$ formulation of the VI problem outperforms
standard EP VI and 2) the two $H_{\infty}$ filters under-perform, but as we
will see are not overconfident so are more trustworthy. There is a small
perturbation on the different filter’s dimension coordinates to improve
readability.  Figure 5: Worst Case Scaled Error (e.g. the largest absolute
error divided by estimated standard deviation) for the five filters. The error
bars are 93% confidence intervals. We see that 1) the variational filters
without added robustness can be significantly overconfident 2) the two
$H_{\infty}$ filters, particularly the $\mathbb{L}^{2}$ $H_{\infty}$ better
capture the worst case estimate uncertainty. These methods trade robustness
and lower bias for slower convergence. There is a small perturbation on the
different filter’s dimension coordinates to improve readability.

We consider comparing 5 different methods for parameter estimation. They are:
the Kalman filter, the variational inference filter using the EP objective,
the augmented VI-H∞ filter, the $\mathbb{L}^{2}$ filter, and the augmented
$\mathbb{L}^{2}$-H∞ filter. Without loss of generality, we assume no system
dynamics and no process noise, focusing only on the sequential update in the
parameter estimate $\hat{\vTh_{t}}$, conditioned on the data $\bm{y}_{t}$. We
consider the case when dimension of y is set to be 1, emulating our systems of
interest, with much lower observation dimension compared with the parameter
dimension. Our observation vector is
$x_{t}\sim\mathcal{N}\left(\vXBar,0.5\bm{I}\right)$. For each problem a
$\vXBar$ was selected from $\mathcal{N}\left(0,\bm{I}\right)$. The
observations had additive measurement noise with variance $0.1$. The prior on
$\theta_{0}$ was $\mathcal{N}\left(0,\bm{I}\right)$. The filter observes
$1000$ time steps. This choice of problem is more difficult than if we had
included process noise or a larger dimensional observation. Process noise
allows estimation errors to become less important over time while larger
dimensional observation spaces allow for the update to include more
information without having to try to sequentially assimilate that same
information over multiple observation steps as in our case.

In fig. 4, we compare the MSE estimate of the 5 filters with varying dimension
of the problem size from 2 to 512 in powers of 2. The mean and range of MSE
was estimated from 32 randomized problems for each dimension. The figure also
shows the uncertainty in those estimates. The Kalman filter response with no
imposed structure on the $\bm{P}_{t}$ matrix, performs the best of all, as
expected since it is optimal. The availability of cross-correlation terms
allows for the lowest MSE estimates. It is interesting to note that the
information pseudo-metric based filter outperforms the traditional VI filter,
with both having diagonal $\bm{P}_{t}$. The other two proposed augmented H∞
filters both under-perform with respect to the MSE estimates. However, it
should be pointed out that the gap between the MSEs of the VI and $H_{\infty}$
filters seems to decrease respect to the increase in dimension of the problem.
We further see that the MSE estimates fluctuates the least among them too. We
assert that this will lead to better convergence characteristics of the
posterior with increase in the problem size at the cost of time of
convergence, which is suitable for many parameter estimation problems. The
diagonal approximation of the $\bm{P}_{t}$ certainly introduces an error with
respect the true posterior. In fig. 5 we show the worst case scaled error
performance of the 5 filters. The scaling with respect to the estimated
standard deviation for each parameter and the worst performance was selected.
When this metric is larger than the Kalman filter, then that filtering
algorithm is unjustifiably confident in its estimate. Conversely when it is
smaller then it is under confident. Among the filters with $\bm{P}_{t}$
diagonal, we see that $\mathbb{L}^{2}$-H∞ outperforms the rest. Without the H∞
correction the VI filter and $\mathbb{L}^{2}$ are expected to performing worse
than the ones with H∞ correction, which is shown in the figure. This
highlights that the $\mathbb{L}^{2}$-H∞ is very robust meaning that it’s
predictions are trustworthy.

##  VII Conclusion

The standard KF becomes exceedingly memory intensive as the dimension of the
underlying state space increases. Variants of the KF have been proposed to
reduce the dimension of the system, thus making implementation in high
dimensions possible. The reduced rank KF project the state/parameter vector of
the model onto a lower dimensional subspace. The success of this approach
depend on meticulously choosing the reduction operator which may not be
possible when the effective problem dimension is large.

In this paper, we propose a reduced memory filter based upon variational
inference with a information psuedometric and H∞ filter in order to resolve
the storage and computational issue related to the error covariance matrix
while retaining robustness. With VI, we approximate the full rank $n\times n$
covarianace matrix, with a matrix characterized by $n$ elements along its
diagonal. The storage and computational complexity of the operations is now
$O\left(n\right)$. In order to test our proposed filters we consider test
cases with increasing problem size to analyze its impact of our method for
large dimensional problems. Our filters, specifically the $\mathbb{L}^{2}$-H∞
exhibits slower learning, but enables information update along all the
directions of the state space, keeping the worst case performance better than
other filters. In our future work we would like to analyze properties of our
proposed filter and experiment with real-world high dimensional problems,
particularly those emerging in Machine Learning.

## References

  * [1] H. Auvinen, J.. Bardsley, H. Haario and T. Kauranne  “The variational Kalman filter and an efficient implementation using limited memory BFGS”  In _International Journal for Numerical Methods in Fluids_ 64.3 Wiley, 2009, pp. 314–335  DOI: 10.1002/fld.2153
  * [2] Harri Auvinen, Johnathan M Bardsley, Heikki Haario and Tuomo Kauranne  “Large-scale Kalman filtering using the limited memory BFGS method”  In _Electronic Transactions on Numerical Analysis_ 35 Institute of Computational Mathematics, 2009, pp. 217–233
  * [3] Johnathan M. Bardsley, Albert Parker, Antti Solonen and Marylesa Howard  “Krylov space approximate Kalman filtering”  In _Numerical Linear Algebra with Applications_ 20.2 Wiley, 2011, pp. 171–184  DOI: 10.1002/nla.805
  * [4] David M Blei, Alp Kucukelbir and Jon D McAuliffe  “Variational inference: A review for statisticians”  In _Journal of the American statistical Association_ 112.518 Taylor & Francis, 2017, pp. 859–877
  * [5] Mark A. Cane et al.  “Mapping tropical Pacific sea level: Data assimilation via a reduced state space Kalman filter”  In _Journal of Geophysical Research: Oceans_ 101.C10 American Geophysical Union (AGU), 1996, pp. 22599–22617  DOI: 10.1029/96jc01684
  * [6] Dick P. Dee  “Simplification of the Kalman filter for meteorological data assimilation”  In _Quarterly Journal of the Royal Meteorological Society_ 117.498 Wiley, 1991, pp. 365–384  DOI: 10.1002/qj.49711749806
  * [7] Jed A. Duersch and Thomas A. Catanach  “Generalizing Information to the Evolution of Rational Belief”  In _Entropy_ 22.1, 2020
  * [8] Nikola B Kovachki and Andrew M Stuart  “Ensemble Kalman inversion: a derivative-free technique for machine learning tasks”  In _Inverse Problems_ 35.9 IOP Publishing, 2019, pp. 095005
  * [9] Dan Simon  “Optimal state estimation: Kalman, H$\infty$, and nonlinear approaches. Hoboken”  In _NJ: John Wiley and Sons_ 10, 2006, pp. 0470045345
  * [10] Christo Kurisummoottil Thomas and Dirk Slock  “Gaussian variational Bayes Kalman filtering for dynamic sparse Bayesian learning”  In _IEEE 5th Intl. Conf. on Time Ser. and Forecast.(ITISE)_ , 2018
  * [11] Arto Voutilainen et al.  “A filtering approach for estimating lake water quality from remote sensing data”  In _International Journal of Applied Earth Observation and Geoinformation_ 9.1 Elsevier BV, 2007, pp. 50–64  DOI: 10.1016/j.jag.2006.07.001
  * [12] Cheng Zhang, Judith Butepage, Hedvig Kjellstrom and Stephan Mandt  “Advances in variational inference”  In _IEEE transactions on pattern analysis and machine intelligence_ 41.8 IEEE, 2018, pp. 2008–2026