Recursive Estimation for Sparse Gaussian Process Regression

TL;DR

This paper introduces a recursive estimation method for sparse Gaussian Process regression that enables efficient online learning with mini-batches, achieving faster convergence and better performance on large datasets.

Contribution

It develops a novel recursive approach linking sparse GP methods with Bayesian recursive estimation, allowing analytic updates in online mini-batch training.

Findings

Faster convergence compared to existing methods

Superior performance on large-scale datasets

Effective online learning with recursive updates

Abstract

Gaussian Processes (GPs) are powerful kernelized methods for non-parameteric regression used in many applications. However, their use is limited to a few thousand of training samples due to their cubic time complexity. In order to scale GPs to larger datasets, several sparse approximations based on so-called inducing points have been proposed in the literature. In this work we investigate the connection between a general class of sparse inducing point GP regression methods and Bayesian recursive estimation which enables Kalman Filter like updating for online learning. The majority of previous work has focused on the batch setting, in particular for learning the model parameters and the position of the inducing points, here instead we focus on training with mini-batches. By exploiting the Kalman filter formulation, we propose a novel approach that estimates such parameters by recursively propagating the analytical gradients of the posterior over mini-batches of the data. Compared to state of the art methods, our method keeps analytic updates for the mean and covariance of the posterior, thus reducing drastically the size of the optimization problem. We show that our method achieves faster convergence and superior performance compared to state of the art sequential Gaussian Process regression on synthetic GP as well as real-world data with up to a million of data samples.

Equations

$$\begin{array}{rl}{\mathbf{\mu }}_{\ast }& ={\mathbf{K}}_{\ast \mathbf{X}}{\left({\mathbf{K}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)}^{-1}\mathbf{y},\\ {\mathbf{\Sigma }}_{\ast }& ={\mathbf{K}}_{\ast \ast }-{\mathbf{K}}_{\ast \mathbf{X}}{\left({\mathbf{K}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)}^{-1}{\mathbf{K}}_{\mathbf{X}\ast },\end{array}$$

$$\log p\left(\mathbf{y}|\mathbf{\varphi }\right)=\log \mathcal{N}\left(\mathbf{y}|0,{\mathbf{K}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right).$$

$$\begin{array}{rl}{\mathbf{\mu }}_{\ast }& ={\mathbf{Q}}_{\ast \mathbf{X}}{\left({\overline{\mathbf{K}}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)}^{-1}\mathbf{y},\\ {\mathbf{\Sigma }}_{\ast }& ={\mathbf{K}}_{\ast \ast }-{\mathbf{Q}}_{\ast \mathbf{X}}{\left({\overline{\mathbf{K}}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)}^{-1}{\mathbf{Q}}_{\mathbf{X}\ast },\end{array}$$

$$\begin{array}{rl}{\mathcal{L}}_{PEP}(\mathbf{\theta })=& \log \mathcal{N}\left(\mathbf{y}|0,{\overline{\mathbf{K}}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)\\ & -\frac{1-\alpha }{2\alpha }\sum _{i=1}^N\log \left(1+\frac{\alpha }{{\sigma }_n^2}{\left[{\mathbf{D}}_{\mathbf{X}}\right]}_{ii}\right),\end{array}$$

$${\mathcal{L}}_{VFE}(\mathbf{\theta })=\log \mathcal{N}\left(\mathbf{y}|0,{\mathbf{Q}}_{\mathbf{X}\mathbf{X}}+{\sigma }_n^2\mathbb{I}\right)-\frac{\text{Tr}\left[{\mathbf{D}}_{\mathbf{X}}\right]}{2{\sigma }_n^2}.$$

$$\begin{array}{rl}& {\mathcal{L}}_{SVGP}(\mathbf{\mu },\mathbf{\Sigma },\mathbf{\theta })=-\text{KL}\left[\mathbf{q}(\mathbf{u})||p\left(\mathbf{u}|\mathbf{\theta }\right)\right]\\ & +K\sum _{k=1}^K\int \mathbf{q}(\mathbf{u})p\left({\mathbf{f}}_k|\mathbf{u},\mathbf{\theta }\right)p\left({\mathbf{y}}_k|{\mathbf{f}}_k,\mathbf{\theta }\right)\mathrm{d}\mathbf{u},\end{array}$$

$$f\left(\mathbf{X}\right)$$

$${\mathbf{y}}_k$$

$${\mathbf{f}}_k={\mathbf{H}}_k\mathbf{u}+{\mathbf{\gamma }}_k$$

$$p\left({\mathbf{y}}_k|{\mathbf{f}}_k\right)$$

$$p\left(\mathbf{u}\right)$$

$$p\left({\mathbf{f}}_k|\mathbf{u}\right)$$

$$p\left({\mathbf{f}}_{\ast }|\mathbf{u}\right)$$

$$p\left(\mathbf{u}|\mathbf{y}\right)\propto p\left(\mathbf{y}|\mathbf{u}\right)p\left(\mathbf{u}\right)$$

$$\begin{array}{r}p\left(\mathbf{u}|{\mathbf{y}}_{1:k}\right)\propto p\left({\mathbf{y}}_k|\mathbf{u}\right)p\left(\mathbf{u}|{\mathbf{y}}_{1:k-1}\right)\propto \mathcal{N}\left(\mathbf{u}|{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k\right),\end{array}$$

$$\begin{array}{rl}{\mathbf{r}}_k& ={\mathbf{y}}_k-{\mathbf{H}}_k{\mathbf{\mu }}_{k-1};\\ {\mathbf{S}}_k& ={\mathbf{H}}_k{\mathbf{\Sigma }}_{k-1}{\mathbf{H}}_k^T+{\mathbf{V}}_k;\\ {\mathbf{G}}_k& ={\mathbf{\Sigma }}_{k-1}{\mathbf{H}}_k^T{\mathbf{S}}_k^{-1};\end{array}\begin{array}{rl}{\mathbf{\mu }}_k& ={\mathbf{\mu }}_{k-1}+{\mathbf{G}}_k{\mathbf{r}}_k;\\ {\mathbf{\Sigma }}_k& ={\mathbf{\Sigma }}_{k-1}-{\mathbf{G}}_k{\mathbf{S}}_k{\mathbf{G}}_k^T.\end{array}$$

$$\begin{array}{rl}p\left({\mathbf{f}}_{\ast }|{\mathbf{y}}_{1:k}\right)& =\int p\left({\mathbf{f}}_{\ast }|\mathbf{u}\right)p\left(\mathbf{u}|{\mathbf{y}}_{1:k}\right)\mathrm{d}\mathbf{u}\\ & =\mathcal{N}\left({\mathbf{f}}_{\ast }|{\mathbf{H}}_{\ast }{\mathbf{\mu }}_k,{\mathbf{H}}_{\ast }{\mathbf{\Sigma }}_k{\mathbf{H}}_{\ast }^T+{\mathbf{V}}_{\ast }\right)\end{array}$$

$${\mathbf{\mu }}_K^{\ast }$$

$${\mathbf{\Sigma }}_K^{\ast }$$

$$\begin{array}{rl}{\mathbf{\mu }}_K^{\ast }& ={\mathbf{Q}}_{\ast \mathbf{X}}{\left({\mathbf{Q}}_{\mathbf{X}\mathbf{X}}+\mathbf{V}\right)}^{-1}\mathbf{y};\\ {\mathbf{\Sigma }}_K^{\ast }& ={\mathbf{Q}}_{\ast \ast }-{\mathbf{Q}}_{\ast \mathbf{X}}{\left({\mathbf{Q}}_{\mathbf{X}\mathbf{X}}+\mathbf{V}\right)}^{-1}{\mathbf{Q}}_{\mathbf{X}\ast }+{\mathbf{V}}_{\ast }.\end{array}$$

$$\begin{array}{rl}\log p\left(\mathbf{y}\right)& =\log \int p\left(\mathbf{y}|\mathbf{u}\right)p\left(\mathbf{u}\right)\mathrm{d}\mathbf{u}\\ & =\log \mathcal{N}\left(\mathbf{y}|0,\mathbf{H}{\mathbf{\Sigma }}_0{\mathbf{H}}^T+\mathbf{V}\right).\end{array}$$

$$p\left({\mathbf{y}}_k|{\mathbf{y}}_{1:k-1}\right)$$

$$=\mathcal{N}\left({\mathbf{r}}_k|0,{\mathbf{S}}_k\right).$$

$$\begin{array}{rl}& \log \prod _{k=1}^{K}p\left({\mathbf{y}}_k|{\mathbf{y}}_{1:k-1}\right)=\sum _{k=1}^K\log \mathcal{N}\left({\mathbf{r}}_k|0,{\mathbf{S}}_k\right)\\ & =-\frac{N}{2}\log 2\pi -\frac{1}{2}\sum _{k=1}^K\log \left|{\mathbf{S}}_k\right|+{\mathbf{r}}_k^T{\mathbf{S}}_k^{-1}{\mathbf{r}}_k.\end{array}$$

$$\begin{array}{r}\psi (\mathbf{\theta })=\sum _{k=1}^K\log \mathcal{N}\left({\mathbf{r}}_k^{\mathbf{\theta }}|0,{\mathbf{S}}_k^{\mathbf{\theta }}\right)-a_k(\mathbf{\theta })\end{array}$$

$${\mathbf{\theta }}^{(t)}={\mathbf{\theta }}^{(t-1)}-{\gamma }^{(t-1)}{\frac{\partial {\psi }_k}{\partial \mathbf{\theta }}}_{|\mathbf{\theta }={\mathbf{\theta }}^{(t-1)}},$$

$${\psi }^{(K)}(\mathbf{\theta })=\sum _{k=1}^K{d}_k(\mathbf{\theta })-a_k(\mathbf{\theta })=\sum _{k=1}^K{\psi }_k(\mathbf{\theta })$$

$$\frac{\partial d_k(\mathbf{\theta })}{\partial \mathbf{\theta }}=-\frac{1}{2}\frac{\partial \log \left|{\mathbf{S}}_k^{\mathbf{\theta }}\right|}{\partial \mathbf{\theta }}-\frac{1}{2}\frac{\partial ({\mathbf{r}}_k^{\mathbf{\theta }}{)}^T({\mathbf{S}}_k^{\mathbf{\theta }}{)}^{-1}{\mathbf{r}}_k^{\mathbf{\theta }}}{\partial \mathbf{\theta }}$$

$$\frac{\partial {\mathbf{u}}_k}{\partial \theta }$$

$$\frac{\partial {\mathbf{\Sigma }}_k}{\partial \theta }$$

$$\frac{\partial {\psi }_k(\mathbf{\theta })}{\partial \theta }=\frac{\partial d_k(\mathbf{\theta })}{\partial \theta }-\frac{\partial a_k(\mathbf{\theta })}{\partial \theta }$$

$$p\left(\mathbf{u}|{\mathbf{y}}_{1:k},{\mathbf{\theta }}^{(e,k)}\right)\text{raise}1.07639ptto0.0pt\propto \text{hss}\text{lower}4.03563pt\sim p\left({\mathbf{y}}_k|\mathbf{u},{\mathbf{\theta }}^{(e,k)}\right)p\left(\mathbf{u}|{\mathbf{y}}_{1:k-1},{\mathbf{\theta }}^{(e,k-1)}\right).$$

$$\frac{d}{dt}h(t)=w_1(t)+w_2(t)-0.2\sqrt{h(t)}$$

$$\frac{d}{dt}{C}_b(t)=(C_{b1}-C_b(t))\frac{w_1(t)}{h(t)}+(C_{b2}-C_b(t))\frac{w_2(t)}{h(t)}-\frac{k_1{C}_b(t)}{(1+k_2{C}_b(t){)}^2},$$

Full text

, , ,

1 Introduction

Gaussian process (GPs) regression is used in many applications, ranging from machine learning, social sciences, natural sciences and engineering, due to its modeling flexibility, robustness to overfitting and availability of well-calibrated predictive uncertainty estimates. In control engineering, for example, GPs have been used in system identification for impulse response estimation [25, 9, 26, 24], nonlinear ARX models [19, 3], learning of ODEs [1, 21], latent force modeling [42] and to learn the state space of a nonlinear dynamical system [12, 23, 35]. However, GPs do not scale to large data sets due to their $\mathcal{O}(N^{2})$ memory and $\mathcal{O}(N^{3})$ computational costs, where $N$ is the number of training samples. For this reason several sparse GP approximations have been proposed in the literature. Often such approximations are based on inducing points methods, where the unknown function is represented by its values at a set of $M\ll N$ pseudo-inputs, called inducing points. Among such methods, Subset of Regressors (SoR/DIC) approximations [32, 39, 33] produces overconfident predictions when leaving the training data. On the other hand, Deterministic Training Conditional (DTC) [11, 31], Fully Independent Training Conditional (FITC) [34], Fully Independent Conditional (FIC) [27] and Partially Independent Training Conditional (PITC) [27] all produce sensible uncertainty estimates. These models differ from each other in the definition of their joint prior over the latent function and test values. Titsias [36], instead, proposed to retain the exact prior but to perform approximate (variational) inference for the posterior, leading to the Variational Free Energy (VFE) method which converges to full GP as $M$ increases. Bui et al. [7] introduced Power Expectation Propagation (PEP), based on the minimization of an $\alpha$-divergence, which unifies most of the previously mentioned models. Typically, inference is achieved in $\mathcal{O}(M^{2}N)$ time and $\mathcal{O}(MN)$ space. In order to find good parameters (inducing input points and kernel hyper-parameters), either the log marginal likelihood of the sparse models or a lower bound are numerically optimized.

The previously mentioned approximations focus on the batch setting, i.e., all data is available at once and can be processed together. For big data, where the number of samples can be millions, keeping all data in memory is not possible, moreover the data might even arrive sequentially. Bui et al. [6] developed an algorithm to update hyper-parameters in an online fashion promising in a streaming setting, but with limited accuracy as each sample is considered only once.

We focus here on the setting where hyper-parameters are learned by reconsidering mini-batches several times. In order to speed up the optimization, we would like to update the parameters more frequently for a subset of data and update the posterior in a sequential way. In this setting, Hensman et al. [15] applied Stochastic Variational Inference (SVI, [17]) to an uncollapsed lower bound of the marginal likelihood. The resulting Stochastic Variational Gaussian Process (SVGP) method allows to optimize the parameters with mini-batches. Although showing high scalability and good accuracy, SVGP has two main drawbacks: i) the (variational) posterior is not given analytically, which leads to $\mathcal{O}(M^{2})$ additional many parameters; ii) the uncollapsed bounds are in practice often less tight than the corresponding collapsed VFE batch bounds because the (variational) posterior is not optimally eliminated. The large number of parameters ($\approx MD+M^{2}$, where $D$ is the input space dimension) leads to a hard-to-tune optimization problem which requires appropriately decaying learning rates. Even for fixed parameters, each sample still needs to be reconsidered many times. An orthogonal direction was pursued by the authors in [14] and [30], where a connection between GPs and State Space models for particular kernels was established for spatio-temporal regression problems, which allows to apply sequential algorithm such as the Kalman Filter, see also [8, 2]. Inspired by this line of research, the authors in [37] focused on efficient implementation and extended the methodology to varying sampling locations over time. These approaches can deal with sequential data and solve the problem of temporal time complexity, however the space complexity is still cubic in $N$. In addition, the hyper-parameters are usually fixed in advance.

In this work we propose a recursive collapsed lower bound to the log marginal likelihood which can be optimized stochastically with mini-batches.

In this respect, the first contribution of this paper is the derivation of a novel Kalman-filter-like (KF) formulation for a generic sparse inducing point method. In particular we show that sparse inducing point models can be seen as a Bayesian kernelized linear regression model with input dependent observation noise, a particular choice of basis functions and noise covariance. Given the model hyper-parameters, KF allows to train sparse GP methods analytically and exactly in an online setting (considering each sample only once, as opposed to the work in [15] and [16]). In this formulation the posterior distribution obtained online is equivalent to full batch methods. This constitutes an interesting technique on its own for applications where hyper-parameters are given, however the analysis above provides a key insight for parameter estimation.

Our second main contribution is a recursive approach to hyper-parameter estimation based on the KF formulation. It is based on recursively exploiting the chain rule for derivatives by recursively propagating the analytical gradients of the posterior which enables us to compute the derivatives of the lower bound sequentially. We show that, when computing the gradients of the recursive collapsed bound in a non-stochastic way, they exactly match the corresponding batch ones. This new Stochastic Recursive Gradient Propagation (SRGP)111Code is available at https://github.com/manuelIDSIA/SRGP. approach constitutes an efficient method to train a very general class of sparse GP regression models with much fewer parameters to be estimated numerically ($\approx MD$) than state of the art sequential GP regression methods ($\approx MD+M^{2}$). Since the number $M$ of inducing points determines the quality of the approximation to full GP, this reduction in number of parameters from $M^{2}$ to $M$ is crucial and results in more accurate and faster convergence than state of the art approaches such as SVGP. For example, in the application to learn the input output behavior of a non-linear plan presented in Sect. 5 the number of parameters estimated by SVGP is $\approx 10500$ while our approach only estimates $\approx 500$ parameters due to the analytical updates.

2 Background on GP Regression

Consider a training set $\mathcal{D}=\left\{y_{i},\bm{x}_{i}\right\}_{i=1}^{N}$ of $N$ pairs of inputs $\bm{x}_{i}\in\mathbb{R}^{D}$ and noisy scalar outputs $y_{i}$ generated by adding independent Gaussian noise to a latent function $f(\bm{x})$, that is $y_{i}=f(\bm{x}_{i})+\varepsilon_{i}$, where $\varepsilon_{i}\sim\mathcal{N}\left(0,\sigma_{n}^{2}\right)$. We denote by $\bm{y}=[y_{1},\ldots,y_{N}]^{T}$ the vector of observations and by $\bm{X}=[\bm{x}_{1}^{T},\ldots,\bm{x}_{N}^{T}]^{T}\in\mathbb{R}^{N\times D}$ the input points.

We model $f$ with a Gaussian Process (GP), a stochastic process defined by its mean function $m(\bm{x})$ and covariance kernel $k(\bm{x},\bm{x}^{\prime})$. The kernel $k$ is a positive definite function [28], such as, for instance, the squared exponential (SE) kernel with individual lengthscales $l_{i}$ for each dimension, that is $k(\bm{x},\bm{x}^{\prime})=\sigma_{0}^{2}\exp\left(-\frac{1}{2}\left(\bm{x}-\bm{x}^{\prime}\right)^{T}\text{Diag}\left[l_{1}^{2},\ldots,l_{D}^{2}\right]\left(\bm{x}-\bm{x}^{\prime}\right)\right)$. We assume $m(x)\equiv 0$ for the sake of simplicity and we use the SE kernel throughout this paper however all methods work with any positive definite kernel. Given the training values $\bm{f}=f\left(\bm{X}\right)=\left[f(\bm{x}_{1}),\ldots,f(\bm{x}_{N})\right]^{T}$ and a test latent function value $f_{*}=f(\bm{x}_{*})$ at a test point $\bm{x}_{*}\in\mathbb{R}^{D}$, then the joint distribution $p(\bm{f},f_{*})$ is Gaussian. Our likelihood is Gaussian, $p\left(\bm{y}|\bm{f}\right)=\mathcal{N}\left(\bm{y}|\bm{f},\sigma_{n}^{2}\mathbb{I}\right)$, and with Bayes theorem (see e.g. [28]) we obtain analytically the posterior predictive distribution $p\left(f_{*}|\bm{y}\right)=\mathcal{N}\left(f_{*}|\bm{\mu}_{*},\bm{\Sigma}_{*}\right)$ with

[TABLE]

where $\left[\bm{K}_{\bm{A}\bm{B}}\right]_{ij}=k(\bm{a}_{i},\bm{b}_{j})$ for any $\bm{A}\in\mathbb{R}^{{M_{1}}\times D}$ and $\bm{B}\in\mathbb{R}^{{M_{2}}\times D}$ with the corresponding rows $\bm{a}_{i},\bm{b}_{j}$. For brevity, we use $*$ to indicate $\bm{x}_{*}$. The GP depends via the kernel matrices on the hyper-parameters $\bm{\phi}=\{\sigma_{0},l_{1},\ldots,l_{D},\sigma_{n}\}$ typically estimated by maximizing the log marginal likelihood

[TABLE]

Note that the computations for inference require the inversion of the matrix in Eq. (1) which scales as $\mathcal{O}(N^{3})$ in time and $\mathcal{O}(N^{2})$ for memory (given $\bm{\phi}$).

2.1 Batch Sparse GP Regression

Sparse GP regression methods based on inducing points approaches reduce the computational complexity by introducing $M\ll N$ inducing points $\bm{u}\in\mathbb{R}^{M}$ that optimally summarize the dependency of the whole training data. The inducing inputs $\bm{R}\in\mathbb{R}^{M\times D}$ are in the $D$-dimensional input data space and the inducing outputs $\bm{u}:=f(\bm{R})$ are the corresponding GP-function values, see also Fig. 1. The GP prior over $\bm{f}$ and $f_{*}$ is augmented with the inducing outputs $\bm{u}$, leading to a joint $p\left(\bm{f},f_{*},\bm{u}\right)$ and marginal $p\left(\bm{u}\right)=\mathcal{N}\left(\bm{0},\bm{K}_{\bm{R}\bm{R}}\right)$ prior. By marginalizing out the inducing points, the original prior $p\left(\bm{f},f_{*}\right)=\int p\left(\bm{f},f_{*}|\bm{u}\right)p\left(\bm{u}\right)\mathop{}\!\mathrm{d}\bm{u}$ is recovered. The fundamental approximation in all sparse GP models is that given the inducing outputs $\bm{u}$, $\bm{f}$ and $f$ are conditionally independent. Consequently, inference in these models can be done in $\mathcal{O}(M^{2}N)$ time and $\mathcal{O}(MN)$ space [34].

We briefly recall here the sparse predictive distribution and the variational lower bound to the log marginal likelihood for the Power Expectation Propagation (PEP) model [7] because it unifies the main sparse inducing points approaches. The variational lower bound is used for optimizing the parameters $\bm{\theta}:=\{\bm{\phi},\bm{R}\}$. In the following, we denote $\bm{Q}_{\bm{A}\bm{B}}=\bm{K}_{\bm{A}\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}\bm{K}_{\bm{R}\bm{B}}$ and $\bm{D}_{\bm{A}}=\bm{K}_{\bm{A}\bm{A}}-\bm{Q}_{\bm{A}\bm{A}}$ for any $\bm{A},\bm{B}$. The predictive distribution $p\left(f_{*}|\bm{y}\right)=\mathcal{N}\left(f_{*}|\bm{\mu}_{*},\bm{\Sigma}_{*}\right)$ of PEP is given by

[TABLE]

where $\overline{\bm{K}}_{\bm{X}\bm{X}}=\bm{Q}_{\bm{X}\bm{X}}+\alpha\text{Diag}\left[\bm{D}_{\bm{X}}\right]$. A lower bound to the sparse log marginal likelihood is analytically available

[TABLE]

where we omit the explicit dependency on $\bm{\theta}$ via $\overline{\bm{K}}_{\bm{X}\bm{X}}$ and $\bm{D}_{\bm{X}}$ for the sake of brevity. This bound can be used to learn the parameters $\bm{\theta}$, similarly to Eq. (2) for full GP.

The special case $\alpha\rightarrow 0$ was originally introduced in [36] where the author proposed to maximize a variational lower bound to the true GP marginal likelihood, obtaining the Variational Free Energy (VFE) or the collapsed lower bound

[TABLE]

In (5) the variational distribution over the inducing points is optimally eliminated and analytically available. The rightmost term in (5) acts as a regularizer that prevents overfitting and has the effect that the sparse GP predictive distribution (3) converges [36] to the exact GP predictive distribution (1) as the number of inducing points increases, when optimizing $\bm{\theta}$ with (5). See also [7, 20, 27, 28] for recent reviews on the subject.

2.2 Sequential Sparse GP Regression

The optimization for $\bm{\theta}$ of the collapsed lower bound (5) requires to process the whole dataset, which is very inefficient and not feasible for large $N$. We would like to update the parameters more frequently, therefore, we split the data $\mathcal{D}=\left\{\bm{y}_{k},\bm{X}_{k}\right\}_{k=1}^{K}$ into $K$ mini-batches of size $B$ and denote $\bm{f}_{k}$ the corresponding sparse GP value. Stochastic Variational Gaussian Process (SVGP) [15] achieves this result by applying stochastic optimization to an uncollapsed lower bound to the log marginal likelihood

[TABLE]

where the variational distribution $q(\bm{u})$ is part of the bound and explicitly parametrized as $\bm{q}(\bm{u})=\mathcal{N}\left(\bm{u}|\bm{\mu},\bm{\Sigma}\right)$. This uncollapsed bound satisfies $\mathcal{L}_{SVGP}(\bm{\mu},\bm{\Sigma},\bm{\theta})\leq\mathcal{L}_{VFE}(\bm{\theta})$ with equality when inserting the optimal mean and covariance of the variational distribution of VFE. The key property of this bound is that it can be written as a sum of $K$ terms, which allows Stochastic Variational Inference (SVI, [17]). Note that collapsing the bound, i.e. inserting the optimal distribution, reintroduces dependencies between the observations, and eliminates the global parameter $\bm{u}$ which is needed for SVI. For this reason, all variational parameters are numerically estimated by following the noisy gradients of a stochastic estimate of the lower bound $\mathcal{L}_{SVGP}$. By passing through the training data a sufficient number of times, the variational distribution converges to the batch solution of VFE method. This approach, however, requires a large number of parameters: in addition to the parameters $\bm{\theta}$, all entries in the mean vector $\bm{\mu}$ and the covariance matrix $\bm{\Sigma}$ have to be estimated numerically, which is in order $\mathcal{O}(M^{2})$.

3 Recursive Sparse GP Regression

In this section we establish the connection between Bayesian recursive estimation and sparse inducing point GP models. We recall the weight-space view for a large class of sparse inducing point GP models, which we present here as a particular kernelized version of a Bayesian linear regression model. See [28, Ch. 2.1], for an analogous discussion on the full GP model. Here, however, we show how to exploit the KF to train many sparse methods analytically either in an online setting for fixed hyper-parameters. This allows us to introduce a recursive log marginal likelihood with a model specific regularization term for parameter estimation.

3.1 Weight-Space View of Generic Sparse GP

For a mini-batch $\bm{X}\in\mathbb{R}^{B\times D}$ of size $B$, consider the generic sparse GP model

[TABLE]

where the sparse GP value $f\left(\bm{X}\right)$ is modeled by a linear combination of basis-functions $H\left(\bm{X}\right)\in\mathbb{R}^{B\times M}$, (stochastic) weights $\bm{u}\in\mathbb{R}^{M}$ with a prior $p\left(\bm{u}\right)=\mathcal{N}\left(\bm{0},\bm{\Sigma}_{0}\right)$ and an input dependent error term $\ \gamma\left(\bm{X}\right)\sim\mathcal{N}\left(0,V\left(\bm{X}\right)\right)$ that takes into account the sparse approximation. For $k=1,\ldots,K$, the noisy observations $\bm{y}_{k}$ are obtained by adding independent noise $\bm{\varepsilon}_{k}\sim\mathcal{N}\left(0,\sigma_{n}^{2}\mathbb{I}\right)$ to $f\left(\bm{X}_{k}\right)$, yielding the model

[TABLE]

where we distinguish the training $\bm{f}_{k}=f\left(\bm{X}_{k}\right)$ and test $\bm{f}_{*}=f\left(\bm{X}_{*}\right)$ cases depending on the input $\bm{X}_{k}$ and $\bm{X}_{*}$. Assuming $\bm{\gamma}_{k}$, $\bm{\gamma}_{*}$ and $\bm{\varepsilon}_{k}$ are independent, by linearity and Gaussianity we can compactly write

[TABLE]

Combining (10) and (12) and by integrating out $\bm{f}_{k}$ we obtain the likelihood $p\left(\bm{y}_{k}|\bm{u}\right)=\mathcal{N}\left(\bm{y}_{k}|\bm{H}_{k}\bm{u},\bm{V}_{k}\right)$ where $\bm{V}_{k}=\overline{\bm{V}}_{k}+\sigma_{n}^{2}\mathbb{I}$. This shows that a generic sparse GP regression model can be seen as a Bayesian non-linear regression model with additional input dependent observation noise, a particular choice of basis functions $\bm{H}_{k}$ and covariance structures $\overline{\bm{V}}_{k}$, $\bm{V}_{*}$. For inducing inputs $\bm{R}\in\mathbb{R}^{M\times D}$, we have $\bm{\Sigma}_{0}=\bm{K}_{\bm{R}\bm{R}}$, $\bm{H}_{k}=\bm{K}_{\bm{X}_{k}\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$ and $\bm{H}_{*}=\bm{K}_{*\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$. Different choices of the quantities $\overline{\bm{V}}_{k}$ and $\bm{V}_{*}$ lead to a range of sparse GP models summarized in the bottom table in Fig. 2.

3.2 Training

Given the prior $p\left(\bm{u}\right)$ and the likelihood $p\left(\bm{y}_{k}|\bm{u}\right)$, the posterior over the weights $\bm{u}$ conditioned on the data $\bm{y}_{1:k}$ can be computed either in a batch or in a recursive manner.

3.2.1 Batch Estimation

The batch likelihood is $p\left(\bm{y}|\bm{u}\right)=\prod_{k=1}^{K}p\left(\bm{y}_{k}|\bm{u}\right)=\mathcal{N}\left(\bm{y}|\bm{H}\bm{u},\bm{V}\right)$ with $\bm{H}=\left[\bm{H}_{1}^{T},\ldots,\bm{H}_{K}^{T}\right]^{T}\in\mathbb{R}^{N\times M}$ and $\bm{V}$ a block-diagonal matrix with blocks $\bm{V}_{k}$. The posterior over $\bm{u}$ given the data $\bm{y}$ can be obtained by Bayes’ rule, i.e.

[TABLE]

with $\bm{\Sigma}_{K}=\left(\bm{\Sigma}_{0}^{-1}+\bm{H}^{T}\bm{V}^{-1}\bm{H}\right)^{-1}$ and $\bm{\mu}_{K}=\bm{\Sigma}_{K}\bm{H}^{T}\bm{V}^{-1}\bm{y}.$

3.2.2 Recursive Estimation

An equivalent solution can be obtained by propagating recursively $p\left(\bm{u}|\bm{y}_{1:k-1}\right)$. By interpreting this previous posterior as the prior, the updated posterior can be recursively computed by

[TABLE]

where $\bm{\mu}_{k}=\bm{\Sigma}_{k}\left(\bm{H}_{k}^{T}\bm{V}_{k}^{-1}\bm{y}_{k}+\bm{\Sigma}_{k-1}^{-1}\bm{\mu}_{k-1}\right)$ and $\bm{\Sigma}_{k}=\left(\bm{\Sigma}_{k-1}^{-1}+\bm{H}_{k}^{T}\bm{V}_{k}^{-1}\bm{H}_{k}\right)^{-1}$.

Kalman Filter like updating

the KF constitutes an efficient way to update the mean and covariance of $p\left(\bm{u}|\bm{y}_{1:k}\right)$. Applying the Woodbury identity to $\bm{\Sigma}_{k}$ in Eq. (15) and introducing temporary variables yields

[TABLE]

Starting the recursion with $\bm{\mu}_{0}=\bm{0}$ and $\bm{\Sigma}_{0}$, the posterior distribution at step $K$ is equivalent to (14) independent of the order of the data. We want to emphasize that the only difference in the estimation part between the sparse GP models is the form of the additional noise $\bm{V}_{k}=\overline{\bm{V}}_{k}+\sigma_{n}^{2}\mathbb{I}$.

Transformation

instead of running a KF with $\bm{\Sigma}_{0}=\bm{K}_{\bm{R}\bm{R}}$, $\bm{H}_{k}=\bm{K}_{\bm{X}_{k}\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$ and $\bm{H}_{*}=\bm{K}_{*\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$, an equivalent predictive distribution is also obtained when using $\tilde{\bm{\Sigma}}_{0}=\bm{K}_{\bm{R}\bm{R}}^{-1}$ and $\tilde{\bm{H}}_{k}=\bm{K}_{\bm{X}_{k}\bm{R}}$ together with $\tilde{\bm{H}}_{*}=\bm{K}_{*\bm{R}}$. For any $k$ we then propagate a transformed posterior distribution $\tilde{\bm{\mu}}_{k}=\bm{K}_{\bm{R}\bm{R}}^{-1}\bm{\mu}_{k}$, $\tilde{\bm{\Sigma}}_{k}=\bm{K}_{\bm{R}\bm{R}}^{-1}\bm{\Sigma}_{k}\bm{K}_{\bm{R}\bm{R}}^{-1}$ and $\bm{\mu}_{k}=\bm{K}_{\bm{R}\bm{R}}\tilde{\bm{\mu}}_{k}$, $\bm{\Sigma}_{k}=\bm{K}_{\bm{R}\bm{R}}\tilde{\bm{\Sigma}}_{k}\bm{K}_{\bm{R}\bm{R}}$, respectively. This parametrization constitutes a computational shortcut, since the basis functions are very easy to interpret and do not include any matrix multiplication. Note that also the log marginal likelihood discussed below is not affected by this transformation.

3.3 Prediction

Given a new $\bm{X}_{*}\in\mathbb{R}^{A\times D}$, the predictive distribution after seeing $\bm{y}_{1:k}$ of the sparse GP methods can be obtained by

[TABLE]

using $p\left(\bm{f}_{*}|\bm{u}\right)=\mathcal{N}\left(\bm{f}_{*}|\bm{H}_{*}\bm{u},\bm{V}_{*}\right)$ with $\bm{H}_{*}=\bm{K}_{*\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$ and $\bm{V}_{*}$ the model specific prediction covariance. The predictions for $\bm{y}^{*}$ are obtained by adding $\sigma_{n}^{2}\mathbb{I}$ to the covariance of $\bm{f}_{*}|\bm{y}_{1:k}$. At step $K$, by applying the Woodbury identity to the batch covariance $\bm{\Sigma}_{K}$ in (14), we get for the predictive distribution in (17)

[TABLE]

where $\bm{\Sigma}=\left(\bm{H}\bm{\Sigma}_{0}\bm{H}^{T}+\bm{V}\right)^{-1}$. Inserting the particular choices for $\bm{\Sigma}_{0}$, $\bm{H}$ and $\bm{H}_{*}$ yields the usual formulation for the sparse predictive distribution

[TABLE]

Depending on the choice of the covariances $\overline{\bm{V}}$ and $\bm{V}_{*}$, we obtain for instance (3) for PEP, or the analogous predictions for VFE and FITC, respectively.

3.4 Online learning

This connection between sparse GP models and recursive estimation allows us to train the sparse GP models analytically online for streaming data for fixed $\bm{\theta}$.

As an illustrative example consider $N=100$ data samples in $D=1$, we are interested in training a PEP model with $\alpha=0.5$ with $M=15$ inducing points with fixed $\bm{\theta}$. Let’s assume that the data samples $\{\bm{x}_{k},y_{k}\}$ arrive sequentially in a stream. Thus we have $B=1$ and $K=100$. Here we use the transformation and we apply, for each data sample $k$, the recursion in (16). We note that here $\tilde{\bm{r}}_{k}$, $\tilde{\bm{V}}_{k}$ and $\tilde{\bm{S}}_{k}$ are numbers as $B=1$ and $\tilde{\bm{H}}_{k},\tilde{\bm{G}}_{k}\in\mathbb{R}^{15}$.

We obtain predictions for new data $\bm{X}_{*}$, by applying (17) with $\tilde{\bm{H}}_{*}=\bm{K}_{*\bm{R}}$ and $\bm{V}_{*}=\bm{K}_{*\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}\bm{K}_{\bm{R}*}$. Note that there is no need to transform back the posterior over the inducing points, since it is already taken into account in the prediction step. After processing all $N$ samples, the predictive distribution and the cumulative bound of log marginal likelihood correspond to the batch version, as shown in Fig. 3.

3.5 Marginal Likelihood

In the batch setting, the log marginal likelihood $\log p\left(\bm{y}\right)$ can be computed by marginalizing out $\bm{u}$, that is

[TABLE]

In the recursive setting, $p\left(\bm{y}\right)$ can be factorized into $\prod_{k=1}^{K}p\left(\bm{y}_{k}|\bm{y}_{1:k-1}\right),$ where

[TABLE]

The log of the joint marginal likelihood involving all terms of (3.5) can be explicitly written as

[TABLE]

The iterative maximization of a lower bound of the recursive factorized marginal likelihood in (3.5) leads to the recursive KF updates in (16) for the posterior and to the lower bound

[TABLE]

which includes a model specific regularization term $a_{k}$ (see the right-most column in Fig. 2). We refer to this as the recursive collapsed bound and a detailed derivation for the VFE model is given in App. B. Using the model specific quantities $\overline{\bm{V}}_{k}$ and $a_{k}$, this recursive computation of the lower bound of the marginal likelihood are equivalent to the batch counterparts for all sparse models, for instance (5) and (4) for VFE and PEP, respectively.

4 Hyper-parameters Estimation

The previous section presented an online procedure for training sparse GP models at fixed hyper-parameters $\bm{\theta}$. Here we show that, by exploiting the connections highlighted before, we can optimize $\bm{\theta}$ sequentially. The recursive collapsed bound in (23) decomposes into a recursive sum over the mini-batches which allows to optimize the hyper-parameters $\bm{\theta}$ sequentially as opposed to the collapsed bound (5). Our bound (23) enables the application of stochastic optimization without needing to estimate all entries in the posterior mean vector and covariance matrix as in SVGP. Compared to the uncollapsed bound in (LABEL:eq:uncollapsed_bound), the variational distribution is recursively and analytically eliminated, thus reducing the number of parameters to be numerically estimated drastically from $\mathcal{O}(MD+M^{2})$ to $\mathcal{O}(MD)$.

Finding a maximizer $\bm{\theta}\in\bm{\Theta}$ of an objective function $\Psi(\bm{\theta})=\sum_{k=1}^{K}\psi_{k}(\bm{\theta})$ can be achieved by applying Stochastic Gradient Descent(SGD), with the update

[TABLE]

where $\gamma^{(t-1)}$ might be a sophisticated function of $\bm{\theta}^{(0)},\ldots,\bm{\theta}^{(t-1)}$ (for instance using ADAM [18], where also a bias correction term is included). We call one pass over the $K$ mini-batches an epoch. We denote $\bm{\theta}^{(e,k)}\in\bm{\Theta}^{(e,k)}$ the estimate of $\bm{\theta}$ in epoch $e\in E$ for mini-batch $k$.

4.1 Recursive Gradient Propagation (RGP)

We rewrite the recursive collapsed bound in (23) as

[TABLE]

where $d_{k}(\bm{\theta})=\log\mathcal{N}\left(\bm{r}_{k}^{\bm{\theta}}|\bm{0},\bm{S}_{k}^{\bm{\theta}}\right)$ and $a_{k}(\bm{\theta})$ the model specific regularization term. Since $\psi^{(K)}(\bm{\theta})$ decomposes into a (recursive) sum over the mini-batches, we directly compute the derivative of $\psi_{k}(\bm{\theta})$ w.r.t. $\bm{\theta}\in\bm{\Theta}$. The derivative of $a_{k}$ is straightforward, for $d_{k}$ we have

[TABLE]

with $\bm{r}_{k}^{\bm{\theta}}=\bm{y}_{k}-\bm{H}_{k}^{\bm{\theta}}\bm{\mu}_{k-1}^{\bm{\theta}}$ and $\bm{S}_{k}^{\bm{\theta}}=\bm{H}_{k}^{\bm{\theta}}\bm{\Sigma}_{k-1}^{\bm{\theta}}(\bm{H}_{k}^{\bm{\theta}})^{T}+\bm{V}_{k}^{\bm{\theta}}$. It is important to note that ignoring naively the dependency of $\bm{\theta}$ through $\bm{\mu}_{k-1}^{\bm{\theta}}$ and $\bm{\Sigma}_{k-1}^{\bm{\theta}}$ completely forgets the past and thus results in overfitting the current mini-batch. In order to compute the derivatives of $\bm{\mu}_{k}^{\bm{\theta}}$ and $\bm{\Sigma}_{k}^{\bm{\theta}}$, we exploit the chain rule for derivatives and recursively propagate the gradients of the mean and the covariance over time, that is

[TABLE]

where $\frac{\partial\bm{G}_{k}}{\partial\theta}$, $\frac{\partial\bm{r}_{k}}{\partial\theta}$ and $\frac{\partial\bm{S}_{k}}{\partial\theta}$ are computed recursively according to (16). Computing the derivatives of $d_{k}$ as explained in (26) and (27), the stochastic gradient

[TABLE]

can be computed for each mini-batch $k$.

Proposition 1**.**

*Consider a fixed parameter vector $\bar{\bm{\theta}}\in\bm{\Theta}$, the gradients of the full batch lower bound $\mathcal{L}_{PEP}(\bar{\bm{\theta}})$ in (4) with respect to $\bm{\theta}$ are equal to the cumulative partial derivatives $\frac{\partial\psi^{(K)}}{\partial\bm{\theta}}$ of the recursive collapsed bound with the above learning procedure. That is, it holds $\frac{\partial\mathcal{L}_{PEP}(\bar{\bm{\theta}})}{\partial\bm{\theta}}=\frac{\partial\psi^{(K)}(\bar{\bm{\theta}})}{\partial\bm{\theta}}.$ *

This shows that, when the gradients are cumulated over all data samples, each gradient step of our recursive procedure is equivalent to a gradient step for $\mathcal{L}_{PEP}$ and, therefore, follows the gradients of an optimal collapsed lower bound.

The toy example in Fig. 3 also shows this equivalence. The numbers in the bottom left and right corners show the cumulative recursive collapsed bound $\psi^{(k)}$ and its cumulative derivative $\frac{\partial\psi^{(k)}}{\partial\l}$ (abbreviated as $\dot{\psi}^{(k)}$) with respect to the lengthscale. The lower bound of the marginal likelihood as well as its derivatives are exactly the same value as the corresponding batch counterpart in Fig. 1.

For Stochastic Recursive Gradient Propagation (SRGP), in each epoch $e$ and mini-batch $k$, we interleave the update step of the inducing points in Eq. (15) with the SGD update (24) of the parameters $\bm{\theta}^{(e,k)}$, i.e.

[TABLE]

More concretely, we update after each mini-batch $k$ the parameters $\bm{\theta}^{(e,k)}$ with (28),(24) and propagate recursively the posterior with (16) and its derivative (27). In order to compute all the derivatives with respect to $\bm{\theta}\in\bm{\Theta}$, we exploit several matrix derivative rules which simplify the computation significantly, see App. A. Finally note that the form of the gradients in eq. (26) and (28) implies that the noise in the stochastic gradient at step $k$ depends on the noise at step $k$-1, thus excluding standard convergence proofs such as [5]. Recent results on non-convex optimization problems [10, 40, 41] show convergence proofs for function classes that include our objective, nonetheless the Markov noise of our problem still excludes it from this general theory. For this reason we leave a convergence analysis to future work.

4.2 Computational complexity

In the following, we assume that the batch size $B$ is larger than the number of inducing points $M$. For one mini-batch, the time complexity to update the posterior is dominated by matrix multiplications of size $B$ and $M$, thus $\mathcal{O}(B^{2}M)$. In order to propagate the gradients of the posterior and to compute the derivative of the bound needs $\mathcal{O}(BM^{2})$ for a mini-batch and a parameter $\bm{\theta}\in\bm{\Theta}$. Thus, updating a mini-batch including all $\mathcal{O}(MD)$ parameters costs $\mathcal{O}(BM^{3}D+B^{2}M)$ for the SRGP method. Since SRGP stores the gradients of the posterior, it requires $\mathcal{O}(M^{3}D+BM)$ storage.

On the other hand, SVGP needs $\mathcal{O}(M^{2}+BM)$ storage and $\mathcal{O}(BM^{3}+B^{2}M)$ time per mini-batch, where the latter can be broken down into once $\mathcal{O}(B^{2}M)$ and $\mathcal{O}(BM)$ for each of the $\mathcal{O}(M^{2})$ parameters. This means, for moderate dimensions, our algorithm has the same time complexity as state of the art method SVGP. However, due to the analytic updates of the posterior we achieve an higher accuracy and less epochs are needed as shown in Fig. 4 and in Sect. 5 empirically.

Fig. 4 shows the convergence of SRGP on a 1-D toy example with $N=1000$ data samples and $M=15$ inducing points. The parameters are sequentially optimized with our recursive approach (blue) and as comparison with SVGP (green) with a mini-batch size of $B=100$ over several epochs. The root-mean-squared-error (RMSE) computed on test points, the bound of the log marginal likelihood (LML) as well as the hyper-parameters converge in a few iterations to the corresponding batch values of VFE (red). Due to the analytic updates of the posterior, the accuracy is higher and SRGP needs much less epochs until convergence.

4.3 Mini-batch size

The size of the mini-batches has an impact on the speed of convergence of the algorithm. Proposition 1 tells us that if we use a full batch, our algorithm requires the same number of gradient updates as a full batch method to converge. On the other hand smaller batches should require more updates and should lead to a higher variance in the results. Fig. 5 shows a comparison of different mini-batch sizes on a 1-D toy example with $N=10^{\prime}000$ data samples generated with the same parameters as in Sect. 5.1. The convergence to the full batch value is slower as the batch size decreases. Moreover the variance of the error, over the repetitions, is much larger for smaller batch sizes: in the last 10 normalized gradient updates, the standard deviation of the error is on average $3.8\times 10^{-3}$ for $B=100$ and $8.1\times 10^{-4}$ for $B=5^{\prime}000$, denoting a more stable procedure for higher batch sizes. As the mini-batch size increases the computational cost for each gradient update also increases. In this example one gradient update requires on average $4.9\times 10^{-3}$ sec and $5.8\times 10^{-2}$ sec with $B=100$ and $B=5^{\prime}000$ respectively.222The times are measured on a laptop with a Intel i5-7300U CPU $@$ 2.6 GHz. These considerations suggest that a reasonable choice is a large mini-batch size within the computational and time budgets.

5 Experiments

We first benchmark our method with $N=100^{\prime}000$ synthetic data samples generated by a GP in several dimensions. Next, we apply our approach to the Airline data used in [15] with a million of data samples. Finally a more realistic setup is presented where we use up to a million data samples to train a nonlinear plant. We compare our SRGP method to full GP and sparse batch method VFE for a subset of data (using the implementation in GPy [13]) and to the state of the art stochastic parameter estimation method SVGP implemented in GPflow [22]. Our algorithm works also for many other sparse models, however, only large-scale implementations of standard SVGP are available (corresponding to the VFE model), thus we restrict the investigation to this model.

5.1 GP Simulation

In this section we test our proposed learning procedure on simulated GP data. We generate $N=100^{\prime}000$ data samples from a zero-mean (sparse) GP with SE covariance kernel with hyper-parameters $\sigma_{0}=1,\sigma_{n}=0.1$ and $l=\{0.1,0.2,0.5\}$ in $D=\{1,2,5\}$ dimensions. The initial $M=\{20,50,100\}$ inducing points are randomly selected points from the data and the hyper-parameters of a SE kernel with individual lengthscales for each dimension are initialized to the same values for both algorithms ($\sigma_{0}=1,\sigma_{n}=1,l_{1},\ldots,l_{D}=1$). All parameters are sequentially optimized with our recursive approach and with SVGP with a mini-batch size of $B=5000$. The stochastic gradient descent method ADAM [18] is employed for both methods with learning rates $\{0.001,0.005,0.005\}$ for SVGP and $\{0.0001,0.001,0.005\}$ for SRGP (based on some preliminary experiments). Each experiment is replicated $10$ times.

Fig. 6 shows the bound to the log marginal likelihood, the RMSE and the coverage of $10^{\prime}000$ test points for the data dimensions $D=\{1,2,5\}$ of both methods over $50$ epochs. The shaded lines indicates the $10$ repetitions and the thick line correspond to the mean. The recursive propagation of the gradients achieves faster convergence and more accurate performance regarding mean RMSE and smaller values for the log marginal likelihood. The higher accuracy and faster convergence can possibly be explained by the analytic updates of the posterior mean and covariance which leads to less parameters to be optimized numerically.

5.2 Airline Data

For the second example we apply our recursive method to the Airline Data used in [15]. It consists of flight arrival and departure times for more than 2 millions flights in the USA from January 2008 and April 2008. We preprocessed the data as similar as possible as described in [15] resulting in 8 variables: age of the aircraft, distance that needs to be covered, airtime, departure time, arrival time, day of the week, day of the month and month. We trained our recursive method as well as SVGP with an SE kernel on $N=1^{\prime}000^{\prime}000$ data samples with $M=500$ inducing points randomly selected from the data and a mini-batch size of $B=10^{\prime}000$. The ADAM learning rates are set to $0.005$ for both methods and the size of the test set is $50^{\prime}000$. For $5$ different repetitions, the RMSE as a function of epochs is depicted in Fig. 7. The mean coverage on test data (at $95\%$) is comparable for both methods with values of $0.92$ and $0.97$ for SVGP and SRGP respectively. The overall performance of SRGP is superior to SVGP.

5.3 Non-Linear Plant

GPs are a powerful way to model complex functions in a non-parametric way, thus they are suitable to learn the complex input output behavior of a non-linear plant. However, with full or even sparse batch GP methods the use is restricted to a few thousands of samples. With our sequential learning method, we are able to exploit the huge amount of available data by training with up to a million of samples.

We consider a Continuous Stirred Tank Reactor (CSTR). The dynamic model of the plant is

[TABLE]

where $C_{b}(t)$ is the product concentration at the output of the process, $h(t)$ is the liquid level, $w_{1}(t)$ is the flow rate of concentrated feed $C_{b1}$, and $w_{2}(t)$ is the flow rate of the diluted feed $C_{b2}$. The input concentrations are $C_{b1}=24.9$ and $C_{b2}=0.1$. The constants associated with the rate of consumption are $k_{1}=k_{2}=1$. The objective of the controller is to maintain the product concentration by changing the flow $w_{1}(t)$. To simplify the example, we assume that $w_{2}(t)=0.1$ and that the level of the tank $h(t)$ is not controlled. We denote the controlled outputs $C_{b}(t),C_{b}(t-1),\dots,C_{b}(t-p)$ as $f_{t},f_{t-1},\ldots,f_{t-p}$ and the control variables as $w_{t},w_{t-1},\ldots,w_{t-p}$. Therefore, the plant identification problem can be shaped into the problem of estimating the non-linear function $f_{t}=g(f_{t-1},\ldots,f_{t-p},w_{t},w_{t-1},\ldots,w_{t-p})$ which depends on the $p$ previous values as well as on the current and the $p$ past control values $w$.

However, we can only observe a noisy version of the controlled response, that is $y_{t}=f_{t}+\varepsilon_{t}$ with $\varepsilon\sim\mathcal{N}\left(0,\sigma_{n}^{2}\right)$. Using a sampling rate of $0.2s$, we have generated $1^{\prime}200^{\prime}000$ observations (about $3$ days of observations). The plant input is a series of steps, with random height (in the interval $[0,4]$), occurring at random intervals (in the interval $[5,20]s$). For different numbers $N_{train}$, we use the samples $y_{10^{6}-N_{train}},\ldots,y_{10^{6}}$ for training and the last $200^{\prime}000$ are used as a test set. The goal is to learn a model for the controlled response $y_{t}$ given $\bm{x}_{t}=[y_{t-1},y_{t-2},w_{t},w_{t-1},w_{t-2}]^{T}\in\mathbb{R}^{5}$ for the particular choice of $p=2$. We model the non-linear function $g$ with a GP with a SE kernel. For comparison, we train full GP and sparse batch GP (with $100$ inducing points) on a time horizon $N_{train}$ of up to $10^{\prime}000$ and $50^{\prime}000$ past values, respectively. With the sequential version SVGP and our recursive gradient propagation method SRGP (both with $100$ inducing points and mini-batch size of $1^{\prime}000$), we use a time horizon of up to a million. This situation is depicted in Fig. 8, where for $1500$ training samples $y_{t}$ (red dots), the true (unknown) function $f_{t}$ (green) and the control input $w_{t}$ (grey) is shown together with the predicted values with full GP (red dotted) and recursive GP (blue dotted) trained on a time horizon of $1^{\prime}000$ and $10^{\prime}000$, respectively. In Fig. 9, the RMSE and the median computed on the test set (with $10$ repetitions) is depicted for full GP, sparse GP (VFE), SVGP and SRGP trained with varying time horizons. For small and medium training sizes, when the batch methods are applicable, our recursive method achieves the same performance as the batch counterpart (VFE) and is comparable to full GP. Due to the analytic updates of the posterior, SRGP outperforms SVGP regarding both RMSE and median for all training sizes. By exploiting more than several thousand past values, a significant increase in performance of SRGP can be still observed, thus it constitutes an approach to accurately scale GPs up to a million of past values.

6 Conclusion

In this paper we introduced a recursive inference and parameter estimation method SRGP for a general class of sparse GP approximations. Since the posterior updates are given analytically, one pass through the data is sufficient to compute the posterior for given parameters. For parameter estimation, we proposed a recursive collapsed bound to the log marginal likelihood that matches exactly the batch version but can be used for stochastic estimation. Due to the analytic updates of the posterior our method has much less parameters to be estimated numerically. As a consequence, the experimental section showed that our recursive method needs less epochs and has superior accuracy compared to state of the art, thus constitutes an efficient methodology for scaling GPs to big data problems.

Our approach could be enhanced in several directions. While the proposed method only exploits the update equations of the KF, an interesting direction would be to include a dynamic in a state space model that takes into account the varying hyper-parameters which makes it also applicable for the streaming setting as [6]. Moreover, we further plan to investigate distributed parameter estimation based on an information filter formulation of the problem. Finally, we aim to provide a convergence analysis of the proposed method to recursively learn the hyper-parameters using a similar approach recently employed for ADAM or AMSGrad [29, 38].

{ack}

This work is supported by the Swiss National Research Programme 75 ”Big Data” grant n. 407540_167199 / 1.

Appendix A Details for Recursive Gradient Propagation

We show here the computation for the recursive gradient propagation from Sect. 4.1 for the PEP model. For other models, $a_{k}$ and $\overline{\bm{V}}$ from the Table 2 could be used correspondingly. We will use the following notation: $\text{diag}\left[\bm{A}\right]=\bm{d},d_{i}=a_{ii}$, $\text{Diag}\left[\bm{d}\right]=\bm{A},a_{ii}=d_{i},a_{ij}=0$, $\bm{A}\odot\bm{B}=\bm{C},c_{ij}=a_{ij}b_{ij}$, $\bm{A}\div\bm{B}=\bm{C},c_{ij}=\frac{a_{ij}}{b_{ij}}$, $\bm{A}^{\odot 2}=\bm{C},c_{ij}=a_{ij}^{2}$, $\dot{\bm{A}}=\frac{\partial\bm{A}(\theta)}{\partial\theta},\forall\theta\in\bm{\Theta}$, $\text{sum}\left[\bm{A}\right]=\sum_{i,j}a_{ij}$, $\mathrm{1}_{[z]}=1\text{ if }z=\text{true},0\text{ other.}$

Initialization

[TABLE]

Natural Mean and Precision Updates

[TABLE]

Intermediate Derivatives

[TABLE]

**Derivative Updates

Loop over $\theta_{i}\in\bm{\theta}$:**

[TABLE]

[TABLE]

For the noise $\sigma_{n}$, all kernel derivatives are zero, therefore the calculations simplify significantly.

Appendix B Derivation of Recursive Collapsed Bound

We provide more details and a detailed derivation for the recursive collapsed lower bound (25) for the VFE model from Sect. 4.1. Instead of lower bounding directly the batch log marginal likelihood as done by Titsias [36] in the batch case, our approach relies on the recursive factorization of the joint log marginal likelihood $\log p\left(\bm{y}|\bm{\theta}\right)=\log\prod_{k=1}^{K}p\left(\bm{y}_{k}|\bm{y}_{1:k-1},\bm{\theta}\right)=\sum_{k=1}^{K}\log p\left(\bm{y}_{k}|\bm{y}_{1:k-1},\bm{\theta}\right).$ The properties induced by the sparse augmented inducing point model yield $\log p\left(\bm{y}|\bm{\theta}\right)=$

$\sum_{k=1}^{K}\log\int p\left(\bm{y}_{k}|\bm{f}_{k},\bm{\theta}\right)p\left(\bm{f}_{k}|\bm{u},\bm{\theta}\right)p\left(\bm{u}|\bm{y}_{1:k-1},\bm{\theta}\right)\mathop{}\!\mathrm{d}{\bm{f}_{k}}\mathop{}\!\mathrm{d}\bm{u}.$ We first introduce the variational distributions $q_{k}(\bm{f}_{k},\bm{u})=p\left(\bm{f}_{k}|\bm{u},\bm{\theta}\right)q_{k}(\bm{u})\approx p\left(\bm{f}_{k},\bm{u}|\bm{y}_{1:k},\bm{\theta}\right)$, then by applying Jensen’s inequality to each individual term in the true log marginal likelihood, we obtain the lower bound

[TABLE]

The quantity $p\left(\bm{u}|\bm{y}_{1:k-1},\bm{\theta}\right)$ is unknown, however, we can replace it with $q_{k-1}(\bm{u})$ leading to $\mathcal{L}\left(q_{1},\ldots,q_{K},\bm{\theta}\right)$

[TABLE]

Maximizing this lower bound recursively with respect to the distributions $q_{k}(\bm{u})$ leads to a sequence of optimal variational distributions $q_{k}^{*}(\bm{u})$ for the inducing outputs

[TABLE]

where $\bm{\Sigma}_{k}=\left(\bm{\Sigma}_{k-1}^{-1}+\frac{1}{\sigma_{n}^{2}}\bm{H}_{k}^{T}\bm{H}_{k}\right)^{-1}$ and $\bm{H}_{k}=\bm{K}_{\bm{X}_{k}\bm{R}}\bm{K}_{\bm{R}\bm{R}}^{-1}$. Plugging $q_{k}^{*}(\bm{u})$ into (29) yields

[TABLE]

The recursive bound $\mathcal{L}_{REC}(\bm{\theta})$ is equivalent to the batch collapsed bound $\mathcal{L}_{VFE}(\bm{\theta})$ in (5) and it holds $\mathcal{L}_{SVGP}(\bm{\mu},\bm{\Sigma},\bm{\theta})\leq\mathcal{L}_{REC}(\bm{\theta})$ with equality when inserting the optimal variational posterior. The variational posterior update in (30) has the same form as the recursive updates in (15). Similarly, the recursive collapsed lower bound in (31) is equal to (25).

We provide below a detailed derivation that follows closely the proof in [36] for the recursive collapsed bound (31) as well as the sequence of optimal distributions (30). We assume mini-batches of size $B$, that is, we have training data $\mathcal{D}=\left\{\bm{y}_{k},\bm{X}_{k}\right\}_{k=1}^{K}$ and the corresponding latent function values$\left\{\bm{f}_{k}\right\}_{k=1}^{K}$. We briefly recap the involved quantities and introduce abbreviations:

[TABLE]

Starting from (29), we have the bound

[TABLE]

which can be rearranged to

[TABLE]

where $\log G(\bm{u},\bm{y}_{k})=\int p\left(\bm{f}_{k}|\bm{u},\bm{\theta}\right)\log p\left(\bm{y}_{k}|\bm{f}_{k},\bm{\theta}\right)\mathop{}\!\mathrm{d}{\bm{f}_{k}}$. The integral involving $\bm{f}_{k}$ is computed as $\log G(\bm{u},\bm{y}_{k})=\int p\left(\bm{f}_{k}|\bm{u},\bm{\theta}\right)\log p\left(\bm{y}_{k}|\bm{f}_{k},\bm{\theta}\right)\mathop{}\!\mathrm{d}{\bm{f}_{k}}$, which equals

[TABLE]

Using $\mathbb{E}\left[\bm{x}^{T}\bm{A}\bm{x}\right]=\text{Tr}\left[\bm{A}\bm{\Sigma}\right]+\bm{\mu}^{T}\bm{A}\bm{\mu}$ with $p\left(\bm{x}\right)=\mathcal{N}\left(\bm{x}|\bm{\mu},\bm{\Sigma}\right)$ yields

[TABLE]

Substitute this expression back, the lower bound becomes

[TABLE]

We can now maximize this bound with respect to $q_{k}(\bm{u})$. Here since we have not constrained $q_{b}$ to belong to any fixed family of distributions, we can compute the optimal bound by reversing the Jensen’s inequality leading to

[TABLE]

where we used a linear Gaussian identity (see, e.g., [4] Ch. 2.3) in the last step . This is equal to Eq. (25). The optimal distribution $q_{k}^{*}$ that gives rise to this bound is proportional to $\mathcal{N}\left(\bm{y}_{k}|\bm{H}_{k}\bm{u},\sigma_{n}^{2}\mathbb{I}\right)q_{k-1}(\bm{u})=\mathcal{N}\left(\bm{y}_{k}|\bm{H}_{k}\bm{u},\sigma_{n}^{2}\mathbb{I}\right)\mathcal{N}\left(\bm{u}|\bm{\mu}_{k-1},\bm{\Sigma}_{k-1}\right)$ and can be analytically computed leading to

[TABLE]

where $\bm{\Sigma}_{k}=\left(\bm{\Sigma}_{k-1}^{-1}+\frac{1}{\sigma_{n}^{2}}\bm{H}_{k}^{T}\bm{H}_{k}\right)^{-1}$. This matches the result in Eq. (15) and (16) and completes the proof.

Appendix C Proof of proposition 1

We showed (App. B) how the batch approximate log-likelihood $\mathcal{L}_{PEP}$ in (4) can be recursively computed. In Sect. 4.1, we showed an equivalent way to compute it, that is, we have the cumulative bound $\psi^{(K)}(\bm{\theta})=\sum_{k}^{K}\psi_{k}(\bm{\theta})$ and it satisfies $\mathcal{L}_{PEP}(\bm{\theta})=\sum_{k}^{K}\psi_{k}(\bm{\theta})$. By induction and linearity of derivatives, we therefore get $\frac{\partial\psi^{(K)}(\bm{\theta})}{\partial\theta}=\frac{\partial\psi_{k}(\bm{\theta})}{\partial\theta}+\frac{\partial\psi^{(K-1)}(\bm{\theta})}{\partial\theta}$ and $\frac{\partial\mathcal{L}_{PEP}}{\partial\theta}=\frac{\partial\psi^{(K)}}{\partial\theta}$.