The Use of Gaussian Processes in System Identification

TL;DR

This paper reviews how Gaussian processes are applied in system identification, including modeling techniques like NFIR, NARX, and state-space models, highlighting their role in learning dynamic input-output relationships.

Contribution

It provides a concise overview of Gaussian process-based methods in system identification, emphasizing recent developments and applications.

Findings

Gaussian processes enable flexible nonlinear modeling.

They facilitate probabilistic predictions with uncertainty quantification.

Various models like NFIR, NARX, and state-space are discussed.

Abstract

Gaussian processes are used in machine learning to learn input-output mappings from observed data. Gaussian process regression is based on imposing a Gaussian process prior on the unknown regressor function and statistically conditioning it on the observed data. In system identification, Gaussian processes are used to form time series prediction models such as non-linear finite-impulse response (NFIR) models as well as non-linear autoregressive (NARX) models. Gaussian process state-space models (GPSS) can be used to learn the dynamic and measurement models for a state-space representation of the input-output data. Temporal and spatio-temporal Gaussian processes can be directly used to form regressor on the data in the time domain. The aim of this article is to briefly outline the main directions in system identification methods using Gaussian processes.

Figures

Equations

$$f(\mathbf{z})$$

$$y_k$$

$$k(\mathbf{z},{\mathbf{z}}^{\prime })=s^2\exp \left(-\frac{\Vert \mathbf{z}-{\mathbf{z}}^{\prime }{\Vert }^2}{2{\ell }^2}\right)$$

$$\left(\begin{array}{c}\mathbf{y}\\ f({\mathbf{Z}}^{\ast })\end{array}\right)\sim \mathcal{N}\left(\left(\begin{array}{c}m(\mathbf{Z})\\ m({\mathbf{Z}}^{\ast })\end{array}\right),\left(\begin{array}{cc}\left(\mathbf{K}+{\sigma }_{\mathrm{n}}^{2}I\right)& {\mathbf{k}}^{\top }({\mathbf{Z}}^{\ast })\\ \mathbf{k}({\mathbf{Z}}^{\ast })& k({\mathbf{Z}}^{\ast },{\mathbf{Z}}^{\ast })\end{array}\right)\right),$$

$$\begin{array}{rl}\mathrm{E}[f({\mathbf{Z}}^{\ast })\mid \mathbf{y}]& =m({\mathbf{Z}}^{\ast })\\ & +\mathbf{k}({\mathbf{Z}}^{\ast }){\left[\mathbf{K}+{\sigma }_{\mathrm{n}}^2\mathbf{I}\right]}^{-1}(\mathbf{y}-m(\mathbf{Z})),\\ \mathrm{Cov}[f({\mathbf{Z}}^{\ast })\mid \mathbf{y}]& =k({\mathbf{Z}}^{\ast },{\mathbf{Z}}^{\ast })\\ & -\mathbf{k}({\mathbf{Z}}^{\ast }){\left[\mathbf{K}+{\sigma }_{\mathrm{n}}^2\mathbf{I}\right]}^{-1}{\mathbf{k}}^{\top }({\mathbf{Z}}^{\ast }).\end{array}$$

$$\begin{array}{rl}-\log p(\mathbf{y}\mid \mathbf{\phi })& =\frac{1}{2}\log |2\pi ({\mathbf{K}}_{\mathbf{\phi }}+{\sigma }_{\mathrm{n}}^2\mathbf{I})|\\ & +\frac{1}{2}(\mathbf{y}-m_{\mathbf{\phi }}(\mathbf{Z}){)}^{\top }{\left[{\mathbf{K}}_{\mathbf{\phi }}+{\sigma }_{\mathrm{n}}^2\mathbf{I}\right]}^{-1}\\ & \times (\mathbf{y}-m_{\mathbf{\phi }}(\mathbf{Z})).\end{array}$$

$$p(\mathbf{\phi }\mid \mathbf{y})=\frac{p(\mathbf{y}\mid \mathbf{\phi })p(\mathbf{\phi })}{\int p(\mathbf{y}\mid \mathbf{\phi })p(\mathbf{\phi })\mathrm{d}\mathbf{\phi }},$$

$${\hat{y}}_k=f(u_{k-1},\dots ,u_{k-m}),$$

$$y_k$$

$$y_k=f(y_{k-1},\dots ,y_{k-n},u_{k-1},\dots ,u_{k-m})+{ϵ}_k,$$

$$y_k=f({\hat{y}}_{k-1},\dots ,{\hat{y}}_{k-n},u_{k-1},\dots ,u_{k-m})+{ϵ}_k,$$

$${\mathbf{x}}_{k+1}$$

$$y_k$$

$${\tilde{\mathbf{x}}}_{k+1}$$

$$y_k$$

$${\mathbf{x}}_{k+1}$$

$$y_k$$

$$p({\mathbf{f}}_k\mid {\mathbf{f}}_{1:k-1},{\mathbf{x}}_{0:k-1})$$

$$\begin{array}{rl}p({\mathbf{x}}_{0:N})& =\prod _{k=1}^N\mathcal{N}({\mathbf{x}}_k\mid {\mathbf{\mu }}_k({\mathbf{x}}_{0:k-1}),{\mathbf{\Sigma }}_k({\mathbf{x}}_{0:k-1})),\end{array}$$

$$\begin{array}{rl}& p(f({\mathbf{x}}^{\ast })\mid {\mathbf{x}}^{\ast },y_{1:N})\\ & =\int p(f({\mathbf{x}}^{\ast })\mid {\mathbf{x}}^{\ast },{\mathbf{x}}_{0:N})p({\mathbf{x}}_{0:N}\mid y_{1:N})\mathrm{d}{\mathbf{x}}_{0:N},\end{array}$$

$$f(\mathbf{x},\mathbf{u})=\sum _{i=1}^S{c}_i{\varphi }_i(\mathbf{x},\mathbf{u}),$$

$${\mathbf{x}}_{k+1}$$

$$y_k$$

$$\begin{array}{rl}& p(f({\mathbf{x}}^{\ast })\mid y_{1:N})\\ & =\int p(f({\mathbf{x}}^{\ast })\mid {\mathbf{x}}^{\ast },{\mathbf{x}}_{0:N})p({\mathbf{x}}_{0:N}\mid y_{1:N})\mathrm{d}{\mathbf{x}}_{0:N}\\ & \approx \int p(f({\mathbf{x}}^{\ast })\mid {\mathbf{x}}^{\ast },{\mathbf{f}}^u)p({\mathbf{f}}^u\mid y_{1:N})\mathrm{d}{\mathbf{f}}^u.\end{array}$$

$$\begin{array}{rl}& p(f({\mathbf{x}}^{\ast })\mid y_{1:N})\\ & \approx \int p(f({\mathbf{x}}^{\ast })\mid {\mathbf{x}}^{\ast },{\mathbf{f}}^u)q({\mathbf{f}}^u)\mathrm{d}{\mathbf{f}}^u\end{array}$$

$$\begin{array}{rl}f(t)& \sim \mathcal{G}\mathcal{P}(m(t),k(t,t^{\prime })),\\ y_k& =f(t_k)+{ϵ}_k.\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}\mathbf{x}(t)}{\mathrm{d}t}& =\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{\eta }(t),\\ y_k& =\mathbf{C}\mathbf{x}(t_k)+{ϵ}_k,\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}\mathbf{x}(\mathbf{z},t)}{\mathrm{d}t}& =\mathcal{A}\mathbf{x}(\mathbf{z},t)+\mathbf{B}\mathbf{\eta }(\mathbf{z},t),\\ y_k& =\mathcal{C}\mathbf{x}(\mathbf{z},t_k)+{ϵ}_k,\end{array}$$

$$\begin{array}{r}\frac{{\mathrm{d}}^{2}x(t)}{\mathrm{d}{t}^2}+\gamma \frac{\mathrm{d}x(t)}{\mathrm{d}t}+{\nu }^2=\xi (t),\end{array}$$

Full text

The Use of Gaussian Processes in System Identification

Simo Särkkä

(To appear in Encyclopedia of systems and control, 2nd edition)

1 Abstract

Gaussian processes are used in machine learning to learn input-output mappings from observed data. Gaussian process regression is based on imposing a Gaussian process prior on the unknown regressor function and statistically conditioning it on the observed data. In system identification, Gaussian processes are used to form time series prediction models such as non-linear finite-impulse response (NFIR) models as well as non-linear autoregressive (NARX) models. Gaussian process state-space models (GPSS) can be used to learn the dynamic and measurement models for a state-space representation of the input-output data. Temporal and spatio-temporal Gaussian processes can be directly used to form regressor on the data in the time domain. The aim of this article is to briefly outline the main directions in system identification methods using Gaussian processes.

2 Keywords

Gaussian process regression, non-linear system identification, GP-NFIR model, GP-ARX model, GP-NOE model, Gaussian process state-space model, temporal Gaussian process, state-space Gaussian process

3 Introduction

Gaussian process regression (Rasmussen and Williams, 2006) refers to a statistical methodology where we use Gaussian processes as prior models for regression functions that we fit to observed data. This kind of methodology is particularly popular in machine learning although the origins of the basic ideas can be traced to geostatistics (Cressie, 1993). In geostatistics, the corresponding methodology is called ”kriging” which is named after South African mining engineer D. G. Krige. In system identification, Gaussian processes can be used to identify (or ”learn” in machine learning terms) the input-output relationships from observed data. Even when there is no explicit input in the system, Gaussian processes can be used to identify a model for an observed time series of outputs which is a specific form of a system identification problem. Overviews of the use of Gaussian processes in system identification can be found, for example, in the monograph of Kocijan (2016), and PhD theses of McHutchon (2015) and Frigola (2016).

4 Gaussian processes in system identification

4.1 Gaussian process regression

Gaussian process regression is concerned with the following problem: Given a set of observed (training) input-output data $\mathcal{D}=\{(\mathbf{z}_{k},y_{k}):k=1,\ldots,N\}$ from an unknown function $y=f(\mathbf{z})$, predict the values of the function at new (test) inputs $\{\mathbf{z}_{k}^{*}:k=1,\ldots,M\}$. That is, the problem is a classical regression problem. However, the classical solution to the problem usually amounts to fixing a parametric function class $f(\mathbf{z};\boldsymbol{\theta})$, where $\boldsymbol{\theta}$ is a set of parameters and then fitting the parameters to the observed data. In Gaussian process regression, we take a different route; instead of fixing a parametric class of functions, we put a Gaussian process prior measure on the whole regression function and condition on the observed data using Bayes’ rule (Rasmussen and Williams, 2006).

4.1.1 Gaussian process regression problem

Mathematically, the Gaussian process regression problem can be written as

[TABLE]

where Equation (1a) tells that, a priori, the function is a Gaussian process with mean function $m(\mathbf{z})=\mathrm{E}[f(\mathbf{z})]$ and covariance function (or kernel) $k(\mathbf{z},\mathbf{z}^{\prime})=\mathrm{Cov}[f(\mathbf{z}),f(\mathbf{z^{\prime}})]=\mathrm{E}[(f(\mathbf{z})-m(\mathbf{z}))\,(f(\mathbf{z}^{\prime})-m(\mathbf{z}^{\prime}))]$. Equation (1b) tells that we observe the function values at points $\mathbf{z}_{k}$, $k=1,\ldots,N$ and that they are corrupted by (independent) Gaussian noises with variance $\sigma_{\mathrm{n}}^{2}$.

The mean and covariance functions define the regressor function class and they, or at least their parametric classes, need to be selected a priori. The mean function can typically be selected to be identically zero $m(\mathbf{z})=0$. The covariance function defines the smoothness properties of the functions, and a typical choice in machine learning is the squared exponential covariance function

[TABLE]

which produces infinitely differentiable (i.e., analytic) regressor functions. The parameters $s$ and $\ell$ in the aforementioned covariance function define the magnitude and length scales of the regressor functions, respectively. Other common choices of covariance functions are, for example, the Matérn class of covariance functions (Matérn, 1960; Rasmussen and Williams, 2006).

4.1.2 Gaussian process regression equations

Given the mean and covariance functions as well as the measurements, we can form the Gaussian process regressor. Assuming that the noises are independent of the function values, we can write the joint distribution of the observed values and the unknown function values as follows:

[TABLE]

where $\mathbf{y}=\begin{pmatrix}y_{1}&\ldots&y_{N}\end{pmatrix}^{\top}$, $m(\mathbf{Z})=\begin{pmatrix}m(\mathbf{z}_{1})&\cdots&m(\mathbf{z}_{N})\end{pmatrix}^{\top}$, $m(\mathbf{Z}^{*})=\begin{pmatrix}m(\mathbf{z}_{1}^{*})&\ldots&m(\mathbf{z}_{M}^{*})\end{pmatrix}^{\top}$, and $\mathbf{K}$ and $\mathbf{k}(\mathbf{Z}^{*})$ denote matrices with element $(i,j)$ given as $k(\mathbf{z}_{i},\mathbf{z}_{j})$ and $k(\mathbf{z}_{i}^{*},\mathbf{z}_{j})$, respectively. By conditioning this joint Gaussian distribution on the measurements $\mathbf{y}$ we get that the conditional (i.e., posterior) distribution of the function values $f(\mathbf{Z}^{*})=\begin{pmatrix}f(\mathbf{z}_{1}^{*})&\ldots&f(\mathbf{z}_{M}^{*})\end{pmatrix}^{\top}$ is Gaussian with the mean and covariance

[TABLE]

These are the fundamental equations of Gaussian process regression. An example of Gaussian process regression with squared exponential covariance function is shown in Figure 1.

4.1.3 Hyperparameter learning

Even though Gaussian process regression is a non-parametric method, for which we do not need to fix a parametric class of functions, the mean and covariance functions can have unknown hyperparameters $\boldsymbol{\varphi}$ which can be estimated from data. For example, the squared exponential covariance function in Equation (2) has the hyperparameters $\boldsymbol{\varphi}=(s,\ell)$.

A common way to estimate the parameters is to maximize the marginal likelihood – also called evidence – $p(\mathbf{y}\mid\boldsymbol{\varphi})$ of the measurements, or equivalently, minimize the negative log-likelihood of the measurements

[TABLE]

The gradient of this function with respect to the hyperparameters is also available (see, e.g., Rasmussen and Williams, 2006) which allows for the use of gradient-based optimization methods to estimate the parameters.

Instead of using the maximum likelihood method to estimate the parameters, it is also possible to use a Bayesian approach to the problem and consider the posterior distribution of the hyperparameters

[TABLE]

where $p(\boldsymbol{\varphi})$ is the prior distribution of the hyperparameters. We can, for example, compute the maximum a posteriori estimate of the parameters by finding the maximum of this distribution or use Markov chain Monte Carlo (MCMC) methods (Brooks et al, 2011) to estimate the statistics of the distribution.

In what follows, to avoid notational clutter, we drop out the hyperparameters from the Gaussian process formulations and inference methods although they are commonly estimated as part of the Gaussian process learning.

4.1.4 Reduction of computational complexity

A limitation of Gaussian process regression in its explicit form is that the computational complexities of the regression Equations (4) and likelihood Equation (5) are cubic $\mathcal{O}(N^{3})$ in the number of measurements $N$. This is due to the $N\times N$ matrix inversion appearing in the equations which, even when implemented with Cholesky or LU decompositions, needs a cubic number of computational steps.

Ways of solving the computational complexity problem are, for example, sparse approximations using inducing points (Quiñonero-Candela and Rasmussen, 2005; Rasmussen and Williams, 2006; Titsias, 2009), approximating the problem with a discrete Gaussian random field model (Lindgren et al, 2011), or use of random or deterministic basis/spectral expansions (Quiñonero-Candela et al, 2010; Solin and Särkkä, 2018).

4.2 GP-NFIR, GP-NARX, GP-NOE, and related models

In system identification, we can use Gaussian processes to model unknown input-output relationships in time series. Several different model architectures are available for this purpose. Let us assume that we have a system with input sequence $u_{1},u_{2},\ldots$ and output sequence $y_{1},y_{2},\ldots$ and the aim is to predict the outputs from inputs. We also assume that we have been given a set of training data consisting of known inputs and (noisy) outputs. In the following, we present some typically used architectures that have been proposed for this purpose. More details can be found in the monograph of Kocijan (2016).

4.2.1 GP-NFIR model

The Gaussian process non-linear finite impulse response (GP-NFIR) model (Ackermann et al, 2011; Kocijan, 2016) has the form (see Figure 2)

[TABLE]

where $f(\cdot)$ is an unknown mapping which we model as a Gaussian process, and $\hat{y}_{k}$ denotes the estimate produced by the regressor. In this model, we form a Gaussian process regressor that predicts the current output from a finite number of previous inputs. This model can be identified by reducing it into a Gaussian process regression model

[TABLE]

where $\mathbf{z}_{k}=\begin{pmatrix}u_{k-1}&\ldots&u_{k-m}\end{pmatrix}^{\top}$ and by using standard Gaussian process regression methods on it.

4.2.2 GP-NARX model

The Gaussian process nonlinear autoregressive model with exogenous input (GP-NARX) (Kocijan et al, 2005; Kocijan, 2016) is a model of the form (see Figure 3)

[TABLE]

where $\epsilon_{k}$ is a Gaussian random variable. This model can be reduced to a Gaussian process regression problem by setting $\mathbf{z}_{k}=\begin{pmatrix}y_{k-1}&\cdots&y_{k-n}&u_{k-1}&\cdots&u_{k-m}\end{pmatrix}^{\top}$ in Equation (8a).

4.2.3 GP-NOE model

In Gaussian process nonlinear output error (GP-NOE) model (Kocijan and Petelin, 2011; Kocijan, 2016) we form a Gaussian process regressor for the problem (see Figure 4)

[TABLE]

where $\hat{y}_{k-1},\ldots,\hat{y}_{k-n}$ are the Gaussian process regressor predictions from the previous steps.

Learning in this kind of model requires further approximations because the predictions of the Gaussian process are directly used as inputs on the next step.

4.2.4 Other model architectures

As discussed in Kocijan (2016), it is also possible extend these architectures to, for example, GP-NARMAX (nonlinear autoregressive and moving average model with exogenous input) models and NJB (nonlinear Box–Jenkins) models.

4.3 Gaussian process state-space (GPSS) models

Another approach to system identification is to form a state-space model where the dynamic and measurement models are identified using Gaussian process regression methods. This leads to so-called Gaussian process state-space models.

4.3.1 General GPSS model

A Gaussian process state-space (GPSS) model (see Figure 5) has the mathematical form (e.g. Kocijan, 2016)

[TABLE]

where the state vector $\mathbf{x}_{k}$, $k=0,1,2,\ldots,N$ contains the current state of the system, $u_{1},u_{2},\ldots$ is the input sequence and $y_{1},y_{2},\ldots$ is the output sequence. In the model, $\mathbf{w}_{k}$ is a Gaussian distributed process noise. The aim is now to learn the functions $f(\mathbf{x}_{k},u_{k})$ and $g(\mathbf{x}_{k},u_{k})$, which are modeled as Gaussian processes, given the input and output sequences, or in some cases, also given direct observations of the state vector.

4.3.2 Learning with fully observed state

When the state vector $\mathbf{x}_{k}$ is fully observed, then both the dynamic model (11a) and measurement model (11b) become standard Gaussian process regression models. In dynamic model (11a), the training set consists of measurements $\mathbf{x}_{k+1}$ with the corresponding inputs $(\mathbf{x}_{k},u_{k})$, and in measurement model (11b) the measurements are $y_{k}$ with the corresponding inputs $(\mathbf{x}_{k},u_{k})$. This kind of fully observed models are important in many applications such as robotics (Deisenroth et al, 2015).

After conditioning on the training data, the functions $f$ and $g$ will still be Gaussian processes and their mean and covariance functions are given by (multivariate generalizations) of Equations (4). State estimation in this kind of models has been considered by Ko and Fox (2009) and Deisenroth et al (2011), and it turns out that it is possible to construct closed-form Gaussian approximation (moment matching) based filters and smoothers for these models (Deisenroth et al, 2011). Control problems related to this kind of models have been considered, for example, by Deisenroth et al (2015).

4.3.3 Marginalization of the GP

When the states $\mathbf{x}_{k}$ are not observed, then we need to treat both the states and the Gaussian processes as unknown. There are a few different ways to cope with the model in that case, and one approach is the marginalization approach of Frigola et al (2013, 2014b, 2014a) and Frigola (2016). First note that if we have a method to learn $f$, we can learn both $f$ and $g$ using a state-augmentation trick (Frigola, 2016): we define an augmented state as $\tilde{\mathbf{x}}=\begin{pmatrix}\mathbf{x}&\gamma\end{pmatrix}^{\top}$, Gaussian process $h(\boldsymbol{\gamma},u)=\begin{pmatrix}f(\mathbf{x},u)&g(\mathbf{x},u)\end{pmatrix}^{\top}$, and augmented process noise $\tilde{\mathbf{w}}=\begin{pmatrix}\mathbf{w}_{k}&\mathbf{0}\end{pmatrix}^{\top}$, which reduces the model to

[TABLE]

This is now a model with an unknown dynamic model, but with a given linear Gaussian measurement model $p(y_{k}\mid\tilde{\mathbf{x}}_{k},u_{k})=\mathcal{N}(y_{k}\mid\gamma_{k},\sigma^{2}_{\mathrm{n}})$.

Thus, without a loss of generality, we can focus on models with an unknown dynamic model $f$, and a known measurement model:

[TABLE]

The aim is to learn the function $f(\mathbf{x})$ from a sequence of measurement data $y_{k}$.

One way to understand (Frigola, 2016) this model is that we could hypothetically generate data from it by sampling the (infinite-dimensional) function $f(\mathbf{x},u)$, and then starting from $\mathbf{x}_{0}$ sequentially produce each $\{\mathbf{f}_{1},\mathbf{x}_{1},\ldots,\mathbf{f}_{N},\mathbf{x}_{N}\}$, where we have denoted $\mathbf{f}_{k}=f(\mathbf{x}_{k-1},u_{k-1})$ for $k=1,2,\ldots,N$. Each of the conditional distributions

[TABLE]

turns out to be Gaussian. Note that above, we have introduced the short-hand notation $\mathbf{f}_{1:k}=(\mathbf{f}_{1},\ldots,\mathbf{f}_{k})$ which we will use also in the rest of this article.

The above observation allows us to integrate out (i.e., marginalize) the Gaussian process from the model in closed form. The result is the following representation (Frigola et al, 2013):

[TABLE]

where the means $\boldsymbol{\mu}_{k}(\mathbf{x}_{0:k-1})$ and covariances $\boldsymbol{\Sigma}_{k}(\mathbf{x}_{0:k-1})$ are (quite complicated) functions of the prior mean and covariance functions of the Gaussian process $f$, which are evaluated on the whole previous state history. The above equation defines a non-Markovian prior model for the state sequence $\mathbf{x}_{k}$, $k=1,\ldots,N$.

For a given $\mathbf{x}_{0:N}$, the distribution $p(f(\mathbf{x}^{*})\mid\mathbf{x}^{*},\mathbf{x}_{0:N})$ for a test point $\mathbf{x}^{*}$ can be computed by using conventional Gaussian process prediction equation as follow:

[TABLE]

which can be numerically approximated, provided that we use a convenient (e.g. Monte Carlo) approximation for $p(\mathbf{x}_{0:N}\mid y_{1:N})$.

Given the model (15) it is then possible to use, for example, particle Markov chain Monte Carlo methods (Frigola et al, 2013) to sample state trajectories from the posterior distribution $p(\mathbf{x}_{0:N}\mid y_{1:N})$ jointly with the parameters of the model, which provides a Monte Carlo approximation to the above integral. Other proposed approaches are, for example, particle stochastic approximation expectation–maximization (EM, Frigola et al, 2014b) which uses a Monte Carlo approximation to the EM algorithm aiming at computing the maximum likelihood estimates of the parameters while handling the states as missing data.

4.3.4 Approximation of the GP

Another way to approach the problem where both the states and Gaussian processes are unknown is to approximate the Gaussian process as finite-dimensional parametric model and use conventional parameter estimation methods to the model.

One possible approximation considered in Svensson et al (2016) is to employ a Karhunen–Loeve type of basis function expansion of the Gaussian process as follows:

[TABLE]

where $\phi_{i}(\mathbf{x},\mathbf{u})$ are deterministic basis functions (e.g. sinusoidals) and $c_{i}$ are Gaussian random variables. With this approximation, the model in Equation (13) becomes

[TABLE]

where learning of the Gaussian process $f$ reduces to estimation of the finite number of parameters $\mathbf{c}=\begin{pmatrix}c_{1}&\cdots&c_{S}\end{pmatrix}^{\top}$ in the state-space model. The states and parameters in this model can now be determined using, for example, particle Markov chain Monte Carlo (PMCMC) methods (see Svensson et al, 2016).

Another possibility is to use inducing points (typically denoted with $\mathbf{u}$, but here we denote them with $\mathbf{f}^{u}$ to avoid confusion with the input sequence). In those approaches the idea is to first perform Gaussian process inference on the inducing points alone, that is, compute $p(\mathbf{f}^{u}\mid y_{1:N})$ and then compute the (approximate) predictions by conditioning on the inducing points instead of the original data. This leads to the approximation

[TABLE]

In the Turner et al (2010) the inducing points are learned using expectation–maximization (EM) algorithm. Frigola et al (2014a) propose a method for variational learning (or integration over) the inducing points by forming a variational Bayesian approximation $q(\mathbf{f}^{u})\approx p(\mathbf{f}^{u}\mid y_{1:N})$, which further results in

[TABLE]

and turns out to be analytically tractable as the optimal variational distribution $q(\mathbf{f}^{u})$ is Gaussian.

4.4 Spatio-temporal Gaussian process models

4.4.1 Temporal Gaussian processes

Another way of modeling time series using Gaussian processes is by considering them as functions of time (Hartikainen and Särkkä, 2010) which are sampled at certain time instants $t_{1},t_{2},\ldots$:

[TABLE]

That is, instead of attempting to form a predictor from the previous measurements or inputs, the idea is to condition the temporal Gaussian process on its observed values and use the conditional Gaussian process for predicting values at new time points.

Unfortunately, due to the cubic computational scaling of the Gaussian process regression, this quickly becomes in tractable when time series length increases. However, temporal Gaussian process regression is closely related to the classical Kalman filtering and Rauch-Tung-Striebel smoothing problems (e.g. Särkkä, 2013), which can be used to reduce the required computations. It turns out that provided that the Gaussian process is stationary, that is, the covariance function only depends on the time difference $k(t,t^{\prime})=k(t-t^{\prime})$, then, under certain restrictions, the Gaussian process regression problem is essentially equivalent to state-estimation in a model of the form

[TABLE]

where $\boldsymbol{\eta}(t)$ is a white noise process and the matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are selected suitably to match the original covariance function. For example, the Matérn covariance functions with half-integer smoothness parameters have exact representations as state-space models (Hartikainen and Särkkä, 2010).

The Gaussian process regression problem can be now solved by applying a Kalman filter and Rauch-Tung-Striebel smoother on this problem. These methods have the fortunate property, that their complexity is linear $\mathcal{O}(N)$ with respect to the number of measurements $N$ as opposed to the cubic complexity of direct Gaussian process regression solution.

4.4.2 Spatio-temporal Gaussian processes

A similar state-space approach also works for spatio-temporal Gaussian process models with a covariance function of a stationary form $k(\mathbf{z},t;\mathbf{z}^{\prime},t^{\prime})=k(\mathbf{z},\mathbf{z}^{\prime};t-t^{\prime})$. In that case, the state-estimation problem becomes infinite-dimensional, that is, a distributed parameter system

[TABLE]

where $\mathcal{A}$ is a matrix of operators and $\mathcal{C}$ is a matrix of functionals (Särkkä et al, 2013). The solution of the Gaussian process regression problem in this form requires the use of methods from partial differential equations, but in many cases we can obtain an exact $\mathcal{O}(N)$ inference procedure from this route.

4.4.3 Latent force models

In so-called latent force models (Álvarez et al, 2013) the idea is to infer latent force $\xi(t)$ in a differential equation model such as

[TABLE]

where $\xi(t)$ is an unknown function which is modeled as a Gaussian process. The inference in this model can be recast as Gaussian process regression with a modified covariance function. The idea can be further generalized to partial differential equation models and non-linear models when approximation methods such as Laplace approximation are used.

The inference in this kind of models can be further re-stated in the state-space form (Särkkä et al, 2019), which also allows for the study of control problems on latent force models. This formulation also allows the analysis of observability and controllability properties of latent force models.

5 Summary and Future Directions

In this article, we have briefly outlined the main directions in system identification using Gaussian processes. The methods can be divided into three classes: (1) GP-NFIR, GP-NARX, and GP-NOE type of models which directly aim at learning the function from the previous inputs and outputs to the current output, (2) Gaussian process state-space models which aim to learn the dynamic and measurement models in the state-space model, and (3) Gaussian process regression methods for the spatio-temporal time series by direct or state-space Gaussian process regression.

Active problems in research in this area appear to be, for example, joint learning and control in all the model types outlined here. Another direction of further study is the analysis of observability, identifiability, and controllability. In practice, the Gaussian process-based system identification methods are very similar to classical methods, and hence they can be expected to inherit many limitations and theoretical properties of the classical methods.

An emerging research area in Gaussian process-based models are so-called deep Gaussian processes (Damianou and Lawrence, 2013) which borrow the idea of deep neural networks by forming hierarchies of Gaussian processes. This kind of models could also turn out to be useful in system identification. Furthermore, Gaussian processes are also easily combined with first-principles-models (cf. latent force models described above) which allows for flexible gray box modeling using Gaussian processes.

One of the main obstacles in Gaussian process regression is still the computational scaling in the number of measurements, which is also inherited by all the new developments. Although several good approaches to tackle this problem have been proposed, the problem still is that they inherently replace the original model with an approximation. New better approaches to this problem are likely to appear in the near future.

6 Cross References