# Improved Calibration of Numerical Integration Error in Sigma-Point Filters

Jakub Prüher, Toni Karvonen, Chris J. Oates, Ondřej Straka and Simo Särkkä J.
Prüher and O. Straka are with Department of Cybernetics, University of West
Bohemia, Czech Republic.T. Karvonen and S. Särkkä are with Department of
Electrical Engineering and Automation, Aalto University, Finland.C. J. Oates
is with School of Mathematics, Statistics & Physics, Newcastle University and
the Alan Turing Institute, United Kingdom.

###### Abstract

The sigma-point filters, such as the UKF, which exploit numerical quadrature
to obtain an additional order of accuracy in the moment transformation step,
are popular alternatives to the ubiquitous EKF. The classical quadrature rules
used in the sigma-point filters are motivated via polynomial approximation of
the integrand, however in the applied context these assumptions cannot always
be justified. As a result, quadrature error can introduce bias into estimated
moments, for which there is no compensatory mechanism in the classical sigma-
point filters. This can lead in turn to estimates and predictions that are
poorly calibrated. In this article, we investigate the Bayes–Sard quadrature
method in the context of sigma-point filters, which enables uncertainty due to
quadrature error to be formalised within a probabilistic model. Our first
contribution is to derive the well-known classical quadratures as special
cases of the Bayes–Sard quadrature method. Based on this, a general-purpose
moment transform is developed and utilised in the design of novel sigma-point
filter, which explicitly accounts for the additional uncertainty due to
quadrature error. Numerical experiments on a challenging tracking example with
misspecified initial conditions show that the additional uncertainty
quantification built into our method leads to better-calibrated state
estimates with improved RMSE.

###### Index Terms:

Kalman filters, Bayesian quadrature, quantification of uncertainty, sigma-
points, Gaussian processes.

##  I Introduction

This article is concerned with quantification of uncertainty associated with
sigma-point approximations, which are widely employed in nonlinear local
filtering algorithms, such as the unscented Kalman filter (UKF). The goal of
filtering algorithms is to estimate the state of a dynamical stochastic system
based on all measurements obtained until the present. The applications of
filters are manifold, ranging from global positioning [1], object tracking [2,
3], simultaneous localization and mapping [4] to weather forecasting [5] and
finance [6].

Nonlinear filtering algorithms can be categorised either as local or global.
The global filters, such as particle filters [7] or point-mass filters [8],
keep track of the whole, potentially multi-modal, state posterior, which comes
with increased computational demands. The local filters, on the other hand,
have lower computational load at the cost of introducing more restrictive
assumptions. Instead of keeping track of the whole state posterior, the local
filters only work with the first two statistical moments of the state and the
measurement.

For nonlinear systems and/or measurements, the moments are defined by
intractable integrals that have to be approximated using numerical
quadratures, also known as the sigma-point rules (which is where the filters
get their name). The classical quadrature rules, such as the Gauss–Hermite
rule, are designed with the assumption that the nonlinear integrand is well-
approximated with a polynomial of a given maximal degree. Since these
assumptions are almost never met in practice, there will always be a
quadrature error involved. Standard sigma-point filters do not attempt to
compensate for this source of error, and in practice this can lead to
estimates and predictions that are biased and over-confident [9, 10].

The idea to perform formal probabilistic uncertainty quantification in the
numerical integration context can be traced back to [11], with the name
_Bayesian quadrature_ (BQ) being coined later in [12]. The BQ has in recent
years received much attention in the probabilistic numerics community [13, 14,
15, 16]. The BQ approach posits that the integrand can be modelled by a
stochastic process defined on the domain of integration. This model is
subsequently refined by conditioning on point-wise evaluations of the
integrand which induces a posterior distribution over the value of the
integral. The posterior mean of this distribution is point estimate of the
value of the integral while the posterior variance expresses the integration
error.

Applications of the BQ in nonlinear filtering have appeared previously in [17,
18] with encouraging results. These BQ-based filters do not generally coincide
with any classical sigma-point filter, such as the UKF or Gauss–Hermite Kalman
filter (GHKF), and tend to be rather sensitive to specification of the
stochastic process model for the integrand. It has been shown that classical
sigma-point rules can be cast as _degenerate_ BQ rules [17, 13] (see [19] for
spline methods). This is to say that the variance associated to the integral
vanishes, being thus of no use in modelling integration error.

In this article we utilise the recently proposed Bayes–Sard quadrature [20]
for the design of novel sigma-point filters, which can be viewed as
probabilistic versions of the well-known sigma-point filters. Namely, under
certain conditions, the Bayes–Sard quadrature allows us to recover the
classical sigma-point rules and at the same time endow the sigma-point rule
with non-degenerate probabilistic output. We thus obtain versions of standard
sigma-point filters that are, to some extent, capable of accounting for
numerical integration error in filtering by (in most cases) inflating the
error covariance. In some cases, such covariance inflation is known to improve
stability of nonlinear Kalman filters; see for instance [21, Remark 1], [22,
Section 3.3], and [23, Section V.C].

The rest of the article is structured as follows. In Section II, we formally
outline the nonlinear filtering problem and the nature of sigma-point
approximations. Section III identifies the moment transformation problem as
the central issue in local filtering and describes the structure of sigma-
point moment transforms. The Bayes–Sard quadrature is formalised in Section
IV, which is later used in Section V to design the Bayes–Sard quadrature
moment transform. The numerical experiments are contained in Section VI and
Section VII concludes the article.

##  II Sigma-Point Filtering

This section is devoted to the sigma-point filters, which are a subset of
local filtering algorithms characterised by their reliance on a Gaussian
approximation together with a numerical quadrature method. Let the stochastic
dynamical system and the process by which its state is observed be described
by the state-space model

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

where the function $\bm{\mathbf{f}}:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{x}}$
is the system dynamics,
$\bm{\mathbf{h}}:\mathbb{R}^{d_{x}}\to\mathbb{R}^{d_{z}}$ is the measurement
model, $\bm{\mathbf{x}}_{k}\in\mathbb{R}^{d_{x}}$ is the latent state vector
and $\bm{\mathbf{z}}_{k}\in\mathbb{R}^{d_{z}}$ is the measurement vector. Both
the process noise
$\bm{\mathbf{q}}_{k-1}\sim\mathrm{N}\pqty{\bm{\mathbf{0}},\,\bm{\mathbf{Q}}missing}$
and the measurement noise
$\bm{\mathbf{r}}_{k}\sim\mathrm{N}\pqty{\bm{\mathbf{0}},\,\bm{\mathbf{R}}missing}$
are zero-mean white Gaussian sequences, independent of each other and
independent of the system initial condition
$\bm{\mathbf{x}}_{0}\sim\mathrm{N}\pqty{\bm{\mathbf{m}}^{x}_{0},\,\bm{\mathbf{P}}^{x}_{0}}$.

The Bayesian formulation of the filtering problem can be summarized by the
following two general relations. The state posterior is

$p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{z}}_{1:k}\right.}\propto
p\pqty{\left.\bm{\mathbf{z}}_{k}\,\middle|\,\bm{\mathbf{x}}_{k}\right.}p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{z}}_{1:k-1}\right.},$
(3)

where the likelihood
$p\pqty{\left.\bm{\mathbf{z}}_{k}\,\middle|\,\bm{\mathbf{x}}_{k}\right.}$ is
obtained from the measurement model 2 and
$\bm{\mathbf{z}}_{1:k}\triangleq\Bqty{\bm{\mathbf{z}}_{1},\,\ldots,\,\bm{\mathbf{z}}_{k}}$.
The predictive density is given by the Chapman–Kolmogorov equation

$p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{z}}_{1:k-1}\right.}=\int\limits\\!p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{x}}_{k-1}\right.}p\pqty{\left.\bm{\mathbf{x}}_{k-1}\,\middle|\,\bm{\mathbf{z}}_{1:k-1}\right.}\,\mathrm{d}\bm{\mathbf{x}}_{k-1},$
(4)

where the transition density
$p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{x}}_{k-1}\right.}$ is
obtained from the system dynamics 1.

A vast majority of the well-known local filters, such as EKF, UKF and GHKF,
can be recovered from the Bayesian formulation under a Gaussian approximation
of the joint density of the state and measurement. That is, when the density
$p\pqty{\left.\bm{\mathbf{x}}_{k},\bm{\mathbf{z}}_{k}\,\middle|\,\bm{\mathbf{z}}_{1:k-1}\right.}=p\pqty{\left.\bm{\mathbf{z}}_{k}\,\middle|\,\bm{\mathbf{x}}_{k}\right.}p\pqty{\left.\bm{\mathbf{x}}_{k}\,\middle|\,\bm{\mathbf{z}}_{1:k-1}\right.}$
is approximated by a Gaussian density of the form

$\mathrm{N}\pqty{\left.\bmqty{undef}\,\middle|\,\bmqty{undef},\,\bmqty{undef}\right.},$
(5)

then the mean and covariance of the state posterior 3 have analytical form,
given by 111Note that,
$\bm{\mathbf{m}}^{x}_{k|k}\triangleq\mathbb{E}_{\bm{\mathbf{x}}}\bqty{\bm{\mathbf{x}}_{k}\mid\bm{\mathbf{z}}_{1:k}}$
and
$\bm{\mathbf{P}}^{x}_{k|k}\triangleq\mathbb{E}_{\bm{\mathbf{x}}}\bqty{(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})^{\top}\mid\bm{\mathbf{z}}_{1:k}}$

$\displaystyle\bm{\mathbf{m}}^{x}_{k|k}$
$\displaystyle=\bm{\mathbf{m}}^{x}_{k|k-1}+\bm{\mathbf{G}}_{k}\pqty{\bm{\mathbf{z}}_{k}-\bm{\mathbf{m}}^{z}_{k|k-1}},$
(6) $\displaystyle\bm{\mathbf{P}}^{x}_{k|k}$
$\displaystyle=\bm{\mathbf{P}}^{x}_{k|k-1}-\bm{\mathbf{G}}_{k}\bm{\mathbf{P}}^{z}_{k|k-1}\bm{\mathbf{G}}_{k}^{\top},$
(7)

where
$\bm{\mathbf{G}}_{k}=\bm{\mathbf{P}}^{xz}_{k|k-1}\pqty\big{\bm{\mathbf{P}}^{z}_{k|k-1}}^{-1}$
is the Kalman gain. The predictive moments of the state,
$\bm{\mathbf{m}}^{x}_{k|k-1}$ and $\bm{\mathbf{P}}^{x}_{k|k-1}$, and the
moments of measurements, $\bm{\mathbf{m}}^{z}_{k|k-1}$,
$\bm{\mathbf{P}}^{z}_{k|k-1}$ and $\bm{\mathbf{P}}^{xz}_{k|k-1}$, are defined
as integrals of the form

$\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}\triangleq\int\limits\\!\bm{\mathbf{g}}({\bm{\mathbf{x}}})\mathrm{N}\pqty{\left.{\bm{\mathbf{x}}}\,\middle|\,\bm{\mathbf{m}},\,\bm{\mathbf{P}}\right.}\,\mathrm{d}{\bm{\mathbf{x}}}.$
(8)

Table I shows which quantities have to be substituted for
$\bm{\mathbf{g}}(\bm{\mathbf{x}})$, $\bm{\mathbf{x}}$, $\bm{\mathbf{m}}$ and
$\bm{\mathbf{P}}$ to obtain any of the above predictive moments. Since the
function $\bm{\mathbf{g}}$ being integrated is nonlinear in each case, these
integrals cannot be typically computed analytically and some type of
approximation needs to be employed222When the model functions
$\bm{\mathbf{f}}$ and $\bm{\mathbf{h}}$ are linear, the integrals can be
computed analytically and the eqs. 6 and 7 become identical to the Kalman
filter update.. Each local nonlinear filter is distinguished solely by the
type of integral approximation it uses. For example, the EKF employs the first
order Taylor expansion to linearise $\bm{\mathbf{g}}$ in the vicinity of
$\bm{\mathbf{m}}$, which in turn facilitates analytic tractability of the
moment integrals. On the other hand, the sigma-point filters, such as the UKF
and the GHKF, leverage numerical quadrature for approximation of the integral.
Since quadratures are typically designed to be used with standard Gaussian,
the integrals of the form 8 need to be converted by employing a stochastic
decoupling substitution
${\bm{\mathbf{x}}}^{(n)}=\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}}^{(n)}$,
which leads to an approximation

$\displaystyle\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}$
$\displaystyle\approx\sum_{n=1}^{N}w_{n}\bm{\mathbf{g}}(\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}}^{(n)})$
$\displaystyle=\sum_{n=1}^{N}w_{n}\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}^{(n)}),$
(9)

where ${\bm{\mathbf{\xi}}}^{(n)}$ denotes the $n$-th unit sigma-point,
$w_{n}\in\mathbb{R}$ is the $n$-th weight, $N$ is the total number of sigma-
points, $\bm{\mathbf{L}}$ is a matrix factor such that
$\bm{\mathbf{P}}=\bm{\mathbf{L}}\bm{\mathbf{L}}^{\top}$ and
$\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})\triangleq\bm{\mathbf{g}}(\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}})$.
Note that various quadrature rules are distinguished by the different weights
and sigma-points they prescribe to satisfy various optimality criteria.

Moment | $\>\bm{\mathbf{g}}(\bm{\mathbf{x}})\,$ | $\>\bm{\mathbf{x}}\,$ | $\>\bm{\mathbf{m}}\,$ | $\>\bm{\mathbf{P}}\,$
---|---|---|---|---
$\,\bm{\mathbf{m}}^{x}_{k|k-1}\>$ | $\>\bm{\mathbf{f}}(\bm{\mathbf{x}}_{k-1})\,$ | $\>\bm{\mathbf{x}}_{k-1}\,$ | $\>\bm{\mathbf{m}}^{x}_{k-1|k-1}\,$ | $\>\bm{\mathbf{P}}^{x}_{k-1|k-1}\,$
$\,\bm{\mathbf{P}}^{x}_{k|k-1}\>$ | $\>\Delta\bm{\mathbf{f}}\Delta\bm{\mathbf{f}}^{\top}\,$ | $\>\bm{\mathbf{x}}_{k-1}\,$ | $\>\bm{\mathbf{m}}^{x}_{k-1|k-1}\,$ | $\>\bm{\mathbf{P}}^{x}_{k-1|k-1}\,$
$\,\bm{\mathbf{m}}^{z}_{k|k-1}\>$ | $\>\bm{\mathbf{h}}(\bm{\mathbf{x}}_{k})\,$ | $\>\bm{\mathbf{x}}_{k}\,$ | $\>\bm{\mathbf{m}}^{x}_{k|k-1}\,$ | $\>\bm{\mathbf{P}}^{x}_{k|k-1}\,$
$\,\bm{\mathbf{P}}^{z}_{k|k-1}\>$ | $\>\Delta\bm{\mathbf{h}}\Delta\bm{\mathbf{h}}^{\top}\,$ | $\>\bm{\mathbf{x}}_{k}\,$ | $\>\bm{\mathbf{m}}^{x}_{k|k-1}\,$ | $\>\bm{\mathbf{P}}^{x}_{k|k-1}\,$
$\,\bm{\mathbf{P}}^{xz}_{k|k-1}\>$ | $\>\Delta{\bm{\mathbf{x}}}\Delta\bm{\mathbf{h}}^{\top}\,$ | $\>\bm{\mathbf{x}}_{k}\,$ | $\>\bm{\mathbf{m}}^{x}_{k|k-1}\,$ | $\>\bm{\mathbf{P}}^{x}_{k|k-1}\,$
Table I:  Quantities that need to be substituted into the Gaussian integral
(8) in order to obtain every predictive moment necessary to compute the
moments of the state posterior. The following shorthand notation is used:
$\Delta\bm{\mathbf{f}}=\bm{\mathbf{f}}(\bm{\mathbf{x}}_{k-1})-\bm{\mathbf{m}}^{x}_{k|k-1}$,
$\Delta\bm{\mathbf{h}}=\bm{\mathbf{h}}(\bm{\mathbf{x}}_{k})-\bm{\mathbf{m}}^{z}_{k|k-1}$,
$\Delta{\bm{\mathbf{x}}}=\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k-1}$.

##  III Sigma-Point Moment Transforms

From the above exposition, it is apparent that the central issue in local
filtering is the design of the so-called moment transformations, which
generate approximations to the moments of a random variable under a possibly
nonlinear transformation.

Let ${\bm{\mathbf{x}}}\in\mathbb{R}^{D}$ be an input Gaussian random variable
and $\bm{\mathbf{y}}\in\mathbb{R}^{E}$ an output random variable defined by

$\bm{\mathbf{y}}=\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\quad{\bm{\mathbf{x}}}\sim\mathrm{N}\pqty{\bm{\mathbf{m}},\,\bm{\mathbf{P}}}.$
(10)

If the transformation $\bm{\mathbf{g}}$ is nonlinear, the joint density
$p({\bm{\mathbf{x}}},\,\bm{\mathbf{y}})$ will be non-Gaussian in general.
However, there are many applied situations where $\bm{\mathbf{g}}$ is
approximately linear in the region where probability mass is concentrated. In
such situations the principal error term in the moment transform is numerical
quadrature error. This error is the focus of our present work and, therefore,
in what follows we proceed under the assumption that the Gaussian
approximation

$\mathrm{N}\pqty{\left.\bmqty{undef}\,\middle|\,\bmqty{undef},\,\bmqty{undef}\right.}$
(11)

of $p({\bm{\mathbf{x}}},\,\bm{\mathbf{y}})$ can be justified. In this setting
the moment transformation then reduces to computing the output mean
$\bm{\mathbf{\mu}}$, covariance $\bm{\mathbf{\Pi}}$ and cross-covariance
$\bm{\mathbf{C}}$ as accurately as possible, when supplied with the input
moments, $\bm{\mathbf{m}}$ and $\bm{\mathbf{P}}$. This is a specific instance
of _uncertainty propagation_ ; see for example [24].

In this article we focus on the sigma-point approximations, exemplified by eq.
9, to the moment integrals in Table I. The well-known classical
approximations, such as the Gauss–Hermite, the spherical-radial and the
unscented transform, are conventionally written in the form

$\displaystyle\bm{\mathbf{\mu}}\approx\hat{\bm{\mathbf{\mu}}}$
$\displaystyle=\sum_{n=1}^{N}w_{n}\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}^{(n)}),$
(12) $\displaystyle\bm{\mathbf{\Pi}}\approx\hat{\bm{\mathbf{\Pi}}}$
$\displaystyle=\sum_{n=1}^{N}w_{n}\pqty\big{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}^{(n)})-\hat{\bm{\mathbf{\mu}}}}\pqty\big{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}^{(n)})-\hat{\bm{\mathbf{\mu}}}}^{\top},$
(13) $\displaystyle\bm{\mathbf{C}}\approx\hat{\bm{\mathbf{C}}}$
$\displaystyle=\bm{\mathbf{L}}\sum_{n=1}^{N}w_{n}{\bm{\mathbf{\xi}}}^{(n)}\pqty\big{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}^{(n)})-\hat{\bm{\mathbf{\mu}}}}^{\top},$
(14)

which, under the assumption that $\sum w_{n}=1$ and $\sum
w_{n}{\bm{\mathbf{\xi}}}^{(n)}=0$, we will prefer to write using the matrix
notation as

$\displaystyle\hat{\bm{\mathbf{\mu}}}$
$\displaystyle=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{w}},$ (15a)
$\displaystyle\hat{\bm{\mathbf{\Pi}}}$
$\displaystyle=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{W}}\bm{\mathbf{Y}}-\hat{\bm{\mathbf{\mu}}}\hat{\bm{\mathbf{\mu}}}^{\top},$
(15b) $\displaystyle\hat{\bm{\mathbf{C}}}$
$\displaystyle=\bm{\mathbf{L}}\bm{\mathbf{\Xi}}\bm{\mathbf{W}}_{c}\bm{\mathbf{Y}},$
(15c)

where
$\bm{\mathbf{\Xi}}=\bqty\big{\matrixquantity{{\bm{\mathbf{\xi}}}^{(1)}&\ldots&{\bm{\mathbf{\xi}}}^{(N)}}}$
and the matrix of integrand evaluations is given by
$\bqty{\bm{\mathbf{Y}}}_{ne}\triangleq\tilde{g}_{e}({\bm{\mathbf{\xi}}}^{(n)})$,
where $e$ indexes outputs of $\tilde{\bm{\mathbf{g}}}$ and
$\bqty{\,\cdot\,}_{ne}$ denotes the matrix element at position $(n,e)$. The
vector $\bm{\mathbf{w}}$ contains the weights and
$\bm{\mathbf{W}}=\bm{\mathbf{W}}_{c}=\mathrm{diag}\pqty{\,\bm{\mathbf{w}}\,}$
for any classical sigma-point moment transform. Each moment transform uses a
different set of sigma-points and weights.

We next discuss the unscented transform and the Gauss–Hermite quadrature in
detail.

###  III-A Unscented Transform

The unscented transform (UT) of $D$-dimensional input uses $N=2D+1$ sigma-
points, which exploit symmetry of the Gaussian distribution, given, for
$d=1,\ldots,D$, by

${\bm{\mathbf{\xi}}}^{(0)}=\bm{\mathbf{0}},\quad{\bm{\mathbf{\xi}}}^{(d)}=\sqrt{c}\,\bm{\mathbf{e}}_{d},\quad{\bm{\mathbf{\xi}}}^{(D+d)}=-\sqrt{c}\,\bm{\mathbf{e}}_{D+d},$
(16)

where $\bm{\mathbf{e}}_{d}$ is the standard unit vector and $c=D+\kappa$ for a
scaling parameter $\kappa$. The weights are defined as

$w_{0}=\frac{\kappa}{c},\quad w_{d}=\frac{1}{2c},\quad w_{D+d}=\frac{1}{2c}.$
(17)

This selection of sigma-points and weights yields a quadrature rule that
integrates exactly all polynomials of (total) degree at most three; the
derivation is essentially contained in the proof of Theorem 1. The _spherical-
radial rule_ , which is used in the cubature Kalman filter (CKF) [25], is
equivalent to the UT with $\kappa=0$; it therefore lacks the central sigma-
point.

###  III-B Gauss–Hermite Rule

From non-singularity of the Vandermonde matrix
$[\bm{\mathbf{V}}]_{nm}=x_{n}^{m-1}$ for any distinct sigma-points
$x_{1},\ldots,x_{p}\in\mathbb{R}$ it follows that there are unique weights
such that $\sum_{n=1}^{p}w_{n}x_{n}^{m}=\int x^{m}\mathrm{N}(0,1)\mathrm{d}x$
for every $m\leq p-1$ (i.e., the rule has a degree of exactness $p-1$).
However, degree of exactness $2p-1$ can be achieved with $p$ sigma-points if
these are selected to be the roots of the $p$-th degree Hermite polynomial
$\mathrm{H}_{p}$. The weights are then given by

$w_{n}=\frac{p!}{p^{2}\mathrm{H}_{p-1}({\xi}^{(n)})^{2}}.$ (18)

This is the Gauss–Hermite (GH) rule [26, 27, 28, 29]. In multivariate
versions, the sigma-points are formed as Cartesian products of the
aforementioned one-dimensional points and the weights are products of $w_{n}$
in eq. 18. The multivariate GH rule exactly integrates functions in the space

$\Pi_{2p-1}^{\text{max}}\triangleq\operatorname{\mathrm{span}}\Bqty\Big{{\bm{\mathbf{x}}}^{\bm{\mathbf{\alpha}}}\,\colon\,\bm{\mathbf{\alpha}}\in\mathbb{N}_{0}^{D},\,\max_{d=1,\ldots,D}\alpha_{d}\leq
2p-1},$ (19)

where
${\bm{\mathbf{x}}}^{\bm{\mathbf{\alpha}}}=\prod_{d=1}^{D}{x}_{d}^{\alpha_{d}}$
denotes multivariate monomial. Because of the Cartesian product design, the
number of points, $N=p^{D}$, in the GH rule grows exponentially with
dimension, which makes it practically unattractive for $D>5$ [27]. The problem
can be partially mitigated by using sparse grids [30].

##  IV Bayesian Quadrature

This section reviews the underlying philosophy of the Bayesian quadrature as
an alternative perspective on numerical integration and describes the
Bayes–Sard quadrature as a necessary stepping stone on the way to building the
Bayes–Sard moment transform proposed in Section V. A general formulation of
the BQ is presented for integrals

$\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{g^{\dagger}({\bm{\mathbf{x}}})}=\int\limits\\!g^{\dagger}({\bm{\mathbf{x}}})p({\bm{\mathbf{x}}})\,\mathrm{d}{\bm{\mathbf{x}}}$
(20)

with arbitrary density function $p$. Vector-valued integrands are discussed in
Section IV-B. The moment transform proposed in Section V then specialises to
the case
$p({\bm{\mathbf{x}}})=\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}$.
Throughout this section, the true integrand will be denoted by $g^{\dagger}$
to distinguish it from the stochastic model of the integrand.

From eq. 9 it is clear that the quadrature approximation of the integral 8 is
based on limited knowledge about the behaviour of the integrand, because it
only relies on finitely many evaluations. The design of classical quadrature
rules typically involves formation of polynomial interpolant passing through
the observed function values, which is then integrated instead of the
intractable integrand. The polynomial interpolation of the integrand
consequently implies that the classical rules are only able to integrate
polynomial integrands exactly. Another downside of the classical rules is that
they are unable to account for the functional uncertainty (interpolation
error), which occurs when the integrand is not a polynomial.

The Bayesian approach to quadrature [11, 12, 15] (see [31, 32, 33] for
equivalent non-probabilistic formulations) aims to address these limitations
by treating the numerical approximation of intractable integrals as a problem
of Bayesian statistical inference, where a prior for the integrand is
specified by a stochastic process model $g({\bm{\mathbf{x}}})$ with user-
defined mean function
$m({\bm{\mathbf{x}}})=\mathbb{E}_{g}\bqty{g({\bm{\mathbf{x}}})}$ and
covariance (or kernel) function
$k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})=\mathbb{C}_{g}\bqty{g({\bm{\mathbf{x}}}),g({\bm{\mathbf{x}}}^{\prime})}$,
where ${\bm{\mathbf{x}}}^{\prime}$ denotes the second argument (not a
transpose). The dataset
$\mathcal{D}=\Bqty{({\bm{\mathbf{x}}}^{(n)},g^{\dagger}({\bm{\mathbf{x}}}^{(n)}))}_{n=1}^{N}$
comprises evaluations of the integrand $g^{\dagger}({\bm{\mathbf{x}}}^{(n)})$
at pre-defined points ${\bm{\mathbf{x}}}^{(n)}$. Conditioning on $\mathcal{D}$
leads to a posterior stochastic process, with mean
$m_{\mathcal{D}}({\bm{\mathbf{x}}})=\mathbb{E}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}})}$
and covariance
$k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})=\mathbb{C}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}}),g({\bm{\mathbf{x}}}^{\prime})}$,
which in turn induces a posterior marginal distribution on the value of the
integral $\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{g({\bm{\mathbf{x}}})}$, with the
first two moments given by [34]

$\displaystyle\mathbb{E}_{g\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{g({\bm{\mathbf{x}}})}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\mathbb{E}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}})}},$
(21)
$\displaystyle\mathbb{V}_{g\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{g({\bm{\mathbf{x}}})}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime}}\bqty{\mathbb{C}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}}),g({\bm{\mathbf{x}}}^{\prime})}}.$
(22)

The mean is a convenient point estimate while the full posterior serves as a
probabilistic model of the integration error. The most common stochastic
process model of the integrand is a Gaussian process (GP), which has been
studied extensively [15, 35]. The standard formulation of the GP regression
model, which is limited to modelling of scalar functions
$g^{\dagger}:\mathbb{R}^{D}\to\mathbb{R}$, has been extended to vector
functions in several ways [36, 37, 38]. Alternatively, one could use the
$t$-process (TP) regression [39, 40], which offers more flexible uncertainty
modelling capabilities and contains the GP as a special case.

In the next section, we specify the Bayes–Sard GP model which is later used to
construct the Bayes–Sard quadrature in Section IV-C.

###  IV-A Bayes–Sard Gaussian Process Model

Let $\pi$ be a linear function space spanned by $Q\leq N$ functions
$\phi_{1},\,\ldots,\,\phi_{Q}\colon\mathbb{R}^{D}\to\mathbb{R}$. Modeling of
the scalar integrand $g^{\dagger}:\mathbb{R}^{D}\to\mathbb{R}$ in Bayes–Sard
quadrature (BSQ) begins by considering a hierarchical GP prior given by

$\displaystyle\bm{\mathbf{\gamma}}$
$\displaystyle\sim\mathrm{N}(\bm{\mathbf{0}},\bm{\mathbf{\Sigma}}_{\pi}),$
(23) $\displaystyle m({\bm{\mathbf{x}}})$
$\displaystyle=\sum_{q=1}^{Q}\gamma_{q}\phi_{q}({\bm{\mathbf{x}}}),$ (24)
$\displaystyle g({\bm{\mathbf{x}}})$
$\displaystyle\sim\mathrm{GP}(m({\bm{\mathbf{x}}}),\,k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}})),$
(25)

where the prior mean function
$m({\bm{\mathbf{x}}}):\mathbb{R}^{D}\to\mathbb{R}$ is composed of basis
functions $\phi_{q}({\bm{\mathbf{x}}})$ of $Q$-dimensional linear space $\pi$
and the prior covariance function (kernel)
$k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}}):\mathbb{R}^{D}\times\mathbb{R}^{D}\to\mathbb{R}$
can be any symmetric positive-definite function parametrized by the vector
$\bm{\mathbf{\theta}}$ (see section IV-E for concrete example). The dependence
on $\bm{\mathbf{\theta}}$ will be tacitly assumed and explicitly denoted only
when required. Discussion about the particular choice of the kernel and its
effects is postponed to Section IV-E. The above model differs from the one
often used in Gaussian process based Bayesian quadrature in that the prior
mean function is non-zero and, furthermore, its coefficients are random.

The next phase in modelling is to consider a flat prior limit on the mean
function coefficients, such that $\bm{\mathbf{\Sigma}}_{\pi}\to\infty$ [41,
Chapter 4]. In order for the GP posterior to be well-defined (i.e., to
guarantee the inverse matrices in sections IV-A to 29 exist), the set
$\mathcal{X}=\\{{\bm{\mathbf{x}}}^{(1)},\,\ldots,\,{\bm{\mathbf{x}}}^{(N)}\\}$
of sigma-points must meet the following condition of $\pi$-unisolvency, which
is related to existence of interpolants formed out of linear combinations of
$\phi_{1},\ldots,\phi_{q}$. For instance, for any $N$ distinct points in one
dimension and any function defined on these points it is possible to construct
a unique polynomial interpolant of degree at most $N-1$. Unfortunately, in
higher dimensions there are point sets for which this does not hold; hence the
need for a unisolvency condition. See [20, Supplement B] for more details.

######  Definition 1 ($\pi$-unisolvency).

Let $\pi$ be a $Q$-dimensional linear space spanned by
$\Bqty{\phi_{1},\ldots,\phi_{q}}$. A point set $\mathcal{X}$ is said to be
$\pi$-unisolvent if and only if the $N\times Q$ alternant matrix
$\bqty{\bm{\mathbf{\Phi}}}_{qn}\triangleq\phi_{q}({\bm{\mathbf{x}}}^{(n)})$ is
of full-rank.

Following the Bayesian paradigm, the final step is to condition the GP on the
set of sigma-points $\mathcal{X}$ and the corresponding integrand evaluations,
collectively denoted as
$\mathcal{D}=\Bqty{({\bm{\mathbf{x}}}^{(n)},g^{\dagger}({\bm{\mathbf{x}}}^{(n)}))}_{n=1}^{N}$,
to arrive at the GP posterior. The GP posterior in the flat prior limit, for
$\pi$-unisolvent $\mathcal{X}$, becomes [20]

$\displaystyle\mathbb{E}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}})}$
$\displaystyle=\bm{\mathbf{k}}({\bm{\mathbf{x}}})^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{y}}$
$\displaystyle\mathrel{\phantom{=}}-\bm{\mathbf{\psi}}({\bm{\mathbf{x}}})^{\top}\bqty\big{\bm{\mathbf{\Phi}}^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{\Phi}}}^{-1}\bm{\mathbf{\Phi}}^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{y}},$
(26)
$\displaystyle\mathbb{C}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}}),g({\bm{\mathbf{x}}}^{\prime})}$
$\displaystyle=k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})-\bm{\mathbf{k}}({\bm{\mathbf{x}}})^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{k}}({\bm{\mathbf{x}}}^{\prime})$
$\displaystyle\mathrel{\phantom{=}}+\bm{\mathbf{\psi}}({\bm{\mathbf{x}}})^{\top}\bqty\big{\bm{\mathbf{\Phi}}^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{\Phi}}}^{-1}\bm{\mathbf{\psi}}({\bm{\mathbf{x}}}^{\prime}),$
(27)

where
$\bm{\mathbf{\psi}}({\bm{\mathbf{x}}})=\bm{\mathbf{k}}({\bm{\mathbf{x}}})\bm{\mathbf{K}}^{-1}\bm{\mathbf{\Phi}}-\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})$,
$\bqty{\bm{\mathbf{k}}({\bm{\mathbf{x}}})missing}_{n}\triangleq
k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{(n)})$,
$\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})missing}_{q}\triangleq\phi_{q}({\bm{\mathbf{x}}})$
and $\bqty{\bm{\mathbf{y}}}_{n}\triangleq
g^{\dagger}({\bm{\mathbf{x}}}^{(n)})$, where $\bqty{\,\cdot\,}_{n}$ denotes
the $n$-th element of the given vector. It is worth pointing out that all
sigma-point sets in the established classical filters are $\pi$-unisolvent.

We further restrict the model to the case when $N=Q$, which means the
alternant matrix $\bm{\mathbf{\Phi}}$ is square and, due to $\pi$-unisolvency
of $\mathcal{X}$, invertible. This leads to the final form of the posterior
mean and covariance of the Bayes–Sard GP model:

$\displaystyle\mathbb{E}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}})}$
$\displaystyle=\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{y}},$
(28)
$\displaystyle\mathbb{C}_{g\mid\mathcal{D}}\bqty{g({\bm{\mathbf{x}}}),g({\bm{\mathbf{x}}}^{\prime})}$
$\displaystyle=k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})-2\bm{\mathbf{k}}({\bm{\mathbf{x}}})^{\top}\bm{\mathbf{\Phi}}^{-\top}\bm{\mathbf{\phi}}({\bm{\mathbf{x}}}^{\prime})$
$\displaystyle\mathrel{\phantom{=}}+\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}\bqty\big{\bm{\mathbf{\Phi}}^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{\Phi}}}^{-1}\bm{\mathbf{\phi}}({\bm{\mathbf{x}}}^{\prime}).$
(29)

Note that the posterior mean now only depends on the choice of the function
space $\pi$ and the kernel affects only the posterior covariance.

###  IV-B Vector-Valued Integrands

Until now, we have only considered scalar-valued integrands. The model
specified by eqs. 28 and 29 can be straightforwardly extended to vector-valued
integrands $\bm{\mathbf{g}}^{\dagger}:\mathbb{R}^{D}\to\mathbb{R}^{E}$ that
comply with the specification of the moment transformation problem in eq. 10.
Noticing that we can decompose the integrand as
$\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}})=\bmqty{undef}^{\top}$, the
simplest solution is to use eqs. 28 and 29 to model each coordinate function
independently, either using a common kernel parameter for all outputs, which
is accomplished by

$g_{e}({\bm{\mathbf{x}}})\mid\mathcal{D}\sim\mathrm{GP}\pqty{m_{\mathcal{D}}({\bm{\mathbf{x}}}),\,k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}})},$
(30)

or using a different kernel parameter values for each output, so that

$g_{e}({\bm{\mathbf{x}}})\mid\mathcal{D}\sim\mathrm{GP}\pqty{m_{\mathcal{D}}({\bm{\mathbf{x}}}),\,k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}}_{e})}$
(31)

for all $e=1,\,\ldots,\,E$. In both cases, the GP posterior mean function is
given as

$\bm{\mathbf{m}}_{\mathcal{D}}({\bm{\mathbf{x}}})\triangleq\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{\phi}}({\bm{\mathbf{x}}}),$
(32)

where $\bqty{\bm{\mathbf{Y}}}_{ne}=g^{\dagger}_{e}({\bm{\mathbf{x}}}^{(n)})$.
For the single-parameter model 30 the posterior covariance becomes

$\bm{\mathbf{K}}_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})\triangleq\mathbb{C}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\bm{\mathbf{g}}({\bm{\mathbf{x}}}^{\prime})}=k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}})\,\bm{\mathbf{I}}_{E}$
(33)

and for the multi-parameter model 31, we get

$\vphantom{\Big{(}}\bm{\mathbf{K}}_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})=\mathrm{diag}\pqty{\,\bqty\big{\matrixquantity{k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}}_{1})&\ldots&k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}}_{E})}}\,}.$
(34)

Both of these modelling choices assume that the outputs are conditionally
independent given the inputs. Alternatively, the use of multi-output Gaussian
processes [37, 38] would make it possible to model correlations between
coordinate functions and use coordinate-dependent sigma-points at the expense
of increased computational cost.

###  IV-C Bayes–Sard Quadrature

The advantage of using a GP for modelling the integrand is that as it gets
transformed by the integral, which is a linear operator, the resulting
distribution over the value of the integral is also Gaussian333Analogous
invariance result holds for the TP model as well.. The Bayes–Sard quadrature
[20, 42, 12] enables enforcing exactness conditions of the form

$\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}=\int\limits\\!\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}})p({\bm{\mathbf{x}}})\,\mathrm{d}{\bm{\mathbf{x}}}$

for all functions
$\bm{\mathbf{g}}^{\dagger}\colon\mathbb{R}^{D}\to\mathbb{R}^{E}$ such that
$g^{\dagger}_{e}\in\pi$ for each $e=1,\ldots,E$. As shown in Section IV-D, the
classical quadrature methods can be replicated by judicious choice of the
function space $\pi$. The posterior integral mean and variance under the
Bayes–Sard quadrature are straightforwardly derived by plugging the Bayes–Sard
GP model moments from eqs. 32 and 34 into the general BQ expressions in eqs.
21 and 22.

For the mean of the posterior distribution of the integral, we have

$\displaystyle\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{m}}_{\mathcal{D}}({\bm{\mathbf{x}}})}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{\Phi}}^{-1}\bar{\bm{\mathbf{\phi}}},$
(35)

where

$\bqty{\bar{\bm{\mathbf{\phi}}}}_{q}=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\phi_{q}({\bm{\mathbf{x}}})}=\int\limits\\!\phi_{q}({\bm{\mathbf{x}}})p({\bm{\mathbf{x}}})\,\mathrm{d}{\bm{\mathbf{x}}}.$
(36)

Recognizing that the vector of quadrature weights is
$\bm{\mathbf{w}}=\bm{\mathbf{\Phi}}^{-1}\bar{\bm{\mathbf{\phi}}}$, we see that
the posterior mean of the integral

$\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{w}}=\sum_{n=1}^{N}w_{n}\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}}^{(n)}),$
(37)

takes on the form of weighted sum from eq. 9. The integral covariance becomes

$\vphantom{\Bigg{(}}\mathbb{V}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}=\mathbb{E}_{{\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime}}\bqty{\bm{\mathbf{K}}_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})}=\mathrm{diag}\pqty{\,\bmqty{undef}\,}$
(38)

where

$\displaystyle\bar{k}^{e}_{\mathcal{D}}$
$\displaystyle\triangleq\mathbb{E}_{{\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime}}\bqty{k_{\mathcal{D}}({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime};\bm{\mathbf{\theta}}_{e})}$
$\displaystyle=\bar{k}-2\bar{\bm{\mathbf{k}}}^{\top}\bm{\mathbf{\Phi}}^{-\top}\bar{\bm{\mathbf{\phi}}}+\bar{\bm{\mathbf{\phi}}}^{\top}\bqty\big{\bm{\mathbf{\Phi}}^{\top}\bm{\mathbf{K}}^{-1}\bm{\mathbf{\Phi}}}^{-1}\bar{\bm{\mathbf{\phi}}}.$
(39)

Since the single-parameter model in eq. 33 is a special case of eq. 34, the
posterior integral variance under this model would be a trivial modification
of eq. 38.

###  IV-D Relationship to Classical Sigma-Point Rules

As stated in the previous section, careful selection of $\pi$ (via the basis
functions $\phi_{q}$) allows for recovery of many well-known classical
quadrature rules used in nonlinear filtering. Below, we show that the
unscented transform and the Gauss–Hermite rule are special cases of the BSQ
whenever the space $\pi$ is spanned by suitably selected polynomial basis.
Similar results can be proved for many other sigma-point rules. Note that the
BSQ reports a non-zero integral variance even for $\bm{\mathbf{g}}^{\dagger}$
whose coordinate functions are in $\pi$ (and hence integrated exactly). This
behaviour is desirable because, given only a finite set of function values,
one can never tell with certainty the true nature of the integrand.

######  Theorem 1.

Consider the standard Gaussian distribution,
$p({\bm{\mathbf{x}}})=\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}$.
Select the $2D+1$ dimensional function space

$\pi=\operatorname{\mathrm{span}}\Bqty{1,\,x_{1},\,\ldots,\,x_{D},\,x_{1}^{2},\,\ldots,\,x_{D}^{2}}$
(40)

and the $N=2D+1$ unscented transform points 16. Then, the Bayes–Sard weights
$\bm{\mathbf{w}}=\bm{\mathbf{\Phi}}^{-1}\bar{\bm{\mathbf{\phi}}}$ that
determine the posterior mean 35 coincide with the unscented transform weights
17.

###### Proof:

Because $\dim(\pi)=N$, the Bayes–Sard weights $\bm{\mathbf{w}}$ solve the
linear system $\bm{\mathbf{\Phi}}\bm{\mathbf{w}}=\bar{\bm{\mathbf{\phi}}}$.
That is, they are the unique weights such that

$\sum_{n=0}^{2D}w_{n}v({\bm{\mathbf{\xi}}}^{(n)})=\int\limits\\!v({\bm{\mathbf{x}}})\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}\,\mathrm{d}{\bm{\mathbf{x}}}$
(41)

for every polynomial $v\in\pi$. In the following, let $d=1,\,\ldots,\,D$. We
have
$\int\limits\\!\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}\,\mathrm{d}{\bm{\mathbf{x}}}=1$
and

$\int\limits\\!x_{d}\,\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}\,\mathrm{d}{\bm{\mathbf{x}}}=0,\quad\int\limits\\!x_{d}^{2}\,\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}\,\mathrm{d}{\bm{\mathbf{x}}}=1.$
(42)

Consequently, eq. 41 is equivalent to

$\sum_{n=0}^{2D}w_{n}=1,\quad\sum_{n=0}^{2D}w_{n}x_{n,d}=0,\quad\sum_{n=0}^{2D}w_{n}x_{n,d}^{2}=1$
(43)

Because $\bm{\mathbf{\xi}}_{d}=-\bm{\mathbf{\xi}}_{D+d}$, the second of these
equations implies that $w_{d}=w_{D+d}$, while the third one yields
$w_{d}=w_{D+d}=1/(2c)$. Furthermore, $w_{0}=\kappa/c$ due to the weights
summing up to one. We have thus solved the BSQ weights
$\bm{\mathbf{w}}=\bm{\mathbf{\Phi}}^{-1}\bar{\bm{\mathbf{\phi}}}$ and see that
they are precisely UT weights in eq. 17. ∎

######  Theorem 2.

Consider the standard Gaussian distribution,
$p({\bm{\mathbf{x}}})=\mathrm{N}\pqty{\left.\bm{\mathbf{x}}\,\middle|\,\bm{\mathbf{0}},\,\bm{\mathbf{I}}\right.}$,
and let $p\geq 1$. Select the $p^{D}$ dimensional function space
$\pi=\Pi_{p-1}^{\text{max}}\triangleq\operatorname{\mathrm{span}}\Bqty{{\bm{\mathbf{x}}}^{\bm{\mathbf{\alpha}}}\,\colon\,\bm{\mathbf{\alpha}}\in\mathbb{N}_{0}^{D},\,\max_{d=1,\ldots,D}\alpha_{d}\leq
p-1}$, and the points that constitute the Cartesian product of the roots of
the $p$-th degree Hermite polynomial. Then, the Bayes–Sard weights
$\bm{\mathbf{w}}=\bm{\mathbf{\Phi}}^{-1}\bar{\bm{\mathbf{\phi}}}$ that
determine the posterior mean 35 coincide with the classical Gauss–Hermite
weights from Section III-B.

###### Proof:

Since the Bayes–Sard weights yield, by their definition, a quadrature rule
exact for functions in $\pi$ and it is known that, given the Gauss–Hermite
points, the Gauss–Hermite weights are the unique weights that determine a
quadrature rule exact for this very same function space (see Section III-B),
the result follows. ∎

###  IV-E Choice of Kernel

As already noted, the posterior mean for the integral produced by the BSQ
depends only on $\pi$ and the kernel controls the posterior variance of the
integral. The reasonableness of the BSQ output depends on the reasonableness
of the assumption that $g^{\dagger}$ is “well modelled” by the GP specified by
the kernel $k$. Consequently, selection of the kernel is important in order to
ensure that the integral variance is meaningful in modelling the integration
error. At the same time, the functional form of the kernel is constrained by
the requirement in BSQ to analytically compute the integral of the kernel. To
facilitate analytic tractability of the Bayes–Sard moment transform,
introduced next, we use the radial basis function (RBF) kernel

$k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{\prime})=\alpha^{2}\prod_{d=1}^{D}\exp\pqty{-\frac{({x}_{d}-{x}_{d}^{\prime})^{2}}{2\ell^{2}_{d}}}$
(44)

throughout the remainder. The parameters $\bm{\mathbf{\theta}}$ of this kernel
consist of the scale parameter $\alpha>0$ and dimension-wise lengthscale
parameters $\ell_{1},\ldots,\ell_{D}>0$. A particular modelling assumption
associated with this kernel is that the integrand is infinitely
differentiable. If this is not the case (i.e., there is model
misspecification) the proposed method still works but the uncertainty
quantification for the integral may be rendered less meaningful. For certain
classes of kernels it has been shown that convergence rates to the true
integral as $N\to\infty$ are not much affected by model misspecification [43,
44].

##  V Bayes–Sard Moment Transform

The simplest way to design a moment transform is to use the BSQ directly for
approximation of the moment integrals in eqs. 15a to 15c. However, this design
does not reflect integral uncertainty, which is the key advantage of Bayesian
quadrature, not to mention the fact that we would only obtain the classical
rules as a result. To resolve this issue, we employ the same general
conceptual framework used in the design of the GPQ moment transform in [18],
which can account for the variance of the mean integral 15a.

###  V-A Incorporating Integration Error

First, it is important to realize that the output variable $\bm{\mathbf{y}}$
is now subject to an additional source of uncertainty in $\bm{\mathbf{g}}$
introduced by the model. The key idea is to account for all sources of
uncertainty in the computed moments, which can be achieved with the following

$\displaystyle\bm{\mathbf{\mu}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}})}$
$\displaystyle\approx$
$\displaystyle\hat{\bm{\mathbf{\mu}}}=\mathbb{E}_{{\bm{\mathbf{x}}},\,\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}$
(45) $\displaystyle\bm{\mathbf{\Pi}}$
$\displaystyle=\mathbb{C}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}}),\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}})}$
$\displaystyle\approx$
$\displaystyle\hat{\bm{\mathbf{\Pi}}}=\mathbb{C}_{{\bm{\mathbf{x}}},\,\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\,\bm{\mathbf{g}}({\bm{\mathbf{x}}})}$
(46) $\displaystyle\bm{\mathbf{C}}$
$\displaystyle=\mathbb{C}_{{\bm{\mathbf{x}}}}\bqty{{\bm{\mathbf{x}}},\bm{\mathbf{g}}^{\dagger}({\bm{\mathbf{x}}})}$
$\displaystyle\approx$
$\displaystyle\hat{\bm{\mathbf{C}}}=\mathbb{C}_{{\bm{\mathbf{x}}},\,\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{{\bm{\mathbf{x}}},\,\bm{\mathbf{g}}({\bm{\mathbf{x}}})}$
(47)

Using the law of total expectation and covariance, the approximate moments of
the output can be written as

$\displaystyle\hat{\bm{\mathbf{\mu}}}$
$\displaystyle=\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}},$
(48) $\displaystyle\hat{\bm{\mathbf{\Pi}}}$
$\displaystyle=\mathbb{C}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}+\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{C}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\bm{\mathbf{g}}({\bm{\mathbf{x}}})}},$
(49)
$\displaystyle=\mathbb{C}_{{\bm{\mathbf{x}}}}\bqty{\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}+\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\mathbb{C}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\bm{\mathbf{g}}({\bm{\mathbf{x}}})}},$
(50) $\displaystyle\hat{\bm{\mathbf{C}}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{{\bm{\mathbf{x}}}\,\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}-\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{{\bm{\mathbf{x}}}missing}\mathbb{E}_{\bm{\mathbf{g}}\mid\mathcal{D},{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}.$
(51)

The first equality exposes the fact that integral mean is obtained by
integrating the mean function of the integrand model. The way the integral
uncertainty is incorporated into the output covariance is revealed by eq. 49.
Note that since the model of the integrand has conditionally independent
outputs, the covariance of the integral,
$\mathbb{C}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\mathbb{E}_{{\bm{\mathbf{x}}}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}})}}$,
and the model covariance,
$\mathbb{C}_{\bm{\mathbf{g}}\mid\mathcal{D}}\bqty{\bm{\mathbf{g}}({\bm{\mathbf{x}}}),\bm{\mathbf{g}}({\bm{\mathbf{x}}})}$,
are diagonal matrices. When either of the covariances approaches zero, eqs. 48
to 51 approach their true values. From now on, we will work with the output
covariance in the form 50 because it is easier to analyse and implement.

###  V-B Derivation of Transformed Moments

In the following derivations, explicit conditioning on $\mathcal{D}$ in the
expectations is omitted to reduce notational clutter. We also assume that the
stochastic decoupling substitution has taken place in the integrals, so that
$\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})=\bm{\mathbf{g}}(\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}})$.

Taking the expression for the mean function of the model in eq. 32 and
plugging it into eq. 48, the output mean of the Bayes–Sard quadrature moment
transform (BSQMT) becomes

$\hat{\bm{\mathbf{\mu}}}=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\mathbb{E}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{\Phi}}^{-\top}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{w}},$
(52)

where
$\bm{\mathbf{w}}=\bm{\mathbf{\Phi}}^{-\top}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}$
are the mean weights. The output covariance becomes

$\displaystyle\hat{\bm{\mathbf{\Pi}}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty\big{\mathbb{E}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}\,\mathbb{E}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}^{\top}}-\hat{\bm{\mathbf{\mu}}}\hat{\bm{\mathbf{\mu}}}^{\top}+\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\mathbb{C}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}}),\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}}$
$\displaystyle=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{\Phi}}^{-\top}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{Y}}-\hat{\bm{\mathbf{\mu}}}\hat{\bm{\mathbf{\mu}}}^{\top}+\bar{\sigma}^{2}\bm{\mathbf{I}}_{E}$
$\displaystyle=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{W}}\bm{\mathbf{Y}}-\hat{\bm{\mathbf{\mu}}}\hat{\bm{\mathbf{\mu}}}^{\top}+\bar{\sigma}^{2}\bm{\mathbf{I}}_{E}$
(53)

where the expected model variance is

$\bar{\sigma}^{2}=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{k({\bm{\mathbf{\xi}}},\,{\bm{\mathbf{\xi}}})}-\tr\bqty\big{\bm{\mathbf{D}}^{\top}\bm{\mathbf{\Phi}}^{-\top}+\bm{\mathbf{D}}\bm{\mathbf{\Phi}}^{-1}-\bm{\mathbf{W}}\bm{\mathbf{K}}}.$
(54)

Here
$\bm{\mathbf{D}}=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{k}}({\bm{\mathbf{\xi}}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}$
and the covariance weights are
$\bm{\mathbf{W}}=\bm{\mathbf{\Phi}}^{-\top}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}\bm{\mathbf{\Phi}}^{-1}$.
Finally, the covariance between the input and output becomes

$\displaystyle\hat{\bm{\mathbf{C}}}$
$\displaystyle=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{(\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}})\,\mathbb{E}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}}-\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{m}}+\bm{\mathbf{L}}{\bm{\mathbf{\xi}}}missing}\,\mathbb{E}_{\bm{\mathbf{g}},{\bm{\mathbf{\xi}}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}$
$\displaystyle=\bm{\mathbf{m}}\hat{\bm{\mathbf{\mu}}}+\bm{\mathbf{L}}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{{\bm{\mathbf{\xi}}}\,\mathbb{E}_{\bm{\mathbf{g}}}\bqty{\tilde{\bm{\mathbf{g}}}({\bm{\mathbf{\xi}}})}}-\bm{\mathbf{m}}\hat{\bm{\mathbf{\mu}}}$
$\displaystyle=\bm{\mathbf{L}}\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{{\bm{\mathbf{\xi}}}\,\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{Y}}=\bm{\mathbf{L}}\bm{\mathbf{W}}_{c}\bm{\mathbf{Y}}$
(55)

where the cross-covariance weights are
$\bm{\mathbf{W}}_{c}=\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{{\bm{\mathbf{\xi}}}\,\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}\bm{\mathbf{\Phi}}^{-1}$.

It has now become evident that the output moments depend on the expectations
of the basis functions. In Section IV, we have shown that the classical moment
transforms can be recovered when the basis functions are multivariate
polynomials. When this basis and the RBF kernel eq. 44 are used, the
expectations above are available in closed form. Derivations are confined to
Sections -A to -C so as not to disrupt the flow. The complete algorithm of the
Bayes–Sard moment transform is summarized in Alg. 1.

Input: The mean $\bm{\mathbf{m}}$ and the covariance $\bm{\mathbf{P}}$ of the
input variable ${\bm{\mathbf{x}}}$, the integrand
$\bm{\mathbf{g}}({\bm{\mathbf{x}}})$, the matrix of unit sigma-points
$\bm{\mathbf{\Xi}}$ and the kernel parameters $\bm{\mathbf{\theta}}$.

Output: Approximate mean $\hat{\bm{\mathbf{\mu}}}$ and covariance
$\hat{\bm{\mathbf{\Pi}}}$ of the output variable
$\bm{\mathbf{y}}=\bm{\mathbf{g}}({\bm{\mathbf{x}}})$, and approximate input-
output covariance $\hat{\bm{\mathbf{C}}}$.

Function  _BSQMT(_$\bm{\mathbf{g}}({\bm{\mathbf{x}}})$ , $\bm{\mathbf{m}}$,
$\bm{\mathbf{P}}$, $\bm{\mathbf{\Xi}}$, $\bm{\mathbf{\theta}}$_)_

// form sigma-points

1 $\bm{\mathbf{L}}\leftarrow$ MatrixFactor(_$\bm{\mathbf{P}}$_)

2 $\bm{\mathbf{X}}\leftarrow\bm{\mathbf{m}}+\bm{\mathbf{L}}\bm{\mathbf{\Xi}}$

3

4
$\bar{k}\leftarrow\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{k({\bm{\mathbf{\xi}}},{\bm{\mathbf{\xi}}}\,;\bm{\mathbf{\theta}})}$

5
$\bar{\bm{\mathbf{\phi}}}\leftarrow\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}$

6
$\bm{\mathbf{A}}\leftarrow\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}$

7
$\bm{\mathbf{B}}\leftarrow\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{{\bm{\mathbf{\xi}}}\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}$

8 for  _$n\leftarrow 1$ to $N$_ do

9
$\bqty{\bm{\mathbf{Y}}}_{n*}\leftarrow\bm{\mathbf{g}}({\bm{\mathbf{x}}}^{(n)})$

10 $\bqty{\bm{\mathbf{\Phi}}}_{n*}\leftarrow\bmqty{undef}$

11
$\bqty{\bm{\mathbf{D}}}_{n*}\leftarrow\mathbb{E}_{{\bm{\mathbf{\xi}}}}\bqty{k({\bm{\mathbf{\xi}}},{\bm{\mathbf{\xi}}}^{(n)}\,;\bm{\mathbf{\theta}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}$

12 for  _$m\leftarrow 1$ to $N$_ do

13 $\bqty{\bm{\mathbf{K}}}_{nm}\leftarrow
k({\bm{\mathbf{\xi}}}^{(n)},\,{\bm{\mathbf{\xi}}}^{(m)};\,\bm{\mathbf{\theta}})$

14

15 end for

16

17 end for

// compute BSQ weights

18
$\bm{\mathbf{w}}\leftarrow\bm{\mathbf{\Phi}}^{-\top}\bar{\bm{\mathbf{\phi}}}$

19
$\bm{\mathbf{W}}\leftarrow\bm{\mathbf{\Phi}}^{-\top}\bm{\mathbf{A}}\bm{\mathbf{\Phi}}^{-1}$

20 $\bm{\mathbf{W}}_{c}\leftarrow\bm{\mathbf{B}}\bm{\mathbf{\Phi}}^{-1}$

21

// compute transformed moments

22
$\bar{\sigma}^{2}\leftarrow\bar{k}-\tr\bqty\big{\bm{\mathbf{D}}^{\top}\bm{\mathbf{\Phi}}^{-\top}+\bm{\mathbf{D}}\bm{\mathbf{\Phi}}^{-1}-\bm{\mathbf{W}}\bm{\mathbf{K}}}$

23 $\hat{\bm{\mathbf{\mu}}}\leftarrow\bm{\mathbf{Y}}^{\top}\bm{\mathbf{w}}$

24
$\hat{\bm{\mathbf{\Pi}}}\leftarrow\bm{\mathbf{Y}}^{\top}\bm{\mathbf{W}}\bm{\mathbf{Y}}-\hat{\bm{\mathbf{\mu}}}\hat{\bm{\mathbf{\mu}}}^{\top}+\bar{\sigma}^{2}\bm{\mathbf{I}}$

25
$\hat{\bm{\mathbf{C}}}\leftarrow\bm{\mathbf{L}}\bm{\mathbf{W}}_{c}\bm{\mathbf{Y}}$

26

27 return $\hat{\bm{\mathbf{\mu}}}$, $\hat{\bm{\mathbf{\Pi}}}$,
$\hat{\bm{\mathbf{C}}}$

28 end

Algorithm 1 Bayes-Sard quadrature moment transform

######  Theorem 3.

The BSQ output covariance $\hat{\bm{\mathbf{\Pi}}}$ is positive semi-definite.

###### Proof:

Using the expression for the BSQ mean weights from Alg. 1, we can write the
output covariance as
$\hat{\bm{\mathbf{\Pi}}}=\bm{\mathbf{Y}}^{\top}\bm{\mathbf{\Phi}}^{-\top}(\bm{\mathbf{A}}-\bar{\bm{\mathbf{\phi}}}\bar{\bm{\mathbf{\phi}}}^{\top})\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{Y}}+\bar{\sigma}^{2}\bm{\mathbf{I}}_{E}$.
Define $\bm{\mathbf{Z}}=\bm{\mathbf{\Phi}}^{-1}\bm{\mathbf{Y}}$ and
$\tilde{\bm{\mathbf{A}}}=\bm{\mathbf{A}}-\bar{\bm{\mathbf{\phi}}}\bar{\bm{\mathbf{\phi}}}^{\top}$,
then
$\bm{\mathbf{\Pi}}=\bm{\mathbf{Z}}^{\top}\tilde{\bm{\mathbf{A}}}\bm{\mathbf{Z}}+\bar{\sigma}^{2}\bm{\mathbf{I}}_{E}$.
We recognize that
$\vphantom{\Big{(}}\tilde{\bm{\mathbf{A}}}=\mathbb{V}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}=\mathbb{E}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})^{\top}}-\mathbb{E}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}\mathbb{E}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{\xi}}})}^{\top}\succeq
0$, which follows from the properties of covariance matrices. This implies
that $\bm{\mathbf{Z}}^{\top}\tilde{\bm{\mathbf{A}}}\bm{\mathbf{Z}}\succeq 0$
for any matrix $\bm{\mathbf{Z}}$. Because $\bar{\sigma}^{2}\geq 0$, we have
that $\hat{\bm{\mathbf{\Pi}}}\succeq 0$. ∎

###  V-C Relationship to the Gaussian Process Quadrature MT

The recently proposed Gaussian process quadrature moment transform (GPQMT)
[18], together with the BSQMT, are both instances of the general BQ framework.
The GPQMT uses a zero-mean GP prior model of the integrand as opposed to the
more sophisticated hierarchical prior in eqs. 23 and 25. As a result, the
GPQMT weights are affected by the choice of kernel and its parameter values,
which is not the case in the BSQMT, where the kernel only affects the last
term of the transformed covariance and the weights depend only on the sigma-
points and the choice of the function space $\pi$. Consequently, this makes
BSQMT much less sensitive to misspecification of the kernel parameters, which
is a notorious problem plaguing GPQMT. Discussion of the choice of kernel
parameters can be found in the original publication [18].

Compared to the zero-mean GP employed in GPQMT, the Bayes–Sard GP is a
stronger prior, which means it can provide better fit to the integrand when
conditioned on smaller datasets, such as the UT sigma-points, which are
especially attractive in local filtering applications.

###  V-D BSQ Moment Transform in Sigma-Point Filtering

As outlined in Section II, the local filtering algorithms use the moment
transformations for computing the predictive moments of the system state and
measurement. Alg. 2 summarizes the Bayes–Sard quadrature Kalman filter
(BSQKF), which employs the proposed BSQ moment transform for computing the
predictive moments from Table I. The BSQKF takes two different kernel
parameter values, $\bm{\mathbf{\theta}}_{f}$ and $\bm{\mathbf{\theta}}_{h}$,
because there are two different functions that need to be integrated (see eqs.
1 and 2).

Input: Sequence of measurements $\Bqty\big{\bm{\mathbf{z}}_{k}}_{k=1}^{K}$,
initial conditions $\bm{\mathbf{m}}^{x}_{0|0},\,\bm{\mathbf{P}}^{x}_{0|0}$,
kernel parameters $\bm{\mathbf{\theta}}_{f}$ and $\bm{\mathbf{\theta}}_{h}$,
unit sigma-points $\bm{\mathbf{\Xi}}$

Output: Sequence of state estimates and covariances
$\Bqty\big{\bm{\mathbf{m}}^{x}_{k|k},\,\bm{\mathbf{P}}^{x}_{k|k}missing}_{k=1}^{K}$

for  _$k\leftarrow 1$ to $K$_ do

// predictive state moments

$\bm{\mathbf{m}}^{x}_{k|k-1},\ \bm{\mathbf{P}}^{x}_{k|k-1}\leftarrow$

BSQMT(_$\bm{\mathbf{f}}({\bm{\mathbf{x}}}_{k-1})$ ,
$\bm{\mathbf{m}}^{x}_{k-1|k-1},\ \bm{\mathbf{P}}^{x}_{k-1|k-1},\
\bm{\mathbf{\Xi}},\ \bm{\mathbf{\theta}}_{f}$_)

$\bm{\mathbf{P}}^{x}_{k|k-1}\leftarrow\bm{\mathbf{P}}^{x}_{k|k-1}+\bm{\mathbf{Q}}$

// predictive measurement moments

$\bm{\mathbf{m}}^{z}_{k|k-1},\ \bm{\mathbf{P}}^{z}_{k|k-1},\
\bm{\mathbf{P}}^{xz}_{k|k-1}\leftarrow$

BSQMT(_$\bm{\mathbf{h}}({\bm{\mathbf{x}}}_{k})$ ,
$\bm{\mathbf{m}}^{x}_{k|k-1},\ \bm{\mathbf{P}}^{x}_{k|k-1},\
\bm{\mathbf{\Xi}},\ \bm{\mathbf{\theta}}_{h}$_)

$\bm{\mathbf{P}}^{z}_{k|k-1}\leftarrow\bm{\mathbf{P}}^{z}_{k|k-1}+\bm{\mathbf{R}}$

// measurement update (filtering)

$\bm{\mathbf{G}}_{k}\leftarrow\bm{\mathbf{P}}^{xz}_{k|k-1}\pqty\big{\bm{\mathbf{P}}^{z}_{k|k-1}}^{-1}$

$\bm{\mathbf{m}}^{x}_{k|k}\leftarrow\bm{\mathbf{m}}^{x}_{k|k-1}+\bm{\mathbf{G}}_{k}\pqty\big{\bm{\mathbf{z}}_{k}-\bm{\mathbf{m}}^{z}_{k}}^{\top}$

$\bm{\mathbf{P}}^{x}_{k|k}\leftarrow\bm{\mathbf{P}}^{x}_{k|k-1}-\bm{\mathbf{G}}_{k}\bm{\mathbf{P}}^{z}_{k|k-1}\bm{\mathbf{G}}_{k}^{\top}$

end for

Algorithm 2 Bayes-Sard quadrature Kalman filter.

###  V-E Kernel-Agnostic BSQ Moment Transform

The BSQ moment transform stated in Alg. 1 depends on the choice of kernel and
the values of its parameters, which in general GP regression are
conventionally estimated from data by marginal likelihood maximization [35].
Considering the fact that the classical sigma-point sets are inherently sparse
relative to dimensionality of the input space, maximum likelihood estimation
is not expected to work well in this context. The last resort in practice is
therefore manual parameter tuning. When the moment transform is employed on
its own, like in Section VI-A, tuning can be informed by the analytic form of
integrand with the goal of producing good model fit, which typically yields
good results. However, this tactic has not been observed to work well in
filtering applications (see Sections VI-B and VI-C), because good filter
performance may not be necessarily achieved by parameters that maximize model
fit to the integrands $\bm{\mathbf{f}}$ and $\bm{\mathbf{h}}$ from eqs. 1 and
2.

With all these considerations in mind, in filtering applications it is easier
to re-parametrise the BSQMT in terms of the expected model variance (EMV)
$\bar{\sigma}^{2}$ directly rather than through the complicated dependence on
kernel parameters as given by eq. 54. Note that this re-parametrisation
renders the BSQMT independent of the choice of kernel and thus removes the
need for computation of the EMV-related terms in Alg. 1 and Algs. 1 to 1 in
Alg. 1.

A more systematic approach than manual tuning would be to treat
$\bar{\sigma}^{2}$ as an unknown static parameter, which could be estimated
using the energy function minimization [45]. In case $\bar{\sigma}^{2}$ is
assumed to be time-variant, the state augmentation approach [46] might prove
useful. Exploration of the various parameter estimation techniques is beyond
the scope of this article.

##  VI Numerical Experiments

This section contains a static example and two numerical filtering experiments
that demonstrate performance of the Bayes–Sard quadrature moment transform and
the Bayes–Sard quadrature Kalman filter.

###  VI-A Mapping from Polar to Cartesian Coordinates

To gain insight into the performance of our proposed method, we first consider
the moment transform in isolation, outside of the filtering context. As an
example of ubiquitous nonlinearity, we use the conversion from polar to
Cartesian coordinates, which appears in radars and laser range finders and is
given by

$\bmqty{undef}=\bmqty{undef},$ (56)

where $r$ is radius and $\theta$ is azimuth. When setting the RBF kernel
parameters, we exploited the conditional linearity of the mapping, which is
reflected in the setting $\ell_{1}=60$ and $\ell_{2}=6$, where the range
coordinate $r$ is given relatively large value. The kernel scaling was set to
$\alpha=1$.

The BSQ and GPQ transforms with UT points ($\kappa=2$) were compared with the
classical unscented transform. Since the input space is given in polar
coordinates, the input means were placed on a spiral and several covariances
were considered for each input mean. More specifically, 10 different positions
on a spiral were chosen as input means $\bm{\mathbf{m}}_{i}=\bmqty{undef}$ and
for each $\bm{\mathbf{m}}_{i}$ we considered 10 different input covariance
matrices $\bm{\mathbf{P}}_{j}=\mathrm{diag}\pqty{\,\bmqty{undef}\,}$, where
$\sigma_{r}=$0.5\text{\,}\mathrm{m}$$ and
$\sigma_{\theta,j}\in\bqty{$6\text{\,}\mathrm{\SIUnitSymbolDegree}$,\,$36\text{\,}\mathrm{\SIUnitSymbolDegree}$}$
for $j=1,\ldots,10$. Agreement between the approximate moments
$(\hat{\bm{\mathbf{\mu}}},\ \hat{\bm{\mathbf{\Pi}}})$ and the ground-truth
moments $(\bm{\mathbf{\mu}},\ \bm{\mathbf{\Pi}})$, computed by Monte Carlo
method with $10\,000$ samples, was measured by the symmetrized KL-
divergence444The advantage of SKL over KL-divergence is that it is symmetric;
that is, it holds
$\mathbb{SKL}\bqty{\left.\bm{\mathbf{\mu}},\bm{\mathbf{\Pi}}\,\middle\|\,\hat{\bm{\mathbf{\mu}}},\hat{\bm{\mathbf{\Pi}}}\right.}=\mathbb{SKL}\bqty{\left.\hat{\bm{\mathbf{\mu}}},\hat{\bm{\mathbf{\Pi}}}\,\middle\|\,\bm{\mathbf{\mu}},\bm{\mathbf{\Pi}}\right.}$.
of two Gaussian densities given by

$\displaystyle\mathbb{SKL}\bqty{\left.\hat{\bm{\mathbf{\mu}}},\hat{\bm{\mathbf{\Pi}}}\,\middle\|\,\bm{\mathbf{\mu}},\bm{\mathbf{\Pi}}\right.}$
$\displaystyle=\frac{1}{4}\left\\{(\bm{\mathbf{\mu}}-\hat{\bm{\mathbf{\mu}}})^{\top}\bm{\mathbf{\Pi}}^{-1}(\bm{\mathbf{\mu}}-\hat{\bm{\mathbf{\mu}}})+\right.$
$\displaystyle\mathrel{\phantom{=}}\phantom{\frac{}{4}\left\\{\right.}(\hat{\bm{\mathbf{\mu}}}-\bm{\mathbf{\mu}})^{\top}\hat{\bm{\mathbf{\Pi}}}^{-1}(\hat{\bm{\mathbf{\mu}}}-\bm{\mathbf{\mu}})+$
$\displaystyle\mathrel{\phantom{=}}\phantom{\frac{}{4}\left\\{\right.}\left.\tr(\bm{\mathbf{\Pi}}^{-1}\hat{\bm{\mathbf{\Pi}}})+\tr(\hat{\bm{\mathbf{\Pi}}}^{-1}\bm{\mathbf{\Pi}})-2E\right\\},$
(57)

where $E=\dim(\bm{\mathbf{\mu}})$.

The plots in the left pane of Fig. 1 show the SKL score averaged over the
input covariances $\bm{\mathbf{P}}_{j}$, while the right pane shows SKL
averaged over the input means $\bm{\mathbf{m}}_{i}$. The proposed Bayes–Sard
quadrature transform with UT points (BSQ-UT) outperforms the classical UT,
while achieving almost the same performance as the GPQ-UT. What is especially
noteworthy here, is the fact that BSQ-UT computes the transformed mean
$\hat{\bm{\mathbf{\mu}}}$ identically to the UT and the only difference is in
the manner in which the transformed covariance $\hat{\bm{\mathbf{\Pi}}}$ is
computed. The Fig. 1 shows that the BSQ-UT, which can be viewed as a
probabilistic equivalent of the UT, offers improved performance, because it
accounts for the numerical integration error of the UT. For completeness we
also included the GPQ-UT transform, which dominates this example, but does not
have such a clear theoretical connection to the UT.

$\displaystyle 5$$\displaystyle 10$Position index$\displaystyle
0.00$$\displaystyle 0.02$$\displaystyle 0.04$$\displaystyle
0.06$SKL$\displaystyle 10$$\displaystyle 20$$\displaystyle 30$Azimuth STD
[$\displaystyle\circ$]$\displaystyle 0.00$$\displaystyle 0.02$$\displaystyle
0.04$$\displaystyle 0.06$BSQ-UTGPQ-UTUT Figure 1: Symmetrised KL divergence
between the MC baseline Gaussian and the approximate Gaussian computed by the
BSQ, GP and UT moment transforms.

###  VI-B Univariate Non-stationary Growth Model

As the first experiment with sigma-point filters based on the BSQ moment
transform, we consider the univariate non-stationary growth model (UNGM),
which is often used to benchmark particle filters [47, 48]. The system
dynamics and the observation model are given by

$\displaystyle x_{k}$
$\displaystyle=\frac{1}{2}x_{k-1}+\frac{25x_{k-1}}{1+{x}^{2}_{k-1}}+8\cos(1.2\,k)+q_{k-1},$
(58) $\displaystyle z_{k}$ $\displaystyle=\frac{1}{20}{x}^{2}_{k-1}+r_{k},$
(59)

with the state noise $q_{k-1}\sim\mathrm{N}\pqty{0,\,10}$, measurement noise
$r_{k}\sim\mathrm{N}\pqty{0,\,1}$ and initial conditions
$x_{0}=x_{0|0}\sim\mathrm{N}\pqty{0,\,5}$. Kernel scaling used in the BSQ with
for the UT ($\kappa=2$) and the 7-th order GH points was set to $\alpha=3$ and
the lengthscales to $\ell=0.3$ and $\ell=0.4$, respectively. For the 5-th
order GH points the kernel parameters were set to $\alpha=5$ and $\ell=0.6$.

The RMSE was used to measure the tracking performance. The inclination
indication (INC) [49], given by

$\mathrm{INC}=\frac{10}{K}\sum_{k=1}^{K}\log_{10}\frac{(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})^{\top}(\bm{\mathbf{P}}^{x}_{k|k})^{-1}(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})}{(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})^{\top}\bm{\mathbf{\Sigma}}_{k}^{-1}(\bm{\mathbf{x}}_{k}-\bm{\mathbf{m}}^{x}_{k|k})},$
(60)

where $\bm{\mathbf{\Sigma}}_{k}$ is the mean-squared error matrix of the
state, was used to measure the credibility of the estimates. A perfectly
balanced estimate has $\mathrm{INC}=0$. For $\mathrm{INC}>0$, the estimate is
said to be optimistic, which is to say the covariance is smaller than it
should be, while negative values indicate pessimism. Refer to [49] for other
related credibility measures. We simulated the model for $K=500$ time steps
and averaged the performance scores over $100$ simulations. The variance of
the average scores was estimated by bootstrapping. The parentheses in Tables
II and III contain the uncertainty as the least significant digits of 2
standard deviations.

The BSQ filters with classical points were tested against the well-known
sigma-point filters as well as the GPQ filters from [18]. As seen in Table II,
the filters based on the BSQ outperform the classical sigma-point filters in
terms of RMSE. Assuming the GH points are used, BSQKFs can outperform the GPQ
filters as well. In comparison with the classical filters, the proposed BSQ
filters also provide much more balanced estimates as evidenced by the values
of the INC in Table III.

| Classical | GPQ | BSQ
---|---|---|---
UT | 10.81 (0.14) | 10.37 (0.08) | 9.70 (0.12)
GH-5 | 10.03 (0.14) | 9.01 (0.11) | 8.82 (0.08)
GH-7 | 9.74 (0.13) | 8.80 (0.10) | 8.61 (0.09)
Table II: Filter RMSE for the UNGM example. | Classical | GPQ | BSQ
---|---|---|---
UT | 12.17 (0.06) | 4.87 (0.01) | 4.57 (0.03)
GH-5 | 10.33 (0.07) | 5.26 (0.03) | 1.85 (0.02)
GH-7 | 9.27 (0.07) | 4.95 (0.03) | 2.52 (0.03)
Table III: Filter INC for the UNGM example.

###  VI-C Radar Tracking of Reentry Vehicle

Radar tracking of an object entering Earth’s atmosphere is a very challenging
tracking scenario for nonlinear filters. The following model has been used to
demonstrate superiority of the UKF over the EKF in [50].

The dataset of the ground-truth system state trajectories comprised of $100$
MC simulations of the SDE

$\dot{\bm{\mathbf{x}}}(t)=\bmqty{undef}$ (61)

using the Euler-Maruyama method for the duration of $200\text{\,}\mathrm{s}$
with time step $\Delta t=$0.05\text{\,}\mathrm{s}$$. The state
$\bm{\mathbf{x}}(t)=\bmqty{undef}$ consists of Cartesian position coordinates
$(p^{x},p^{y})$ in $\mathrm{k}\mathrm{m}$, velocity $(v^{x},v^{y})$ in
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$ and an aerodynamic parameter
$\theta$. The functions $D(t)$ and $G(t)$ are the drag-related and the
gravity-related force terms, given by

$\displaystyle D(t)$
$\displaystyle=\beta(t)\exp\pqty{\frac{R_{0}-R(t)}{H_{0}}}V(t),$ (62)
$\displaystyle G(t)$ $\displaystyle=\frac{Gm_{0}}{r^{3}(t)},$ (63)
$\displaystyle\beta(t)$ $\displaystyle=\beta_{0}\exp(\theta(t)),$ (64)

with $R(t)=\sqrt{(p^{x}(t))^{2}+(p^{y}(t))^{2}}$ being the distance from the
centre of the Earth and $V(t)=\sqrt{(v^{x}(t))^{2}+(v^{y}(t))^{2}}$ the speed.
The parameter values were set identically to [50]; thus $R_{0}=6374$,
$H_{0}=13.406$, $\beta_{0}=$-0.597\,83$$ and
$Gm_{0}=$3.9860\text{\times}{10}^{5}$$. The system initial condition
$\bm{\mathbf{x}}(0)\sim\mathrm{N}\pqty{\bm{\mathbf{m}}^{x}_{0},\,\bm{\mathbf{P}}^{x}_{0}}$
was set to

$\displaystyle\bm{\mathbf{m}}^{x}_{0}$ $\displaystyle=\bmqty{undef},$ (65)
$\displaystyle\bm{\mathbf{P}}^{x}_{0}$
$\displaystyle=\mathrm{diag}\pqty{\,\bmqty{undef}\,},$ (66)

where
$\sigma^{2}_{p}=$1\text{\times}{10}^{-6}\text{\,}\mathrm{k}\mathrm{m}^{2}$$
and
$\sigma^{2}_{v}=$1\text{\times}{10}^{-6}\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{s}^{-2}$$.
The vector $\bm{\mathbf{q}}(t)=\bmqty{undef}$ is the white noise process with
covariance function

$\mathbb{C}\bqty{\bm{\mathbf{q}}(t),\bm{\mathbf{q}}(s)}=\delta(t-s)\mathrm{diag}\pqty{\,\bmqty{undef}\,}$
(67)

where $\delta(\cdot)$ is the Dirac delta distribution and
$\sigma^{2}_{v}=$2.4\text{\times}{10}^{-5}\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{s}^{-2}$$
and $\sigma^{2}_{\theta}=$0$$. For every generated trajectory we simulated the
range and bearing measurements of the radar using the measurement model

$\bm{\mathbf{z}}_{k}=\bmqty{undef}+\bm{\mathbf{r}}_{k}$ (68)

with the measurement noise covariance

$\bm{\mathbf{R}}=\mathrm{diag}\pqty{\,\bqty{\matrixquantity{$1\text{\times}{10}^{-6}\text{\,}\mathrm{k}\mathrm{m}^{2}$&$0.17\text{\times}{10}^{-6}\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}^{2}$}}\,}.$
(69)

The filters assumed the discrete-time state-space model

$\bm{\mathbf{x}}_{k}=\bm{\mathbf{x}}_{k-1}+\bmqty{undef}$ (70)

obtained by Euler-Maruyama approximation of eq. 61 with step size $\Delta
t=$0.1\text{\,}\mathrm{s}$$. The resulting process covariance was
$\bm{\mathbf{Q}}=\Delta t\cdot\mathrm{diag}\pqty{\,\bmqty{undef}\,}$ where we
set $\sigma^{2}_{\theta}=$1\text{\times}{10}^{-6}$$. We considered a scenario
with misspecified initial conditions, where the filter assumed that the
initial mean and covariance of the system are

$\displaystyle\bm{\mathbf{m}}^{x}_{0|0}$ $\displaystyle=\bmqty{undef},$ (71)
$\displaystyle\bm{\mathbf{P}}^{x}_{0|0}$
$\displaystyle=\mathrm{diag}\pqty{\,\bmqty{undef}\,},$ (72)

which, compared to the true conditions in eqs. 65 and 66, indicates that the
initial filtered covariance is not reflective of our ignorance about the true
initial velocity.

In the experiments, we compared the BSQKF with UT sigma-points with the
standard UKF. The previously proposed GPQKF is not included in the results,
because we were unable to find kernel parameter values yielding numerically
stable estimates.555Experimental evidence strongly suggests that the GPQKF is
generally less suitable for problem where dimension of the state vector is ¿
3. This is were the advantages of the BSQKF, described in Section V-C, become
practically evident. For this experiment, the BSQKF employed the kernel-
agnostic BSQMT, described in Section V-E, with
$\bar{\sigma}^{2}_{f}=$2\text{\times}{10}^{-4}$\bm{\mathbf{I}}$ for the
dynamics and to $\bar{\sigma}^{2}_{h}=\mathrm{diag}\pqty{\,\bmqty{undef}\,}$
for the measurement model. Both filters used identical UT sigma-points with
scaling parameter $\kappa=0$.

Figs. 2 to 4 compare the development of the RMSE and INC in time for the
position, the velocity and the aerodynamic parameter. Overall, the proposed
BSQKF clearly provides more credible estimates for the position and the
aerodynamic parameter, especially during the initial stages, while the
velocity estimates are more pessimistic, which is generally preferable to
optimism in safety-critical applications. In terms of RMSE, the BSQKF gives
comparable estimate of the velocity and superior estimates of position and
parameter. Note, in Fig. 4, how the parameter RMSE for the UKF stays constant,
while for the BSQKF it is much lower. Table IV summarizes the results.

$\displaystyle 0.0$$\displaystyle 0.5$RMSE [$\mathrm{k}\mathrm{m}$]BSQKFUKF$\displaystyle 0$$\displaystyle 25$$\displaystyle 50$$\displaystyle 75$$\displaystyle 100$$\displaystyle 125$$\displaystyle 150$$\displaystyle 175$$\displaystyle 200$time [$\mathrm{s}$]$\displaystyle 0$$\displaystyle 20$$\displaystyle 40$INC [-] Figure 2: Position RMSE and INC. $\displaystyle 0.0$$\displaystyle 0.5$$\displaystyle 1.0$RMSE [$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$]BSQKFUKF$\displaystyle 0$$\displaystyle 50$$\displaystyle 100$$\displaystyle 150$$\displaystyle 200$time [$\mathrm{s}$]$\displaystyle 0$$\displaystyle 25$INC [-] Figure 3: Velocity RMSE and INC. $\displaystyle 0$$\displaystyle 10$RMSE [-]BSQKFUKF$\displaystyle 0$$\displaystyle 50$$\displaystyle 100$$\displaystyle 150$$\displaystyle 200$time [$\mathrm{s}$]$\displaystyle-200$$\displaystyle 0$INC [-] Figure 4: Aerodynamic parameter RMSE and INC. | RMSE |  | INC
---|---|---|---
| BSQKF | UKF |  | BSQKF | UKF
Position [$\mathrm{k}\mathrm{m}$] | 0.018 | 0.025 |  | -0.967 | 1.204
Velocity [$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$] | 0.020 | 0.026 |  | -14.464 | 0.045
Parameter [-] | 0.137 | 15.580 |  | -17.213 | 23.792
Table IV: RMSE and INC averaged over time and simulations for position,
velocity and the aerodynamic parameter.

##  VII Conclusions and Discussion

In this article, we designed a general-purpose moment transformation based on
Bayes–Sard quadrature, which allows for explicit modelling of numerical
integration error through the use of a stochastic process model. The
hierarchical GP prior was shown to be key in developing probabilistic models
whose means correspond to the classical quadrature rules used in the sigma-
point filters and whose variance is statistically meaningful. We designed the
BSQ Kalman filter by leveraging the proposed BSQ moment transform for
computation of the predictive moments.

The BSQMT could easily be applied in nonlinear sigma-point smoothing as well.
Further improvement in performance could be achieved by the use of more
expressive multi-output GP models.

[Computation of Certain Expectations] This appendix derives analytical
expressions for a number of expectations that are needed for implementation of
the Bayes–Sard moment transform. In the following, let
$\bm{\mathbf{x}}\sim\mathrm{N}\pqty{\bm{\mathbf{0}},\,\bm{\mathbf{I}}}$ and
let $\bm{\mathbf{\phi}}({\bm{\mathbf{x}}}):\mathbb{R}^{D}\to\mathbb{R}^{Q}$
consist of monomials defined by the multi-indices
$\bm{\mathbf{\alpha}}_{1},\,\ldots,\,\bm{\mathbf{\alpha}}_{Q}\in\mathbb{N}_{0}^{D}$,
such that $\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})=\bmqty{undef}^{\top}$. Given
a sigma-point set
$\mathcal{X}=\Bqty{{\bm{\mathbf{x}}}_{1},\,\ldots,\,{\bm{\mathbf{x}}}_{N}}\subset\mathbb{R}^{D}$,
let $\bm{\mathbf{k}}({\bm{\mathbf{x}}})$ denote a vector of kernel
evaluations, such that
$\bqty{\bm{\mathbf{k}}({\bm{\mathbf{x}}})missing}_{n}=k({\bm{\mathbf{x}}},{\bm{\mathbf{x}}}^{(n)})$,
where the kernel is the RBF kernel in eq. 44.

###  -A Computing
$\mathbb{E}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}}$

Elements of the $Q\times Q$ matrix
$\bm{\mathbf{A}}\triangleq\mathbb{E}\bqty{\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}}$
are

$[\bm{\mathbf{A}}]_{qr}=\prod_{d=1}^{D}\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\\!{x}_{d}^{\alpha_{q,d}+\alpha_{r,d}}\mathrm{e}^{-{x}_{d}^{2}/2}\,\mathrm{d}{x}_{d}.$

From the standard Gaussian moment formula

$\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\\!{x}^{n}\mathrm{e}^{-{x}^{2}/2}\,\mathrm{d}x=\begin{cases}(n-1)!!&\text{if
$n$ is even,}\\\ 0&\text{if $n$ is odd}\end{cases}$

we obtain

$[\bm{\mathbf{A}}]_{qr}=\prod_{d=1}^{D}(\alpha_{q,d}+\alpha_{r,d}-1)!!$

if none of $\alpha_{q,d}+\alpha_{r,d}$ for $d=1,\,\ldots,\,D$ is odd and
$[\bm{\mathbf{A}}]_{qr}=0$ otherwise.

###  -B Computing
$\mathbb{E}\bqty{\bm{\mathbf{k}}({\bm{\mathbf{x}}})\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}}$

Elements of the $N\times Q$ matrix
$\bm{\mathbf{D}}\triangleq\mathbb{E}\bqty{\bm{\mathbf{k}}({\bm{\mathbf{x}}})\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}}$
are

$[\bm{\mathbf{D}}]_{nq}=\prod_{d=1}^{D}\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}{x}_{d}^{\alpha_{q,d}}\exp\pqty{-\frac{({x}_{d}-{x}_{n,d})^{2}}{2\ell_{d}^{2}}}\mathrm{e}^{-{x}_{d}^{2}/2}\,\mathrm{d}{x}_{d}.$

Completion of squares and a change of variables yield

$\begin{split}\int_{\mathbb{R}}&x_{d}^{\alpha_{q,d}}\exp\pqty{\\!-\frac{(x_{d}-x_{n,d})^{2}}{2\ell_{d}^{2}}}\mathrm{e}^{-x_{d}^{2}/2}\,\mathrm{d}x_{d}\\\
&=\int_{\mathbb{R}}x_{d}^{\alpha_{q,d}}\exp\pqty{\\!-\frac{\pqty{x_{d}-\frac{1}{1+\ell_{d}^{2}}x_{n,d}}^{2}}{2\frac{\ell_{d}^{2}}{1+\ell_{d}^{2}}}}\\\
&\hskip
56.9055pt\times\exp\pqty{\\!-\frac{x_{n,d}^{2}}{2(1+\ell_{d}^{2})}}\,\mathrm{d}x_{d}\\\
&=\pqty{\frac{\ell_{d}^{2}}{1+\ell_{d}^{2}}}^{1/2}\exp\pqty{\\!-\frac{x_{n,d}^{2}}{2(1+\ell_{d}^{2})}}\\\
&\hskip
14.22636pt\times\int_{\mathbb{R}}\bqty{\pqty{\frac{\ell_{d}^{2}}{1+\ell_{d}^{2}}}^{1/2}x_{d}+\frac{x_{n,d}}{1+\ell_{d}^{2}}}^{\alpha_{q,d}}\mathrm{e}^{-x_{d}^{2}/2}\,\mathrm{d}x_{d}\\\
&=\ell_{d}\pqty{\frac{\ell_{d}^{2}}{1+\ell_{d}^{2}}}^{(1+\alpha_{q,d})/2}\exp\pqty{\\!-\frac{x_{n,d}^{2}}{2(1+\ell_{d}^{2})}}\\\
&\hskip
14.22636pt\times\int_{\mathbb{R}}\pqty{\ell_{d}x_{d}+\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}}\mathrm{e}^{-x_{d}^{2}/2}\,\mathrm{d}x_{d}.\end{split}$

The binomial theorem yields

$\begin{split}&\pqty{\ell_{d}x_{d}+\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}}\\\
&\phantom{\bigg{(}\ell_{d}x_{d}}=\sum_{s=0}^{\alpha_{q,d}}\frac{\alpha_{q,d}!}{s!(\alpha_{q,d}-s)!}\ell_{d}^{s}x_{d}^{s}\pqty{\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}-s}.\end{split}$

Therefore

$\begin{split}&\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\pqty{\ell_{d}x_{d}+\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2})}}}^{\alpha_{q,d}}\mathrm{e}^{-x_{d}^{2}/2}\,\mathrm{d}x_{d}\\\
&=\sum_{s=0}^{\alpha_{q,d}}\frac{\alpha_{q,d}!}{s!(\alpha_{q,d}-s)!}\ell_{d}^{s}\pqty{\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}-s}\\\
&\hskip
28.45274pt\times\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}x_{d}^{s}\mathrm{e}^{-x_{d}^{2}/2}\,\mathrm{d}x_{d}\\\
&=\sum_{s=0}^{\lfloor\alpha_{q,d}/2\rfloor}\frac{\alpha_{q,d}!}{(2s)!(\alpha_{q,d}-2s)!}\ell_{d}^{2s}\pqty{\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}-2s}\frac{(2s)!}{2^{s}s!}\\\
&=\sum_{s=0}^{\lfloor\alpha_{q,d}/2\rfloor}\frac{\alpha_{q,d}!}{2^{s}s!(\alpha_{q,d}-2s)!}\ell_{d}^{2s}\pqty{\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}}^{\alpha_{q,d}-2s}.\end{split}$

By combining all the equations above, we conclude that

$\begin{split}[\bm{\mathbf{D}}]_{nq}={}&\prod_{d=1}^{D}\Bigg{[}\ell_{d}\bigg{(}\frac{\ell_{d}^{2}}{1+\ell_{d}^{2}}\bigg{)}^{(1+\alpha_{q,d})/2}\exp\bigg(\\!-\frac{x_{n,d}^{2}}{2(1+\ell_{d}^{2})}\bigg{missing})\Bigg{.}\\\
&\times\Bigg{.}\sum_{s=0}^{\lfloor\alpha_{q,d}/2\rfloor}\frac{\alpha_{q,d}!}{2^{s}s!(\alpha_{q,d}-2s)!}\ell_{d}^{2s}\bigg{(}\frac{x_{n,d}}{\sqrt{1+\ell_{d}^{2}}}\bigg{)}^{\alpha_{q,d}-2s}\Bigg{]}.\end{split}$

###  -C Computing
$\mathbb{E}\bqty{{\bm{\mathbf{x}}}\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})missing}$

This reduces to the computation in Section -A. The $D\times Q$ matrix
$\bm{\mathbf{B}}\triangleq\mathbb{E}\bqty{{\bm{\mathbf{x}}}\bm{\mathbf{\phi}}({\bm{\mathbf{x}}})^{\top}}$
has the elements

$\bqty{\bm{\mathbf{B}}}_{dq}=\alpha_{q,d}\prod_{d\neq q}(\alpha_{q,d}-1)!!$

if $\alpha_{q,d}$ is odd and none of $\alpha_{q,d}$ are and
$\bqty{\bm{\mathbf{B}}}_{dq}=0$ otherwise.

## Acknowledgment

JP and OS were supported by the Czech Ministry of Education, Youth and Sports,
project LO1506. TK was supported by the Aalto ELEC Doctoral School. CJO was
supported by the Lloyd’s Register Foundation programme on Data-Centric
Engineering at the Alan Turing Institute, UK. SS was supported by the Academy
of Finland project 313708.

## References

  * [1] M. S. Grewal, L. R. Weill, and A. P. Andrews, _Global Positioning Systems, Inertial Navigation, and Integration_. Wiley, 2007.
  * [2] A. J. Haug, _Bayesian Estimation and Tracking: A Practical Guide_. Wiley-Blackwell, 2012.
  * [3] L. D. Stone, R. L. Streit, T. L. Corwin, and K. L. Bell, _Bayesian Multiple Target Tracking_. Artech House, 2014.
  * [4] H. Durrant-Whyte and T. Bailey, “Simultaneous localization and mapping: part I,” _IEEE Robotics & Automation Magazine_, vol. 13, no. 2, pp. 99–110, jun 2006.
  * [5] S. Gillijns, O. Mendoza, J. Chandrasekar, B. De Moor, D. Bernstein, and A. Ridley, “What is the ensemble Kalman filter and how well does it work?” in _American Control Conference, 2006_ , June 2006, p. 6.
  * [6] R. Bhar, _Stochastic filtering with applications in finance_. World Scientific, 2010.
  * [7] F. Gustafsson, “Particle Filter Theory and Practice with Positioning Applications,” _IEEE Aerospace and Electronic Systems Magazine_ , vol. 25, no. 7, pp. 53–81, 2010.
  * [8] M. Šimandl, J. Královec, and T. Söderström, “Advanced point-mass method for nonlinear state estimation,” _Automatica_ , vol. 42, no. 7, pp. 1133–1145, 2006.
  * [9] S. Holmes, G. Klein, and D. W. Murray, “A Square Root Unscented Kalman Filter for visual monoSLAM,” in _2008 IEEE International Conference on Robotics and Automation_ , May 2008, pp. 3710–3716.
  * [10] K. Lowrey, J. Dao, and E. Todorov, “Real-time State Estimation with Whole-body Multi-contact Dynamics: A Modified UKF Approach,” in _IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids)_ , Nov 2016, pp. 1225–1232.
  * [11] F. Larkin, “Gaussian measure in hilbert space and applications in numerical analysis,” _The Rocky Mountain Journal of Mathematics_ , pp. 379–421, 1972\.
  * [12] A. O’Hagan, “Bayes–Hermite quadrature,” _Journal of Statistical Planning and Inference_ , vol. 29, no. 3, pp. 245–260, 1991.
  * [13] T. Karvonen and S. Särkkä, “Classical quadrature rules via Gaussian processes,” in _IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_ , 2017.
  * [14] C. J. Oates, S. Niederer, A. Lee, F.-X. Briol, and M. Girolami, “Probabilistic Models for Integration Error in the Assessment of Functional Cardiac Models,” in _Advances in Neural Informations Processing Systems (NIPS 2017)_ , 2017, pp. 110–118.
  * [15] F.-X. Briol, C. J. Oates, M. Girolami, M. A. Osborne, and D. Sejdinovic, “Probabilistic Integration: A Role in Statistical Computation? (with discussion and rejoinder),” _Statistical Science_ , 2018, to appear.
  * [16] T. Karvonen, S. Särkkä, C. Oates _et al._ , “Symmetry Exploits for Bayesian Cubature Methods,” _ArXiv e-prints_ , 2018, arXiv:1809.10227.
  * [17] S. Särkkä, J. Hartikainen, L. Svensson, and F. Sandblom, “On the relation between Gaussian process quadratures and sigma-point methods,” _Journal of Advances in Information Fusion_ , vol. 11, pp. 31–46, June 2016\.
  * [18] J. Prüher and O. Straka, “Gaussian Process Quadrature Moment Transform,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 9, pp. 2844–2854, Sep. 2018.
  * [19] P. Diaconis, “Bayesian numerical analysis,” in _Statistical Decision Theory and Related Topics IV_ , S. S. Gupta and J. O. Berger, Eds. Springer, 1988, pp. 163–175.
  * [20] T. Karvonen, C. J. Oates, and S. Särkkä, “A Bayes-Sard Cubature Method,” _Advances in Neural Informations Processing Systems (NIPS 2018)_ , pp. 5882–5893, 2018.
  * [21] K. Xiong, H. Zhang, and C. W. Chan, “Performance evaluation of UKF-based nonlinear filtering,” _Automatica_ , vol. 42, no. 2, pp. 261–270, 2006.
  * [22] K. Xiong, L. D. Liu, and H. Y. Zhang, “Modified unscented Kalman filtering and its application in autonomous satellite navigation,” _Aerospace Science and Technology_ , vol. 13, no. 4, pp. 238–246, 2009.
  * [23] T. Karvonen, S. Bonnabel, E. Moulines, and S. Särkkä, “On Stability of a Class of Filters for Non-Linear Stochastic Systems,” _ArXiv e-prints_ , Sep. 2018, arXiv:1809.05667v1.
  * [24] A. A. Clifford, _Multivariate error analysis: a handbook of error propagation and calculation in many-parameter systems_. John Wiley & Sons, 1973.
  * [25] I. Arasaratnam and S. Haykin, “Cubature Kalman Filters,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 6, pp. 1254–1269, 2009.
  * [26] I. Arasaratnam, S. Haykin, and R. Elliott, “Discrete-Time Nonlinear Filtering Algorithms Using Gauss-Hermite Quadrature,” _Proceedings of the IEEE_ , vol. 95, no. 5, pp. 953–977, 2007.
  * [27] Y. Wu, D. Hu, M. Wu, and X. Hu, “A Numerical-Integration Perspective on Gaussian Filters,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 8, pp. 2910–2921, 2006.
  * [28] K. Ito and K. Xiong, “Gaussian Filters for Nonlinear Filtering Problems,” _IEEE Transactions on Automatic Control_ , vol. 45, no. 5, pp. 910–927, 2000\.
  * [29] G. H. Golub and J. H. Welsch, “Calculation of Gauss quadrature rules,” _Mathematics of Computation_ , vol. 23, no. 106, pp. 221–221, may 1969.
  * [30] B. Jia, M. Xin, and Y. Cheng, “Sparse-grid quadrature nonlinear filtering,” _Automatica_ , vol. 48, no. 2, pp. 327–341, 2012.
  * [31] F. M. Larkin, “Optimal approximation in Hilbert spaces with reproducing kernel functions,” _Mathematics of Computation_ , vol. 24, no. 112, pp. 911–921, 1970.
  * [32] A. Sommariva and M. Vianello, “Numerical cubature on scattered data by radial basis functions,” _Computing_ , vol. 76, no. 3–4, pp. 295–310, 2006.
  * [33] J. Oettershagen, “Construction of optimal cubature algorithms with applications to econometrics and uncertainty quantification,” Ph.D. dissertation, Institut für Numerische Simulation, Universität Bonn, 2017\.
  * [34] C. E. Rasmussen and Z. Ghahramani, “Bayesian Monte Carlo,” in _Advances in Neural Information Processing Systems 15_ , S. Becker, S. Thrun, and K. Obermayer, Eds. MIT Press, 2003, pp. 505–512.
  * [35] C. E. Rasmussen and C. K. Williams, _Gaussian Processes for Machine Learning_. The MIT Press, 2006.
  * [36] M. P. Deisenroth, M. F. Huber, and U. D. Hanebeck, “Analytic moment-based Gaussian process filtering,” in _Proceedings of the 26th Annual International Conference on Machine Learning - ICML ’09_. ACM Press, 2009, pp. 1–8.
  * [37] M. Álvarez, L. Rosasco, and N. Lawrence, “Kernels for Vector-valued Functions: A Review,” _Foundations and Trends in Machine Learning_ , vol. 4, no. 3, pp. 195–266, 2012.
  * [38] X. Xi, F.-X. Briol, and M. Girolami, “Bayesian Quadrature for Multiple Related Integrals,” in _Proceedings of the 35th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, vol. 80. PMLR, 10–15 Jul 2018, pp. 5373–5382.
  * [39] A. Shah, A. G. Wilson, and Z. Ghahramani, “Student-t Processes as Alternatives to Gaussian Processes,” in _Proceedings of the 17th International Conference on Artificial Intelligence and Statistics (AISTATS)_ , vol. 33, 2014, pp. 877–885.
  * [40] J. Prüher, F. Tronarp, T. Karvonen, S. Särkkä, and O. Straka, “Student-$t$ Process Quadratures for Filtering of Non-linear Systems with Heavy-tailed Noise,” in _20th International Conference on Information Fusion (Fusion)_ , July 2017, pp. 1–8.
  * [41] T. J. Santner, B. J. Williams, and W. I. Notz, _The Design and Analysis of Computer Experiments_ , ser. Springer Series in Statistics. Springer, 2003.
  * [42] F. M. Larkin, “Probabilistic error estimates in spline interpolation and quadrature,” in _Information Processing 74 (Proceedings of IFIP Congress, Stockholm, 1974)_ , vol. 74. North-Holland, 1974, pp. 605–609.
  * [43] M. Kanagawa, B. K. Sriperumbudur, and K. Fukumizu, “Convergence guarantees for kernel-based quadrature rules in misspecified settings,” in _Advances in Neural Information Processing Systems 29_. Curran Associates, Inc., 2016, pp. 3288–3296.
  * [44] ——, “Convergence analysis of deterministic kernel-based quadrature rules in misspecified settings,” _ArXiv e-prints_ , 2017, arXiv:1709.00147.
  * [45] S. Särkkä, _Bayesian Filtering and Smoothing_. New York: Cambridge University Press, 2013.
  * [46] S. Haykin, _Kalman Filtering and Neural Networks_. Wiley, 2001.
  * [47] N. J. Gordon, D. J. Salmond, and A. F. M. Smith, “Novel approach to nonlinear/non-Gaussian Bayesian state estimation,” _IEE Proceedings F (Radar and Signal Processing)_ , vol. 140, no. 2, pp. 107–113, 1993.
  * [48] G. Kitagawa, “Monte Carlo Filter and Smoother for Non-Gaussian Nonlinear State Space Models,” _Journal of Computational and Graphical Statistics_ , vol. 5, no. 1, pp. 1–25, 1996.
  * [49] X. R. Li and Z. Zhao, “Measuring Estimator’s Credibility: Noncredibility Index,” in _Information Fusion, 2006 9th International Conference on_ , 2006, pp. 1–8.
  * [50] S. J. Julier and J. K. Uhlmann, “Unscented Filtering and Nonlinear Estimation,” _Proceedings of the IEEE_ , vol. 92, no. 3, pp. 401–422, 2004\.

|  Jakub Prüher is a Ph.D. student in the Identification and Decision Making
research group at NTIS - New Technologies for the Information Society. He
received his bachelor’s degree in computer science in 2011 and his master’s
degree in cybernetics in 2013, both from the University of West Bohemia. His
research interests include nonlinear state estimation methods, use of Bayesian
reasoning in numerical algorithms and machine learning in general.
---|---
|  Toni Karvonen received his Master of Science degree in applied mathematics
from the University of Helsinki in 2015. Since 2016, he has been a doctoral
student at the Department of Electrical Engineering and Automation, Aalto
University, Finland. His research interests are in stochastic state estimation
and numerical analysis, in particular numerical integration and probabilistic
numerics.
---|---
|  Chris J. Oates received his Ph.D. degree in statistics and complexity
science from the University of Warwick, Coventry, UK, in 2013. Since 2018, he
has been a Professor in the School of Mathematics, Statistics and Physics,
Newcastle University, Newcastle-upon-Tyne, UK. He is also a Group Leader for
the Lloyd’s Register Foundation programme on Data-Centric Engineering at the
Alan Turing Institute, London, UK, where he focuses on computational
statistical techniques and the development of probabilistic numerical methods.
Prof. Oates was a recipient of the biennial Research Prize of the Royal
Statistical Society in 2017.
---|---
|  Ondřej Straka (M’16) received the master’s degree in cybernetics and
control engineering and the Ph.D. degree in cybernetics from the University of
West Bohemia, Pilsen, Czech Republic, in 1998 and 2004, respectively. Since
2015, he has been an Associate Professor with the Department of Cybernetics,
University of West Bohemia. He is the Head of the Identification and Decision
Making Research Group (IDM), NTIS—New Technologies for the Information
Society, where he focuses on nonlinear state estimation, stochastic systems,
and system identification. He has participated in a number of projects of
fundamental research and in several project of applied research (e.g., GNSS-
based safe train localization and attitude and heading reference system). He
was involved in development of several software frameworks for nonlinear state
estimation and system identification. He has published over 70 journal and
conference papers in journals, such as Automatica, the IEEE Transactions on
Automatic Control, the IEEE Transactions on Aerospace and Electronic Systems,
the IEEE Transactions on Cybernetics, and Signal Processing and at
international conferences such as American Control Conference, World
Congresses and Symposia of the IFAC, and FUSION Conferences. His current
research interests include local and global nonlinear state estimation
methods, noise CM estimation in state-space models, performance evaluation,
and fault detection in navigation systems. Dr. Straka was a recipient of
Werner von Siemens Excellence Award in 2014 for most important result in the
basic research.
---|---
|  Simo Särkkä received his Doctor of Science (Tech.) degree in electrical and
communications engineering from Helsinki University of Technology, Espoo,
Finland, in 2006. From 2000 to 2010 he worked with Nokia Ltd., Indagon Ltd.,
and Nalco Company in various industrial research projects related to
telecommunications, positioning systems, and industrial process control. From
2010 to 2015 he worked as a Senior Researcher with the Department of
Biomedical Engineering and Computational Science (BECS) of Aalto University,
Finland. Currently, he is an Associate Professor with the Department of
Electrical Engineering and Automation of Aalto University. His research
interests are in multi-sensor data processing systems with applications in
location sensing, health and medical technology, machine learning, inverse
problems, and brain imaging. He is a Senior Member of IEEE and serving as an
Associate Editor of IEEE Signal Processing Letters.
---|---