Splitting Gaussian Process Regression for
Streaming Data

In preparation for the proceeding material, we define some notation. We assume some familiarity of readers with the theory of GPs and kernel methods [24].

The input data matrix and the response vector associated with the $i^{th}$ local model are denoted by
$X_{i}=[\mathbf{x}_{i1}\,\mathbf{x}_{i2}\,\ldots\,\mathbf{x}_{in_{i}}]^{T}\in\mathbb{R}^{n_{i}\times M}$ and $Y_{i}=\left(y_{i1},y_{i2},\ldots,y_{in_{i}}\right)\in\mathbb{R}^{n_{i}}$, respectively, where $\mathbf{x}_{ij}\in\mathcal{X}:=\mathbb{R}^{M}$ is a column vector and $n_{i}$ is the number of observations associated with the $i^{th}$ local model. We call $\mathcal{X}$ the input space.

When creating the $i^{th}$ local GP, its center is defined to be the centroid of $X_{i}$ as follows:

$$\mathbf{c}_{i}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\mathbf{x}_{ij}.$$

Centers are critical to the assignment of observations to local GPs, as well as prediction, in the splitting GP model. Each time a new observation is assigned to a local model, its center is recomputed.

The kernel function $k\!\!:\!\mathcal{X}\times\mathcal{X}\mapsto\mathbb{R}$ is a positive-definite, symmetric function. A vector $\boldsymbol{\theta}$ parameterizes the kernel, and is called the hyperparameter of the GP model. For simplicity, we omit $\boldsymbol{\theta}$ from our notation. We use the term feature space to refer to the reproducing kernel Hilbert space in which $k$ implicitly computes an inner product.

We write $f\sim\mathcal{GP}(\mu,k)$ to say that the function $f$ is distributed as a GP with mean $\mu$ and kernel function $k$, and implicitly assume that the domain of $f$ is $\mathcal{X}$. The notation $f|X_{i},Y_{i}$ denotes that $f$ is conditioned on the data $X_{i},Y_{i}$. We use $\mathbf{x}^{*}$ to denote a test input at which a prediction is to be made and $f^{*}|\mathbf{x}^{*}$ to be a GP conditioned on the test input, i.e. the posterior distribution.

To address modeling a potentially non-stationary latent function, each local GP model is split in a manner which centers the resulting child GPs on regions of different means. This is done by performing principal component analysis (PCA) on the training inputs associated with the parent model. The first principal component gives the direction of most variance of the training inputs, which implies that the corresponding orthogonal hyperplane is the minimizer of within-cluster variance among all linear bisections, leading to Principal Direction Divisive Partitioning (PDDP) [1].

Two new GPs, which we call child GPs, are then created, each of which is assigned the data of the parent model from one side of the hyperplane, as in Fig. 1(a). This heuristic is based on the idea that, by fitting separate GPs to the subsets of training inputs which are maximally different, the model may best adapt to non-stationary behavior of the latent function. We utilize PDDP since the splitting procedure is computed efficiently in closed form, i.e. without the use of a convergent algorithm such as k-means. It also allows for different choices in how the principal direction is computed. For example, in a setting where observations are fully sequential, the principal direction can be efficiently approximated using Oja’s rule [19] to avoid a singular value decomposition of the data matrix.

Critically, the prior of each child GP is taken to be the posterior of the parent model when conditioned on the $m$ training observations, reflecting the Bayesian belief that the latent function’s local covariance structure is not too dissimilar from its covariance structure in the larger subspace prior to splitting.

Formally, in the function space inference view of GP regression,

$$f_{child}\sim f_{parent}|X_{parent},Y_{parent},$$

(1)

where $f_{child}$ is the prior of the child and $f_{parent}|X_{parent},Y_{parent}$ is the posterior of the parent given the training data $X_{parent},Y_{parent}$. This assumption of local similarity is implicit in the use of many smooth kernel functions, particularly the radial basis function family, which are infinitely differentiable.

A new data point $\left(\mathbf{x},y\right)$ is assigned to a child GP whose center is most similar to the predictor $\mathbf{x}$ in the feature space, as determined by the kernel function. That is, for $C$ child GPs, the data point is assigned to the child GP indexed by

$$i_{assign}:=\underset{i=1,...,C}{\arg\!\max}\,k(\mathbf{c}_{i},\mathbf{x}).$$

Similarly, the posterior mean at the test input $\mathbf{x}^{*}$ is computed by weighting the prediction of each child GP by the relative similarity of its center to the test input in the feature space. This idea is based upon an interpretation of the predictive mean as weighting observations by their similarity to the point of prediction in the feature space, according to Shen et al. [25]. The same idea was also used by Nguyen-Tuong et al. [18].

Unlike these previous works, we take a weighted average of all predictions of child GPs, rather than only a subset, as follows:

$$\mathbb{E}[\,f^{*}|\mathbf{x}^{*}\,]=S^{-1}\sum_{i=1}^{C}k(\mathbf{c}_{i},\mathbf{x}^{*})\mathbb{E}[\,f^{*}_{i}|X_{i},Y_{i},\mathbf{x}^{*}\,],$$

(2)

where $S=\sum_{i=1}^{C}k(\mathbf{c}_{i},\mathbf{x}^{*}).$ This may be interpreted as weighting each child GP’s mean prediction by the similarity of its local region to the test input $\mathbf{x}^{*}$. In the next section, we will show that using predictions from all child GPs has important consequences on the theoretical properties of the resulting model.

As a consequence of this weighting procedure and the choice of priors for the child GPs, the splitting GP model has some convenient properties. The proofs of the following propositions may be found in the appendix.

This result is somewhat intuitive, since no additional evidence has been obtained which might alter our Bayesian beliefs about the latent function; we have simply altered the model structure. In agreement with this idea, the equality of the parent and childrens’ predictive distribution no longer holds after any new data is assigned to one of the child GPs. Note that, in general, this relationship does not hold for any split after the first.

An important property of the splitting GP model is the preservation of continuity properties of the child GPs. We provide a result which shows the necessary and sufficient condition for continuity of the splitting GP model’s predictions.

While it may seem obvious that the mean prediction is continuous in the input space, it is reassuring to know that the random field created by the splitting model has the same sufficient condition for mean square continuity as the underlying child GPs. Intuitively, we are computing a continuously weighted average of smooth random fields, so the resulting random field should also be smooth.

Note that it is not necessary to aggregate the predictions of each local model. Instead, one may elicit predictions only from the $C_{0}<C$ local models most similar to the point of prediction, as measured by the kernel. This variation of the prediction method, as used in [18], may result in faster computation of predictions. However, we caution against this practice without careful consideration. This prediction method will result in a loss of continuity of the mean prediction $\mathbb{E}[\,f^{*}|\mathbf{x}^{*}\,]$, and consequently the mean square continuity of the random field $f^{*}|\mathbf{x}^{*}$, since the mean prediction is computed using a maximum function, which is not continuous. The discontinuity will manifest as sharp jumps in the mean predictions, as illustrated in Fig. 1(b).

The algorithm for the splitting GP model has improved complexity in both memory and training/prediction time compared to the full GP regression. Additionally, it maintains some benefit in asymptotic complexity relative to other local GP methods.

The complexity of updating the splitting GP model with a single datum is bounded above by $\mathcal{O}(m^{3})$, corresponding to a matrix inversion of one child GP. The contribution of the PCA-based splitting procedure to the time complexity is negligible. Assuming the PCA is performed via naive singular value decomposition, the procedure would have $\mathcal{O}(m^{2})$ complexity amortized over $m$ sequential observations, which yields a linear additive term $m<m^{3}$.

While each child GP has at most $m$ observations, the average case will clearly be lower. It should be noted that we explicitly chose to not utilize a rank-one update of the Cholesky decomposition of the kernel matrix during the update procedure, a method which is used in other local GP methods such as in [18]. We made this choice to permit updating the model with a batch of observations, in addition to fully sequential updating, which would require an update of rank greater than one. We also observed empirically that repeatedly updating a single child GP using rank-one updates would cause numerical issues if the kernel length-scale parameter is small.

Since the time complexity of the algorithm is characterized primarily by the parameter $m$, the largest number of observations which may be associated with a single child model, it is straightforward to select an appropriate parameter value for applications requiring rapid sequential updating. After empirically determining the necessary wall-clock time needed for an update, the parameter may be adjusted appropriately.

The splitting algorithm also imposes a lower memory complexity in terms of the number of observations, $n$. It is well known that local GP methods effectively store only a block-diagonal approximation of the full covariance matrix [20, 3]. In particular, when the splitting algorithm has $n$ observations, $\lfloor n/m\rfloor$ child GPs have been created, each of which will store a kernel matrix with up to $m^{2}$ entries. The resulting memory complexity of the algorithm is then $\mathcal{O}(mn)<\mathcal{O}(n^{2})$, where $\mathcal{O}(n^{2})$ is the memory complexity of a full GP regression. In contrast, the rBCM by [4] has an asymptotic memory complexity of $\mathcal{O}(n^{2}/E)$, where $E$ is the (constant) number of experts specified as a parameter. Nguyen-Tuong et al. [18] did not present the memory complexity of their local GP model but, under the mild assumption that the input space is bounded, it can be shown that the asymptotic memory complexity is $\mathcal{O}(n^{2})$. Notably, the splitting GP model is, to the best of our knowledge, the only local GP model to achieve a linear memory complexity.

The synthetic data set (see Fig. 2) was constructed to have a non-stationary latent function $f$ of a two-dimensional predictor $(x_{1},x_{2})\in[-1,1]^{2}$. The response function is defined as

	$\displaystyle y$	$\displaystyle=f(x_{1},x_{2})+\epsilon$
		$\displaystyle=5\sin(x_{1}^{2}+x_{2}^{2})+3x_{1}+\epsilon,$		(3)

where $\epsilon\overset{iid}{\sim}\mathcal{N}(0,\sigma)$ with $\sigma=(0.05)\underset{x_{1},x_{2}}{\max}\,f(x_{1},x_{2})$. To construct the synthetic data set, a grid of 10,000 points was constructed in the $x_{1}$-$x_{2}$ plane and the latent function evaluated at each point. A subset of 2500 observations were sampled uniformly, without replacement, from the resulting grid for each replicate of the experiment. This sampling over the grid ensured that no two observations are too close to one another, which may cause numerical issues during training.

For this experiment, we compared three metrics of the models’ performance: prediction error (in mean squared error, or MSE), memory usage (in kilobytes, or kB), and training time (in seconds). Each metric was estimated using a 5-fold cross validation procedure. We utilized common random numbers [28] as a variance reduction technique. Each experiment was replicated 10 times to reduce variability of results. In Fig. 3, the shaded region shows the 95% pointwise confidence region for the mean metric. The experiment was performed with different numbers of observations, ranging from 100 to 2500 in increments of 100, to demonstrate the relative data requirements for each model to converge. Each model was updated sequentially with single observations to simulate a streaming data setting.

For each alternative model compared, the experiment was replicated for a range of parameter values and the most favorable, in terms of MSE, results for each model are reported. For the splitting GP model, we used a splitting limit of $m=500$. For the rBCM, $E=10$ local experts were used. We found the parameter $w_{gen}$ of the local GP model difficult to tune and eventually used $w_{gen}=10^{-3}$ based on an extensive grid search (detailed in the appendix).

In Fig. 3, one can see that the splitting model and full GP required relatively few observations to achieve strong predictive power. The rBCM was much slower to converge, and exhibited high variability in its MSE across replicates of the experiment. We attribute the greater variability in MSE to the rBCM’s random assignment of data to GP experts. On the other hand, the splitting GP model exhibited remarkably low variability in MSE between replicates. It is worth noting that the splitting model’s MSE slowly increased as it splits at intervals of 500 observations (see Fig. 4).

The memory usage of the splitting GP model proved to be significantly lower than the full GP model, and marginally higher than that of the rBCM. However, for ~2500 observations, it can be seen that the memory usage of the rBCM model began to surpass that of the splitting GP model. This is to be expected, since the asymptotic memory complexity of the rBCM is quadratic, as opposed to the linear complexity of the splitting GP model.

The training time of the rBCM and local GP models were found to be comparable, and the splitting GP model took slightly longer. The full GP took significantly longer to train. The change in regime of the training time of the full GP at ~1100 observations is due to specialized numerical methods in the GPyTorch library [6].

The first real-world data set, kin40k, consists of 40,000 records describing the location of a robotic arm as a function of an 8-dimensional control input [17]. This data set was chosen since it exemplifies a task where the fast sequential updating and low memory profile offered by the splitting GP model are desirable. Furthermore, kin40k is a popular benchmark for other GP regression methods [4, 16, 14], and thus facilitates direct comparison.

We used the same training/test split of 10,000/30,000 points as in the previous work [4, 16, 14]. Observations were added to the models in batches, with each batch having a number of observations equal to the number of observations per local model. The results can be seen in Fig. 5(a). The splitting model’s root-mean-square error (RMSE) is comparable with that of the rBCM, which is designed under a stronger assumption (that the latent function is well-modeled by a single GP), which is satisfied in this stationary regression task. In contrast, the local GP model struggled to achieve the same RMSE in this experiment (see the appendix for more detail on the difficulty tuning the parameter $w_{gen}$).

The second real-world data set is powergen, which consists of four predictors describing control inputs of a combined cycle power plant, and a response of the net electrical energy production. This data set is due to Kaya et al. [12] and Tüfekci [27] and is publicly available [5]. We pre-processed the data to remove duplicate observations, leaving a total of 7622 observations.

In the powergen experiment, we compared both the predictive error (in RMSE) and the combined training and prediction time (in seconds) of the splitting GP model and rBCM. We compared only these two models, since they proved to be the most competitive local GP methods in earlier experiments. The data set was randomly divided into a 80%/20% train/test split, and each model evaluated for several different parameter values. Observations were added to the models in batches, with each batch having a number of observations equal to the number of observations per local model. The experiment was replicated 10 times, and common random numbers were used for splitting the data set and training the models for sharper comparisons between the models. The results can be seen in Fig. 5(b).

The splitting GP model’s RMSE decreased with the number of observations per local model, achieving a comparable performance with the rBCM. In terms of runtime, the splitting GP model was significantly faster, while having little variability between replicates. While the rBCM’s runtime decreased with the number of observations per expert, its average runtime was 2.5-6 times longer than that of the splitting GP model, depending on the choice of parameters.