Robust Filtering and Smoothing with Gaussian Processes

I Introduction

Filtering and smoothing in context of dynamic systems refers to a Bayesian methodology for computing posterior distributions of the latent state based on a history of noisy measurements. This kind of methodology can be found, e.g., in navigation, control engineering, robotics, and machine learning . Solutions to filtering and smoothing in linear dynamic systems are well known, and numerous approximations for nonlinear systems have been proposed, for both filtering and smoothing .

In this article, we focus on Gaussian filtering and smoothing in Gaussian process (GP) dynamic systems. GPs are a robust non-parametric method for approximating unknown functions by a posterior distribution over them . Although GPs have been around for decades, they only recently became computationally interesting for applications in robotics, control, and machine learning .

The contribution of this article is the derivation of a novel, principled, and robust Rauch-Tung-Striebel (RTS) smoother for GP dynamic systems, which we call the GP-RTSS. The GP-RTSS computes a Gaussian approximation to the smoothing distribution in closed form. The posterior filtering and smoothing distributions can be computed without linearization or (small) sampling approximations of densities .

We provide numerical evidence that the GP-RTSS is more robust than state-of-the-art nonlinear Gaussian filtering and smoothing algorithms including the extended Kalman filter (EKF) , the unscented Kalman filter (UKF) , the cubature Kalman filter (CKF) , the GP-UKF , and their corresponding RTS smoothers. Robustness refers to the ability of an inferred distribution to explain the “true” state/measurement.

The paper is structured as follows: In Secs. I-A–I-B, we introduce the problem setup and necessary background on Gaussian smoothing and GP dynamic systems. In Sec. II, we briefly introduce Gaussian process regression, discuss the expressiveness of a GP, and explain how to train GPs. Sec. III details our proposed method (GP-RTSS) for smoothing in GP dynamic systems. In Sec. IV, we provide experimental evidence of the robustness of the GP-RTSS. Sec. V concludes the paper with a discussion.

In this article, we consider discrete-time stochastic systems

where ${\boldsymbol{\mathbf{x}}}_{t}\in\mathds{R}^{D}$ is the state, ${\boldsymbol{\mathbf{z}}}_{t}\in\mathds{R}^{E}$ is the measurement at time step $t$, ${\boldsymbol{\mathbf{w}}}_{t}\sim\mathcal{N}({\boldsymbol{\mathbf{0}}},{\mathbf{\Sigma}}_{w})$ is Gaussian system noise, ${\boldsymbol{\mathbf{v}}}_{t}\sim\mathcal{N}({\boldsymbol{\mathbf{0}}},{\mathbf{\Sigma}}_{v})$ is Gaussian measurement noise, $f$ is the transition function (or system function) and $g$ is the measurement function. The discrete time steps $t$ run from 0 to $T$. The initial state ${\boldsymbol{\mathbf{x}}}_{0}$ of the time series is distributed according to a Gaussian prior distribution ${p}({\boldsymbol{\mathbf{x}}}_{0})=\mathcal{N}({\boldsymbol{\mathbf{\mu}}}_{0}^{x},{\mathbf{\Sigma}}_{0}^{x})$. The purpose of filtering and smoothing is to find approximations to the posterior distributions ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:\tau})$, where $1\!\!:\!\!\tau$ in a subindex abbreviates $1,\dotsc,\tau$ with $\tau=t$ during filtering and $\tau=T$ during smoothing. In this article, we consider Gaussian approximations ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:\tau})\approx\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{\mu}}}_{t|\tau}^{x},{\mathbf{\Sigma}}_{t|\tau}^{x})$ of the latent state posterior distributions ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:\tau})$. We use the short-hand notation ${\boldsymbol{\mathbf{a}}}_{b|c}^{d}$ where ${\boldsymbol{\mathbf{a}}}={\boldsymbol{\mathbf{\mu}}}$ denotes the mean ${\boldsymbol{\mathbf{\mu}}}$ and ${\boldsymbol{\mathbf{a}}}={\mathbf{\Sigma}}$ denotes the covariance, $b$ denotes the time step under consideration, $c$ denotes the time step up to which we consider measurements, and $d\in\{x,z\}$ denotes either the latent space ($x$) or the observed space ($z$).

I-B Gaussian RTS Smoothing

Given the filtering distributions ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t})=\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{\mu}}}_{t|t}^{x},{\mathbf{\Sigma}}_{t|t}^{x})$, $t=1,\dotsc,T$, a sufficient condition for Gaussian smoothing is the computation of Gaussian approximations of the joint distributions ${p}({\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$, $t=1,\dotsc,T$ .

In Gaussian smoothers, the standard smoothing distribution for the dynamic system in Eqs. (1)–(2) is always

Depending on the methodology of computing this joint distribution, we can directly derive arbitrary RTS smoothing algorithms, including the URTSS , the EKS , the CKS , a smoothing extension to the CKF , or the GP-URTSS, a smoothing extension to the GP-UKF . The individual smoothers (URTSS, EKS, CKS, GP-based smoothers etc.) simply differ in the way of computing/estimating the means and covariances required in Eqs. (4)–(6) .

To derive the GP-URTSS, we closely follow the derivation of the URTSS . The GP-URTSS is a novel smoother, but its derivation is relatively straightforward and therefore not detailed in this article. Instead, we detail the derivation of the GP-RTSS, a robust Rauch-Tung-Striebel smoother for GP dynamic systems, which is based on analytic computation of the means and (cross-)covariances in Eqs. (4)–(6).

In GP dynamics systems, the transition function $f$ and the measurement function $g$ in Eqs. (1)–(2) are modeled by Gaussian processes. This setup is getting more relevant in practical applications such as robotics and control, where it can be difficult to find an accurate parametric form of $f$ and $g$, respectively . Given the increasing use of GP models in robotics and control, the robustness of Bayesian state estimation is important.

II Gaussian Processes

In the standard GP regression model, we assume that the data $\mathcal{D}\coloneqq\{{\mathbf{X}}\coloneqq[{\boldsymbol{\mathbf{x}}}_{1},\dotsc,{\boldsymbol{\mathbf{x}}}_{n}]^{\top},\,{\boldsymbol{\mathbf{y}}}\coloneqq[y_{1},\dots,y_{n}]^{\top}\}$ have been generated according to $y_{i}=h({\boldsymbol{\mathbf{x}}}_{i})+\varepsilon_{i}$, where $h:\mathds{R}^{D}\to\mathds{R}$ and $\varepsilon_{i}\sim\mathcal{N}(0,\sigma_{\varepsilon}^{2})$ is independent (measurement) noise. GPs consider $h$ a random function and infer a posterior distribution over $h$ from data. The posterior is used to make predictions about function values $h({\boldsymbol{\mathbf{x}}}_{*})$ for arbitrary inputs ${\boldsymbol{\mathbf{x}}}_{*}\in\mathds{R}^{D}$.

Similar to a Gaussian distribution, which is fully specified by a mean vector and a covariance matrix, a GP is fully specified by a mean function $m_{h}(\,\cdot\,)$ and a covariance function

which specifies the covariance between any two function values. Here, $\mathds{E}_{h}$ denotes the expectation with respect to the function $h$. The covariance function $k_{h}(\,\cdot\,,\,\cdot\,)$ is also called a kernel.

Unless stated otherwise, we consider a prior mean function $m_{h}\equiv 0$ and use the squared exponential (SE) covariance function with automatic relevance determination

Although the SE covariance function and the zero-prior mean function are common defaults, they retain a great deal of expressiveness. Inspired by , we demonstrate this expressiveness and show the correspondence of our GP model to a universal function approximator: Consider a function

where $\gamma_{n}\sim\mathcal{N}(0,1)\,,n=1,\dotsc,N$. Note that in the limit $h(x)$ is represented by infinitely many Gaussian-shaped basis functions along the real axis with variance $\lambda^{2}$ and prior (Gaussian) random weights $\gamma_{n}$, for $x\in\mathds{R}$, and for all $i\in\mathds{Z}$. The model in Eq. (10) is considered a universal function approximator. Writing the sums in Eq. (10) as an integral over the real axis $\mathds{R}$, we obtain

where $\gamma(s)\sim\mathcal{N}(0,1)$ is a white-noise process and $\mathcal{K}$ is a Gaussian convolution kernel. The function values of $h$ are jointly normal, which follows from the convolution $\gamma*\mathcal{K}$. We now analyze the mean function and the covariance function of $h$, which fully specify the distribution of $h$. The only random variables are the weights $\gamma(s)$. Computing the expected function of this model (prior mean function) requires averaging over $\gamma(s)$ and yields

since $\mathds{E}_{\gamma}[\gamma(s)]=0$. Hence, the mean function of $h$ equals zero everywhere. Let us now find the covariance function. Since the mean function equals zero, for any $x,x^{\prime}\in\mathds{R}$ we obtain

From Eqs. (13) and (15), we see that the mean function and the covariance function of the universal function approximator in Eq. (10) correspond to the GP model assumptions we made earlier: a prior mean function $m_{h}\equiv 0$ and the SE covariance function in Eq. (9) for a one-dimensional input space. Hence, the considered GP prior implicitly assumes latent functions $h$ that can be described by the universal function approximator in Eq. (11). Examples of covariance functions that encode different model assumptions are given in .

II-B Training via Evidence Maximization

For $E$ target dimensions, we train $E$ GPs assuming that the target dimensions are independent at a deterministically given test input (if the test input is uncertain, the target dimensions covary): After observing a data set $\mathcal{D}$, for each (training) target dimension, we learn the $D+1$ hyper-parameters of the covariance function and the noise variance of the data using evidence maximization : Collecting all $(D+2)E$ hyper-parameters in the vector ${\boldsymbol{\mathbf{\theta}}}$, evidence maximization yields a point estimate $\hat{{\boldsymbol{\mathbf{\theta}}}}\in\operatornamewithlimits{argmax}_{{\boldsymbol{\mathbf{\theta}}}}\log{p}({\boldsymbol{\mathbf{y}}}|{\mathbf{X}},{\boldsymbol{\mathbf{\theta}}})$. Evidence maximization automatically trades off data fit with function complexity and avoids overfitting .

From here onward, we consider the GP dynamics system setup, where two GP models have been trained using evidence maximization: $\mathcal{GP}_{f}$, which models the mapping ${\boldsymbol{\mathbf{x}}}_{t-1}\mapsto{\boldsymbol{\mathbf{x}}}_{t},~{}\mathds{R}^{D}\to\mathds{R}^{D}$, see Eq. (1), and $\mathcal{GP}_{g}$, which models the mapping ${\boldsymbol{\mathbf{x}}}_{t}\mapsto{\boldsymbol{\mathbf{z}}}_{t},~{}\mathds{R}^{D}\to\mathds{R}^{E}$, see Eq. (2). To keep the notation uncluttered, we do not explicitly condition on the hyper-parameters $\hat{{\boldsymbol{\mathbf{\theta}}}}$ and the training data $\mathcal{D}$ in the following.

III Robust Smoothing in Gaussian Process Dynamic Systems

Analytic moment-based filtering in GP dynamic systems has been proposed in , where the filter distribution is given by

for $t=1,\dotsc,T$. Here, we extend these filtering results to analytic moment-based smoothing, where we explicitly take nonlinearities into account (no linearization required) while propagating full Gaussian densities (no sigma/cubature-point representation required) through nonlinear GP models.

In the following, we detail our novel RTS smoothing approach for GP dynamic systems. We fit our smoother in the standard frame of Eqs. (4)–(6). For this, we compute the means and covariances of the Gaussian approximation

to the joint ${p}({\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$, after which the smoother is fully determined . Our approximation does not involve sampling, linearization, or numerical integration. Instead, we present closed-form expressions of a deterministic Gaussian approximation of the joint distribution in Eq. (19).

Using the system Eq. (1) and integrating over all three sources of uncertainties (the system noise, the state ${\boldsymbol{\mathbf{x}}}_{t-1}$, and the system function itself), we apply the law of total expectation and obtain the marginal mean

The expectations in Eq. (20) are taken with respect to the posterior GP distribution ${p}(f)$ and the filter distribution ${p}({\boldsymbol{\mathbf{x}}}_{t-1}|{\boldsymbol{\mathbf{z}}}_{1:t-1})=\mathcal{N}({\boldsymbol{\mathbf{\mu}}}_{t-1|t-1}^{x},{\mathbf{\Sigma}}_{t-1|t-1}^{x})$ at time step $t-1$. Eq. (20) can be rewritten as ${\boldsymbol{\mathbf{\mu}}}_{t|t-1}^{x}=\mathds{E}_{{\boldsymbol{\mathbf{x}}}_{t-1}}[m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})|{\boldsymbol{\mathbf{z}}}_{1:t-1}]$ with $m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})\coloneqq\mathds{E}_{f}[f({\boldsymbol{\mathbf{x}}}_{t-1})|{\boldsymbol{\mathbf{x}}}_{t-1}]$ is the posterior mean function of $\mathcal{GP}_{f}$. Writing $m_{f}$ as a finite sum over the SE kernels centered at all $n$ training inputs , the predicted mean for each target dimension $a=1,\dotsc,D$ is

where ${p}({\boldsymbol{\mathbf{x}}}_{t-1}|{\boldsymbol{\mathbf{z}}}_{1:t-1})=\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t-1}\,|\,{\boldsymbol{\mathbf{\mu}}}_{t-1|t-1}^{x},{\mathbf{\Sigma}}_{t-1|t-1}^{x})$ is the filter distribution at time $t-1$. Moreover, ${\boldsymbol{\mathbf{x}}}_{i}$, $i=1,\dotsc,n$, are the training set of $\mathcal{GP}_{f}$, $k_{f_{a}}$ is the covariance function of $\mathcal{GP}_{f}$ for the $a$th target dimension (GP hyper-parameters are not shared across dimensions), and ${\boldsymbol{\mathbf{\beta}}}_{a}^{x}\coloneqq({\mathbf{K}}_{f_{a}}+\sigma_{w_{a}}^{2}{\mathbf{I}})^{-1}{\boldsymbol{\mathbf{y}}}_{a}\in\mathds{R}^{n}$. For dimension $a$, ${\mathbf{K}}_{f_{a}}$ denotes the kernel matrix (Gram matrix), where ${\mathbf{K}}_{f_{a_{ij}}}=k_{f_{a}}({\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j})$, $i,j=1,\dotsc,n$. Moreover, ${\boldsymbol{\mathbf{y}}}_{a}$ are the training targets, and $\sigma_{w_{a}}^{2}$ is the learned system noise variance. The vector ${\boldsymbol{\mathbf{\beta}}}_{a}^{x}$ has been pulled out of the integration since it is independent of ${\boldsymbol{\mathbf{x}}}_{t-1}$. Note that ${\boldsymbol{\mathbf{x}}}_{t-1}$ serves as a test input from the perspective of the GP regression model.

For the SE covariance function in Eq. (9), the integral in (21) can be computed analytically (other tractable choices are covariance functions containing combinations of squared exponentials, trigonometric functions, and polynomials). The marginal mean is given as

$i=1,\dotsc,n$, being the solution to the integral in Eq. (21). Here, $\alpha_{f_{a}}^{2}$ is the signal variance of the $a$th target dimension of $\mathcal{GP}_{f}$, a learned hyper-parameter of the SE covariance function, see Eq. (9).

III-A2 Marginal Covariance Matrix

We now explicitly compute the entries of the corresponding covariance ${\mathbf{\Sigma}}_{t|t-1}^{x}$. Using the law of total covariance, we obtain for $a,b=1,\dotsc,D$

where we exploited in the last term that the system noise ${\boldsymbol{\mathbf{w}}}$ has mean zero. Note that Eq. (25) is the sum of the covariance of (conditional) expected values and the expectation of a (conditional) covariance. We analyze these terms in the following.

The covariance of the expectations in Eq. (25) is

where we used that $\mathds{E}_{f}[f({\boldsymbol{\mathbf{x}}}_{t-1})|{\boldsymbol{\mathbf{x}}}_{t-1}]=m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})$. With ${\boldsymbol{\mathbf{\beta}}}_{a}^{x}=({\mathbf{K}}_{a}+\sigma_{w_{a}}^{2}{\mathbf{I}})^{-1}{\boldsymbol{\mathbf{y}}}_{a}$ and $m_{f_{a}}({\boldsymbol{\mathbf{x}}}_{t-1})=k_{f_{a}}({\mathbf{X}},{\boldsymbol{\mathbf{x}}}_{t-1})^{\top}{\boldsymbol{\mathbf{\beta}}}^{x}_{a}$, we obtain

Following , the entries of ${\mathbf{Q}}\in\mathds{R}^{n\times n}$ are given as

where we defined ${\mathbf{R}}\coloneqq{\mathbf{\Sigma}}_{t-1|t-1}^{x}({\mathbf{\Lambda}}_{a}^{-1}+{\mathbf{\Lambda}}_{b}^{-1})+{\mathbf{I}}$, ${\boldsymbol{\mathbf{\zeta}}}_{i}\coloneqq{\boldsymbol{\mathbf{x}}}_{i}-{\boldsymbol{\mathbf{\mu}}}_{t-1|t-1}^{x}$, and ${\boldsymbol{\mathbf{z}}}_{ij}\coloneqq{\mathbf{\Lambda}}_{a}^{-1}{\boldsymbol{\mathbf{\zeta}}}_{i}+{\mathbf{\Lambda}}_{b}^{-1}{\boldsymbol{\mathbf{\zeta}}}_{j}$.

The expected covariance in Eq. (25) is given as

since the noise covariance matrix ${\mathbf{\Sigma}}_{w}$ is diagonal. Following our GP training assumption that different target dimensions do not covary if the input is deterministically given, Eq. (29) is only non-zero if $a=b$, i.e., Eq. (29) plays a role only for diagonal entries of ${\mathbf{\Sigma}}_{t|t-1}^{x}$. For these diagonal entries ($a=b$), the expected covariance in Eq. (29) is

Hence, the desired marginal covariance matrix in Eq. (25) is

We have now solved for the marginal distribution ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$ in Eq. (19). Since the approximate Gaussian filter distribution ${p}({\boldsymbol{\mathbf{x}}}_{t-1}|{\boldsymbol{\mathbf{z}}}_{1:t-1})=\mathcal{N}({\boldsymbol{\mathbf{\mu}}}_{t-1|t-1}^{x},{\mathbf{\Sigma}}_{t-1|t-1}^{x})$ is also known, it remains to compute the cross-covariance ${\mathbf{\Sigma}}_{t-1,t|t-1}^{x}$ to fully determine the Gaussian approximation in Eq. (19).

III-B Cross-Covariance

By the definition of a covariance and the system Eq. (1), the missing cross-covariance matrix ${\mathbf{\Sigma}}_{t-1,t|t-1}^{x}$ in Eq. (19) is

where ${\boldsymbol{\mathbf{\mu}}}_{t-1|t-1}^{x}$ is the mean of the filter update at time step $t-1$ and ${\boldsymbol{\mathbf{\mu}}}_{t|t-1}^{x}$ is the mean of the time update, see Eq. (20). Note that we explicitly average out the model uncertainty about $f$. Using the law of total expectations, we obtain

where we used the fact that $\mathds{E}_{f,{\boldsymbol{\mathbf{w}}}_{t}}[f({\boldsymbol{\mathbf{x}}}_{t-1})+{\boldsymbol{\mathbf{w}}}_{t}|{\boldsymbol{\mathbf{x}}}_{t-1}]=m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})$ is the mean function of $\mathcal{GP}_{f}$, which models the mapping ${\boldsymbol{\mathbf{x}}}_{t-1}\mapsto{\boldsymbol{\mathbf{x}}}_{t}$, evaluated at ${\boldsymbol{\mathbf{x}}}_{t-1}$. We thus obtain

Writing $m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})$ as a finite sum over kernels and moving the integration into this sum, the integration in Eq. (37) turns into

for each state dimension $a=1,\dotsc,D$. With the SE covariance function $k_{\text{SE}}$ defined in Eq. (9), we compute the integral analytically and obtain

where we defined $c_{3}^{-1}=(\alpha_{f_{a}}^{2}(2\pi)^{\tfrac{D}{2}}\sqrt{|{\mathbf{\Lambda}}_{a}|})^{-1}$, such that $k_{f_{a}}({\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{x}}}_{i})=c_{3}\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t-1}\,|\,{\boldsymbol{\mathbf{x}}}_{i},{\mathbf{\Lambda}}_{a})$. In the definition of $c_{3}$, $\alpha_{f_{a}}^{2}$ is a hyper-parameter of $\mathcal{GP}_{f}$ responsible for the variance of the latent function in dimension $a$. Using the definition of ${\mathbf{S}}$ in Eq. (24), the product of the two Gaussians in Eq. (38) results in a new (unnormalized) Gaussian $c_{4}^{-1}\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t-1}\,|\,{\boldsymbol{\mathbf{\psi}}}_{i},{\mathbf{\Psi}})$ with

Pulling all constants outside the integral in Eq. (38), the integral determines the expected value of the product of the two Gaussians, ${\boldsymbol{\mathbf{\psi}}}_{i}$. For $a=1,\dotsc,D$, we obtain

Using $c_{3}c_{4}^{-1}=q_{a_{i}}^{x}$, see Eq. (23), and some matrix identities, we finally obtain

With the mean and the covariance of the joint distribution ${p}({\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$ given by Eqs. (22), (33), (39), and the filter step, all necessary components are provided to compute the smoothing distribution ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:T})$ analytically .

IV Simulations

In the following, we present results analyzing the robustness of state-of-the art nonlinear filters (Sec. IV-A) and the performances of the corresponding smoothers (Sec. IV-B).

We consider the nonlinear stochastic dynamic system

which is a modified version of the model used in . The system is modified in two ways: First, Eq. (40) does not contain a purely time-dependent term in the system, which would not allow for learning stationary transition dynamics. Second, we substituted a sinusoidal measurement function for the originally quadratic measurement function used by and . The sinusoidal measurement function increases the difficulty in computing the marginal distribution ${p}({\boldsymbol{\mathbf{z}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$ if the time update distribution ${p}({\boldsymbol{\mathbf{x}}}_{t}|{\boldsymbol{\mathbf{z}}}_{1:t-1})$ is fairly uncertain: While the quadratic measurement function can only lead to bimodal distributions (assuming a Gaussian input distribution), the sinusoidal measurement function in Eq. (41) can lead to an arbitrary number of modes—for a broad input distribution.

The prior variance was set to $\sigma_{0}^{2}=0.5^{2}$, i.e., the initial uncertainty was fairly high. The system and measurement noises (see Eqs. (40)–(41)) were relatively small considering the amplitudes of the system function and the measurement function. For the numerical analysis, a linear grid in the interval $$of mean values$(\mu_{0}^{x})_{i}$,$i=1,\dotsc,100$, was defined. Then, a single latent (initial) state$x_{0}^{(i)}$was sampled from${p}(x_{0}^{(i)})=\mathcal{N}((\mu_{0}^{x})_{i},\sigma_{0}^{2})$,$i=1,\dotsc,100$.

For the dynamic system in Eqs. (40)–(41), we analyzed the robustness in a single filter step of the EKF, the UKF, the CKF, an SIR PF (sequential importance resampling particle filter) with 200 particles, the GP-UKF, and the GP-ADF against the ground truth, closely approximated by the Gibbs-filter . Compared to the evaluation of longer trajectories, evaluating a single filter step makes it easier to analyze the robustness of individual filtering algorithms.

Tab. LABEL:tab:evaluation_nonlinear_system summarizes the expected performances (RMSE: root-mean-square error, MAE: mean-absolute error, NLL: negative log-likelihood) of the EKF, the UKF, the CKF, the GP-UKF, the GP-ADF, the Gibbs-filter, and the SIR PF for estimating the latent state $x$. The results in the table are based on averages over 1,000 test runs and 100 randomly sampled start states per test run (see experimental setup). The table also reports the 95% standard error of the expected performances. The $\star$ indicates a method developed in this paper. Tab. LABEL:tab:evaluation_nonlinear_system indicates that the GP-ADF is the most robust filter and statistically significantly outperforms all filters but the sampling-based Gibbs-filter and the SIR PF. The green color highlights a near-optimal Gaussian filter (Gibbs-filter) and the near-optimal particle filter. Amongst all other filters the GP-ADF is the closest Gaussian filter to the computationally expensive Gibbs-filter . Note that the SIR PF is not a Gaussian filter and is able to express multi-modality in distributions. Therefore, its performance is typically better than the one of Gaussian filters. The difference between the SIR PF and a near-optimal Gaussian filter, the Gibbs-filter, is expressed in Tab. LABEL:tab:evaluation_nonlinear_system. The performance difference essentially depicts how much we lose by using a Gaussian filter instead of a particle filter. The NLL values for the SIR PF are obtained by moment-matching the particles.

The poor performance of the EKF is due to linearization errors. The filters based on small sample approximations of densities (UKF, GP-UKF, CKF) suffer from the degeneracy of these approximations, which is illustrated in Fig. 1. Note that the CKF uses a smaller set of cubature points than the UKF to determine predictive distributions, which makes the CKF statistically even less robust than the UKF.

IV-B Smoother Robustness

where $I$ is the moment of inertia and $g$ the acceleration of gravity. Then, the successor state

Note that the scalar measurement Eq. (44) solely depends on the angle. Thus, the full distribution of the latent state ${\boldsymbol{\mathbf{x}}}$ had to be reconstructed using the cross-correlation information between the angle and the angular velocity.

We experimented with even smaller signal-to-noise ratios. The GP-RTSS remains robust, while the other smoothers remain unstable.

V Discussion and Conclusion

In this paper, we presented GP-RTSS, an analytic Rauch-Tung-Striebel smoother for GP dynamic systems, where the GPs with SE covariance functions are practical implementations of universal function approximators. We showed that the GP-RTSS is more robust to nonlinearities than state-of-the-art smoothers. There are two main reasons for this: First, the GP-RTSS relies neither on linearization (EKS) nor on density approximations (URTSS/CKS) to compute an optimal Gaussian approximation of the predictive distribution when mapping a Gaussian distribution through a nonlinear function. This property avoids incoherent estimates of the filtering and smoothing distributions as discussed in Sec IV-A. Second, GPs allow for more robust “system identification” than standard methods since they coherently represent uncertainties about the system and measurement functions at locations that have not been encountered in the data collection phase. The GP-RTSS is a robust smoother since it accounts for model uncertainties in a principled Bayesian way.

After training the GPs, which can be performed offline, the computational complexity of the GP-RTSS (including filtering) is $\mathcal{O}(T(E^{3}+n^{2}(D^{3}+E^{3})))$ for a time series of length $T$. Here, $n$ is the size of the GP training sets, and $D$ and $E$ are the dimensions of the state and the measurements, respectively. The computational complexity is due to the inversion of the $D$ and $E$-dimensional covariance matrices, and the computation of the matrix ${\mathbf{Q}}\in\mathds{R}^{n\times n}$ in Eq. (28), required for each entry of a $D$ and $E$-dimensional covariance matrix. The computational complexity scales linearly with the number of time steps. The computational demand of classical Gaussian smoothers, such as the URTSS and the EKS is $\mathcal{O}(T(D^{3}+E^{3}))$. Although not reported here, we verified the computational complexity experimentally. Approximating the online computations of the GP-RTSS by numerical integration or grids scales poorly with increasing dimension. These problems already appear in the histogram filter . By explicitly providing equations for the solution of the involved integrals, we show that numerical integration is not necessary and the GP-RTSS is a practical approach to filtering in GP dynamic systems.

Although the GP-RTSS is computationally more involved than the URTSS, the EKS, and the CKS, this does not necessarily imply that smoothing with the GP-RTSS is slower: function evaluations, which are heavily used by the EKS/CKS/URTSS are not necessary in the GP-RTSS (after training). In the pendulum example, repeatedly calling the ODE solver caused the EKS/CKS/URTSS to be slower than the GP-RTSS (with 250 training points) by a factor of two.

The increasing use of GPs for model learning in robotics and control will eventually require principled smoothing methods for GP models. To our best knowledge, the proposed GP-RTSS is the most principled GP-smoother since all computations can be performed analytically exactly, i.e., without function linearization or sigma/cubature point representation of densities, while exactly integrating out the model uncertainty induced by the GP distribution.

Code will be made publicly available at http://mloss.org.

Acknowledgements

This work was partially supported by ONR MURI grant N00014-09-1-1052, by Intel Labs, and by DataPath, Inc.

References