Fast Gaussian Process Regression for Big Data
Abstract
Gaussian Processes are widely used for regression tasks. A known limitation in the application of Gaussian Processes to regression tasks is that the computation of the solution requires performing a matrix inversion. The solution also requires the storage of a large matrix in memory. These factors restrict the application of Gaussian Process regression to small and moderate size data sets. We present an algorithm that combines estimates from models developed using subsets of the data obtained in a manner similar to the bootstrap. The sample size is a critical parameter for this algorithm. Guidelines for reasonable choices of algorithm parameters, based on detailed experimental study, are provided. Various techniques have been proposed to scale Gaussian Processes to large scale regression tasks. The most appropriate choice depends on the problem context. The proposed method is most appropriate for problems where an additive model works well and the response depends on a small number of features. The minimax rate of convergence for such problems is attractive and we can build effective models with a small subset of the data. The Stochastic Variational Gaussian Process and the Sparse Gaussian Process are also appropriate choices for such problems. These methods pick a subset of data based on theoretical considerations. The proposed algorithm uses bagging and random sampling. Results from experiments conducted as part of this study indicate that the algorithm presented in this work can be as effective as these methods.
keywords:
Big Data, Gaussian Process, RegressionMsc:
[2010] 0001, 9900tabular
[cmiauthor1]Chennai Mathematical Institute \fntext[isiauthor]Indian Statistical Institute
1 Introduction
Gaussian Processes (GP) are attractive tools to perform supervised learning tasks on complex datasets on which traditional parametric methods may not be effective. They are also easier to use in comparison to alternatives like neural networks (Rasmussen06gaussianprocesses ()). Gaussian Processes offer some practical advantages over Support Vector Machines (SVM) (ghahramani_mlss_2011 ()). They offer uncertainty estimates with predictions. The kernel and regularization parameters can be learned directly from the data instead of using cross validation. Feature selection can be incorporated into the learning algorithm. For regression, exact inference is possible with Gaussian Processes. To apply Gaussian Processes to classification, we need to resort to approximate inference techniques such as Markov Chain Monte Carlo, Laplace Approximation or Variational Inference. Even though exact inference is possible for Gaussian Process regression, the computation of the solution requires matrix inversion. For a dataset of size , the time complexity of matrix inversion is . The space complexity associated with storing a matrix of size is . This restricts the applicability of the technique to small or moderate sized datasets.
In this paper we present an algorithm that uses subset selection and ideas borrowed from bootstrap aggregation to mitigate the problem discussed above. Parallel implementation of this algorithm is also possible and can further improve performance.
The rest of this paper is organized as follows: In section 2, we discuss the problem context. In section 3, we present our solution to the problem. Our solution is based on combining estimators developed on subsets of the data. The selection of the subsets is based on simple random sampling with replacement (similar to what is done in the bootstrap). The size of the subset selected is a key aspect of this algorithm. This is determined empirically. We present two methods to determine the subset size. We present the motivating ideas leading to the final form of the algorithm. When the model is an additive structure of univariate components, this has attractive implications on the convergence rate (stone1985additive ()). An additive model worked well for the datasets used in this study. Relevant facts from Minimax theory for nonparametric regression, that are consistent with the experimental results reported in this work, are presented. In section 4, we present a brief summary of related work. Applying Gaussian Processes to large datasets has attracted a lot of interest from the machine learning research community. Connecting ideas to research related to the algorithm reported in this work are presented.
Selecting parameters for an algorithm is an arduous task. However this algorithm has only two important parameters, the subset size and the number of estimators. We provide guidelines to pick these parameters based on detailed experiments across a range of datasets. In section 5 we provide experimental results that provide insights into the effect of the parameters associated with the algorithm. In section 6, we illustrate the application of our algorithm to synthetic and real world data sets. We applied the algorithm developed in this work to data sets with over a million instances. We compare the performance of the proposed method to the Sparse Gaussian Process (titsias2009variational ()) and the Stochastic Variational Gaussian Process (hensman2013gaussian ()). The inputs required by these algorithms are similar to the inputs required for the proposed method and therefore are applicable in a similar context.
We compare the estimates obtained from the reported algorithm with two other popular methods to perform regression on large datasets  Gradient Boosted Trees (using XGBoost, chen2016xgboost ()) and the Generalized Additive Model (GAM)(friedman2001elements ()).Results from experiments performed as part of this study show that accuracies from the proposed method are comparable to those obtained from Gradient Boosted Trees or GAM’s. However there are are some distinct advantages to using a Gaussian Process model. A Gaussian Process model yields uncertainty estimates directly whereas methods like Gradient Boosted Trees do not provide this (at least directly). A Gaussian Process model is also directly interpretable in comparison to methods like Gradient Boosted Trees or Neural Networks. Therefore, the proposed method can yield both explanatory and predictive models. It is possible to use stacking (wolpert1992stacked ()) to combine the estimates from the proposed model with those obtained from a competing model (like Gradient Boosted Trees) and obtain higher accuracies. Combining a Gaussian Process solution with XGBoost has been used by lloyd2014gefcom2012 ().
In section 7, we present the conclusions from this work. The contribution of this work is as follows. The proposed method to perform Gaussian Process regression on large datasets has a very simple implementation in comparison to other alternatives, with similar levels of accuracy. The algorithm has two key parameters  the subset size and the number of estimators. Detailed guidelines to pick these parameters are provided. The choice of a method to scale Gaussian Process regression to large datasets depends on the characteristics of the problem. This is discussed in section 4. The proposed method is most effective for problems where the response depends on a small number of features and the kernel characteristics are unknown. In such cases, exploratory data analysis can be used to arrive at appropriate kernel choices duvenaud2014automatic (). Additive models may work well for these problems. Appropriate preprocessing like principal component analysis can be used if needed to create additive models. The rate of convergence for additive models is attractive (stone1982optimal ()). This implies that we can build very effective models with a small proportion of samples. Sparse Gaussian Processes see snelson2005sparse () and Stochastic Variational Gaussian Processes hensman2013gaussian () are also appropriate for such problems. These require a more complex implementation and may require extensive effort to tune the optimization component of the algorithm (see hensman2013gaussian ()). Results of the experiments conducted as part of this study show that the proposed method can match or exceed the performance of these methods.
2 Problem Formulation
A Gaussian Process with additive noise can be represented as:
(1) 
Here :

represents the observed response.

represents an input vector of covariates.

represents the noise. The noise terms are assumed to be identical, independently distributed (IID) random variables drawn from a normal distribution with variance .

represents the function being modeled. It is a multivariate normal with mean function and covariance function .
If we want to make a prediction for the value of the response at a test point , then the predictive equations are given by (see Rasmussen2005 ()):
(2) 
(3) 
(4) 
Here:

is the value of the function at the test point .

represents the covariance between the test point and the training set.

represents the covariance matrix for the training set.

I is the identity matrix.
Equation (3) is the key equation to make predictions. An inspection of equation (3) shows that this equation requires the computation of an inverse. For a training dataset of size , computation of the inverse has time complexity. This is one of the bottle necks in the application of Gaussian Processes. Calculation of the solution also requires the storage of the matrix . This is associated with a space complexity of . This is the other bottle neck associated with Gaussian Processes. The uncertainty associated with our prediction is provided by Equation 4. The covariance is expressed in terms of a kernel that is appropriate for the modeling task. The kernel is associated with a set of hyperparameters. These may be known for example if the modeling task has been well studied or unknown if the problem has not been well studied. In this work we treat these kernel parameters as unknown. When the kernel parameters are unknown, these are estimated using maximum likelihood estimation. The marginal log likelihood is given by (see Rasmussen2005 ()):
(5) 
Using an appropriate optimization technique with Equation (5) as the objective function, the hyperparameters of the kernel () can be determined.
3 Proposed Solution
Since GP regression is both effective and versatile on complex datasets in comparison to parametric methods, a solution to the bottlenecks mentioned above will enable us to apply GP regression to large datasets that are complex. We run into such datasets routinely in this age of big data. Our solution to applying Gaussian Process regression to large data sets is based on developing estimators on smaller subsets of the data. The size of the subset and the desired accuracy are key parameters for the proposed solution. The accuracy is specified in terms of an acceptable error threshold, . We pick smaller subsets of size from the original data set (of size ) in a manner similar to bootstrap aggregation. We develop a Gaussian Process estimator on each subset. The GP fit procedure includes hyperparameter selection using the likelihood as the objective function. We want a subset size such that when we combine the estimators developed on the subsets, we have an estimator that yields a prediction error that is acceptable. The rationale for this approach is based on results from Minimax theory for nonparametric regression that are presented later in this section. To combine estimators, we simply average the prediction from each of the estimators.
The time complexity for fitting a Gaussian Process regression on the entire dataset of size is . The algorithm presented above requires the fitting of Gaussian Process regression models to smaller size datasets (). Therefore, the time complexity is .
We present two methods to determine the subset size:

Estimating the subset size using statistical inference. For a dataset of size , the subset size is expressed as:
(6) is a random variable and is determined based on inference of a proportion using a small sample. The details of this method are presented in the next subsection.

Estimating the subset size using an empirical estimator. We use the following observations to derive this empirical estimator:

The subset size , should be proportional to the size of the dataset, :
We posited that as the size of the dataset increases, there is possibly more detail or variations to account for in the model. Therefore larger datasets would require larger sample sizes.

The sample size , should be a decreasing function of the desired error rate (). This means that decreasing the desired error rate should increase the sample size.
where is an increasing function. Application of the above observations yields the following estimator for sample size:
(7) Here:

is a function that characterizes the fact that an increase in should increase the sample size .

is a function that characterizes the fact that sample size should increase as (desired error level) decreases.

The algorithm is summarized in Algorithm 1
The proposed algorithm is based on model averaging as described in (bishop2006pattern, , Chapter 14). This is similar to combining estimates from regression trees in the Random Forest(breiman2001random ()) algorithm. Given models, the conditional distribution of the response at a point is obtained from:
(8) 
If each of the models is equally probable, then . If each of is a Gaussian , as would be the case when each of the estimators is a GP based on Equation 1, then we have the following:
(9) 
Where:
Model combination is an alternative approach to combining estimates from a set of models. The product of experts model is an example of this approach. In the product of experts approach, each model is considered to be an independent expert. The conditional distribution of the response at a point is obtained from:
(10) 
where is the normalization constant. cao2014generalized () report a study where the individual experts are Gaussian Processes. In this case, each of are . The mean and variance associated with are:
(11)  
(12)  
(13) 
Here is the precision of the expert at point . The results from using a product of experts model are also included in this study (see section 6).
Model averaging and model combination (Equation (9) and Equation (11)) are simply ways to combine estimates from the component estimators used in Algorithm 1. They do not specify any information about the character of the estimator developed using Algorithm 1. To do this, we need the following preliminaries:
Assumption 1.
The function being modeled is in the Reproducing Kernel Hilbert Space of the kernel () of the estimator. The reproducing property of the kernel can be expressed as:
Here represents the index set to select the predictor instances.
Algorithm 1 makes use of estimators, , , , . The reproducing kernels associated with the estimators are , , , . We show that the estimator resulting from combining these estimators is a Gaussian Process.
Lemma 1.
The estimator obtained from the individual estimators , , , using:
is associated with the following kernel:
(14) 
Proof.
From Assumption 1, the kernel associated with each estimator has the reproducing property. So the following equation holds:
Consider the kernel and the inner product . The inner product can be written as:
∎
Using Lemma (1), the estimator developed using Algorithm 1 is a Gaussian Process. This Gaussian Process is associated with a kernel defined by Equation (14). This kernel is defined as a mixture of the kernels used by the individual estimators. The uncertainty estimate at any point can be obtained using Equation (4).
Consistency of estimators is concerned with the asymptotic properties of the developed estimator as we increase the sample size. Consistency of Gaussian Processes has been widely studied (see (Rasmussen06gaussianprocesses, , Chapter 7, Section 7.1)) for details. The convergence rate for the algorithm is a characteristic of great interest. This characteristic tells us the rate at which we can drop the error as we increase the sample size. Minimax theory provides a framework to assess this (see tsybakov2009introduction () or gyorfi2006distribution ()) . A famous result by stone1982optimal () states the minimax rate of convergence for nonparametric regression is , where describes the smoothness of the function being modeled and represents the dimensionality of the feature space. This suggests that the rate of convergence is affected by the curse of dimensionality. However as noted by yang2015minimax (), the following factors usually mitigate the curse of dimensionality:

The data lie in a manifold of low dimensionality even though the dimensionality of the feature space is large. yang2016bayesian () reports a method for such problems.

The function being modeled depends only on a small set of predictors. All datasets reported in this study had this characteristic. Feature selection can be performed in Gaussian Processes using Automatic Relevance Determination (see (Rasmussen06gaussianprocesses, , Chapter 5, Section 5.1)).

The function being modeled admits an additive structure with univariate components. The minimax rate of convergence for this problem is very attractive (see stone1985additive ()). For all datasets reported in this study such an additive model yielded good results.
Feature Selection is therefore a critical first step. This is discussed in section 6. The additive structure and the dependence on a very small set of predictors suggest that we can get reasonable models with a small subset of the data. This was consistent with the experimental results reported in this work. Simple preprocessing like Principal Component Analysis (PCA) could be used to reduce the number of relevant features. PCA also makes the features independent. The data for this study came from public repositories and from very diverse application domains. This suggests that datasets with these characteristics are not uncommon. When data lie in a manifold, methods such as yang2016bayesian () may be more appropriate.
zaytsev2017minimax () investigate the minimax interpolation error under certain conditions (known covariance, stationary Gaussian Process) . As noted in this recent study, theoretical work in estimating the convergence rate using minimax theory is an active area of research. In this work we investigated empirical estimation of the sample size. We describe two methods to determine the sample (subset) size.
3.1 Estimating Subset Size on the Basis of Inference of a Proportion
Since takes values between 0 and 1, it can be interpreted as a proportion. We treat it as a random variable that can be inferred from a small sample. To estimate , we do the following: We pick a small sample of the original dataset by simple random sampling. We start with a small value of and check if the prediction error with this value of is acceptable. If not we increment until we arrive at a that yields an acceptable error. This procedure yields the smallest value that produces an acceptable error on this sample. Since the size of this dataset is small, the above procedure can be performed quickly. This technique yielded reliable estimates of on both synthetic and real world datasets.
3.2 Empirical Estimation of the Subset Size
4 Related Work
(Rasmussen2005, , Chapter 8) provides a detailed discussion of the approaches used to apply Gaussian Process regression to large datasets. quinonero2007approximation () is another detailed review of the various approaches to applying Gaussian Processes to large datasets. The choice of a method appropriate for a regression task is dependent on the problem context. Therefore we discuss the work related to scaling Gaussian Process regression taking the problem context into consideration.
If the data associated with the problem has been well studied and kernel methods have been successfully applied to the problem, then we may have reasonable insights into the nature of the kernel appropriate for the regression task. We may be able to arrive at the kernel hyperparameters quickly from a small set of experiments. On the other hand if the problem and data is new, then we may not have a lot of information about the kernel. In general, scaling Gaussian Process regression to large datasets has two challenges:

Finding a kernel that has good generalization properties for the dataset.

Overcoming the computational hurdles  for training and for storage.
Learning a kernel that has good generalization properties is a related area of research in Gaussian Processes (see wilson2014thesis (), wilson2013gaussian (). When a good kernel representation has been learned, there are many techniques to overcome the computational hurdles. The Nystrom method to approximate the Gram matrix (drineas2005nystrom ()) and the Random Kitchen Sinks (rahimi2008random ()) are probably the most well known. The Random Kitchen Sinks approach maps the data into a low dimensional feature space and learns a linear estimator in this space. It should be noted that these methods work well when the problem needs a stationary kernel for which we know the hyperparameters. Using ”sensible defaults” for hyperparameters and applying these techniques to problems that require a nonstationary kernel may yield poor results. For example with the airline dataset, described later in this study (see section 5.1), the Random Kitchen Sinks cannot be used directly and would require suitable preprocessing (like removing simple trends or using a mean function) so that a stationary kernel would be applicable. Using the Random Kitchen Sinks directly with no preprocessing and using default kernel choices provided with the scikitlearn scikitlearn () implementation yielded poor results (RMSE of 31.49 as opposed to 8.75 with the proposed method).
Using a kernel learning approach to determine a good kernel representation and then solving the computational hurdles independently is one way to approach scaling Gaussian Process regression to large datasets. Another approach to determine the appropriate kernel is to use exploratory data analysis. Guidelines to pick kernels based on exploratory analysis of the data is provided in duvenaud2014automatic (). This is a practical approach when the number of relevant features is not too many, as was the case with the datasets used in this study. It should be noted that hyperparameters for these choices still need to be specified. We may be able to build additive models using this approach. Appropriate preprocessing could help, for example principal component analysis can be applied to make the features independent. Minimax theory for nonparametric regression indicates that the convergence rate for additive models is very attractive. We can build effective models with a small proportion of the data.
The choice of kernel hyperparameters is critical and can affect the performance. When datasets are large and the kernel hyperparameters are unknown, we need algorithms that can address both these issues. Ideally, the algorithm should be able to work with both stationary and nonstationary kernels. The proposed algorithm is one such candidate.
Sparse Gaussian Processes (titsias2009variational ()) and Stochastic Variational Gaussian Processes (hensman2013gaussian ()) are two others. Like the proposed algorithm, these algorithms require the specification of a input size. A subset of points is selected from the dataset for training. The criteria for subset selection is different in each case. These algorithms do not require the specification of the kernel hyperparameters. These are estimated from the data. Stochastic Variational Gaussian Processes can require considerable manual tuning of the optimization parameters. Typical implementations (like gpy2014 ()) for Sparse GP and Stochastic Variational GP use stochastic gradient descent for hyperparameter optimization. This explores the entire dataset in batch size increments. hensman2013gaussian () report the details associated with picking the parameters for the optimization task (learning rates, momentum, batch sizes etc.). In contrast, sample sizes with the proposed algorithm even for datasets with over million instances are typically small (order of few hundred instances). Learning hyperparameters over small datasets is considerably easier. The experiments reported in this work required no tuning effort. We report the performance of Sparse Gaussian Process, Stochastic Variational Gaussian Process and the proposed method on a variety of datasets in section 6
The proposed algorithm uses ideas that have proven effectiveness with other machine learning techniques. Bagging has been used to improve performance using regression trees (Random Forests, breiman2001random ()). Like with Random Forests, the algorithm uses model averaging to combine estimates from component Gaussian Process Regressions. Dropout (srivastava2014dropout ()) is a technique used in neural networks to prevent over fitting. Dropping random units achieves regularization in neural networks. In the proposed algorithm, selecting a sample can be viewed as dropping random instances from the training dataset. Sparse Gaussian Processes and Stochastic Variational Gaussian Processes use theoretical ideas to select a small subset of points to develop a Gaussian Process regression model. This study suggests that combined with model averaging, random selection of the subset can also work well.
For datasets with a Cartesian product structure, by imposing a factorial design of experiment scheme on the dataset and decomposing the covariance matrix as a Kronecker product belyaev2016computationally () discuss an approach to scaling Gaussian Process regression. This approach also uses a prior on the hyperparameters to deal with anisotropy. So the implementation is quite complex. Using a divide an conquer strategy is another theme in scaling Gaussian Process regression to large datasets. The Bayesian Committee Machine (BCM) (tresp2000 ()), is an idea related to the algorithm presented in this work. The BCM proposes a partition of the dataset into parts. An estimator is developed on each partition. The BCM does not choose a subset of the partition. It uses the entire partition for developing the estimator. This is the key difference between the method proposed in this work and the BCM. The estimates from each estimator are assumed to be independent. The BCM assumes that computing a GP solution on the partitions of the dataset is computationally tractable because the partition sizes are small. Datasets encountered today are much larger than those reported in tresp2000 (). In present day datasets the partitions of the dataset based on guidelines provided in tresp2000 () would be very big and computing a full Gaussian Process solution on them may not be computationally tractable. Even when the partition size is not big enough to create computational hurdles, using all the data may result in over fitting. Using a hierarchical model as in deisenroth2015distributed () or park2010hierarchical () are possible ways to work around the size of the partitions, however this requires a complex implementation to partition and recombine computations.
The Locally Approximate Gaussian Process gramacy2015local () fits a local Gaussian Process for a prediction point using a local neighborhood and local isotropy assumption. This method too requires some tuning, the size of neighborhood and method to choose the neighbors are important parameters. For datasets where the local isotropy assumption works well and when the size of the test set is small, this method might be useful. When prediction is required at a large number of test points, this method might be slow if we use a large neighborhood. For example with the airline dataset, described in section 5.1, using the defaults provided with the laGP (laGPR ()) package did not yield good results on the airline delay dataset (RMSE of 24.89 as opposed to 8.75 with the proposed method). The running time for laGP was also considerably longer. Scoring the test set in batches of 15000 rows, the test set prediction took about 50 minutes for the airline delay dataset. The proposed algorithm builds a single model that is used to score the entire test set and completed in about 6.5 minutes.
In summary, there are several ways to scale Gaussian Process regression to large datasets. The choice of a particular method should be guided by the characteristics of the problem. The method proposed in this work is appropriate for large datasets that have a small number of important features for which the kernel characteristics are unknown. In such cases exploratory data analysis can be used to determine appropriate kernel types. Additive models may work well for such datasets. Preprocessing such as principal component analysis can be used if needed to make features independent (so that we can use additive models). Stochastic Variational Gaussian Processes and Sparse Gaussian Processes are also good candidates for such problems. Kernel hyperparameters are learned from the data by these methods. Results of the experiments conducted as part of this study show that the proposed method can match or exceed the performance of the Sparse Gaussian Process or the Stochastic Variational Gaussian Process.
5 Effect of the Parameters
Selection of algorithm parameters appropriate for a machine learning task is an arduous task for all practitioners. To alleviate this difficulty, we provide guidelines for parameter selection based on detailed experimentation. The proposed algorithm has three parameters:

The dataset size

The subset size

The number of estimators
Accordingly, three sets of experiments were performed to capture the effect of each of these parameters on the performance of the algorithm. These experiments are described in this section.
5.1 Datasets
The following datasets were used in this study:

Combined Cycle Power Plant: This dataset was obtained from the UCI Machine Learning repository (Lichman ()). This dataset has 9568 instances. The target variable is the net hourly electrical power output from a power plant. The dataset has four features.

Ailerons: This dataset was obtained from the LIAD(Laboratory of Artificial Intelligence and Decision)(LIACC ()). The target variable for this dataset is the control action associated with the control of a F16 aircraft. The dataset has 40 features and 7154 instances.

Elevators: This dataset was obtained from the LIAD repository (LIACC ()). This dataset is also related to the control of a F16 aircraft. The target for this dataset is the control action variation for the elevators of the aircraft. The dataset has 6 features and 9517 instances.

California Housing: This dataset was obtained from the LIAD repository (LIACC ()). The target variable for this dataset is the median house price. The dataset has 8 features and 20460 instances

Individual Household Electric Power Consumption: This dataset was obtained from the UCI Machine Learning repository (Lichman ()). It captures the electric power consumption in a single household over a four year period at a one minute sampling rate. For the experiments in this study, the Voltage was treated as the target variable. Seasonality and periodicity are important characteristics of this dataset (identified during exploratory data analysis of this dataset). For this reason, minute and hour attributes were created out of the time stamp attribute. Similarly, day of week and day of month attributes were created out of the date attribute. This dataset has over 2 million instances and 12 features.

Airline Delay:This dataset was obtained from the US Department of Transportation’s website (RITA_Delay_Data_Download ()). The data represents arrival delays for US domestic flights during January and February of 2016. This dataset had 12 features and over two hundred and fifty thousand instances. Departure delay is included as one of the predictors while hensman2013gaussian () does not include it. Also the raw data includes a significant amount of early arrivals (negative delays). For all regression methods considered in this study, better models were obtained by limiting the data to the delayed flights only (arrival delays were greater than zero). This suggests that arrival delays and early arrivals are better modeled separately.

The Sinc Function: This is a one dimensional synthetic dataset where the response variable is the sine cardinal function otherwise called the sinc function (noise free). The sinc function is a complex function to learn and is therefore a candidate for this as well as many other machine learning research studies. The dataset had one hundred thousand instances.
5.2 Effect of Dataset Size
These experiments study the effect of the size of the dataset () on the subset size (). For each dataset, we pick a fraction of the data elements and determine the subset size () required to achieve a target accuracy. The target accuracy is an input parameter for the experiment. We repeat this procedure for various settings of the fraction of the dataset selected ( through ). The number of estimators for these experiments was maintained at 30. The rationale for this choice is provided in section 5.4. The results are shown in Figure 4 through 7.
The key observation from these experiments was that the required to maintain a preset accuracy decreases very slowly as increases. This set of experiments was used to identify candidate choices for the parameters of Equation 7. Since we wanted a function , that decreases slowly as increases, we considered . These are slowly decreasing functions of in decreasing order of slowness. After a rigorous empirical study we found that if we choose then it works well on all real world datasets and synthetic data. As discussed in section 3, the time complexity of the proposed algorithm is or . The exponent of decreases monotonically as increases. Since can be expressed as , for , the running time is . The number of estimators , determined empirically, is a constant (about ) for all experiments. The rationale for this choice is explained in section 5.4. The time complexity of the GP computation is therefore . So when N is large enough (such that ), the GP computation is sublinear. For example, when , the time complexity is . is a monotonically increasing function that characterizes the fact that sample size should increase as the acceptable error threshold decreases (e.g. ). An optimal choice of is an area of future work. For the experiments reported in this work, with for low noise datasets () and for noisy datasets (), worked well.
5.3 Effect of Subset Size
These experiments explore the effect of the subset size, captured by the parameter , on the accuracy (), for a fixed dataset size (). For each dataset, a fraction of the dataset is picked for this experiment. This represents the data for the experiment (). We pick such that (i.e, ) is computationally tractable (about 2000). The subset , used for Gaussian Process model development was . We fix the number of estimators to be 30 (see section 5.4 for rationale). For each value in a range of values, we record the RMSE (). The key insight from these experiments was that complex functions, like the Sinc function (see Figure 14) needed a larger subset size to produce a given level of accuracy. Note that we can expect the to drop with the increase in dataset size because of the effect reported in section 5.2. Section 6.3 provides the values associated with the full dataset. The results for these experiments are shown in Figure 9 through 14.
5.4 Effect of Number of Estimators
These experiments capture the effect of the number of estimators () on the accuracy of the algorithm (), for a particular set of and values.
For each dataset the following experiment was performed. The fraction of the dataset to use for the experiment () and the subset size () are selected. For each in a range of values, Algorithm 1 is applied and the error () is recorded. The key insight from this set of experiments was that given a subset size(), there is a point upto which increasing the number of estimators () drops the error, but once a threshold value is reached (around , for all the datasets), there is little benefit in increasing the number of estimators. A plausible explanation for this behavior could be that increasing the number of estimators reduces the variance component of the error in the biasvariance decomposition of the error(). The results for these experiments are shown in Figure 16 through 21.
6 Application of the Algorithm
In this section we describe the results of applying the algorithm reported in this work to the datasets described in section 5.1.
6.1 Independent Performance Assessments
It may of interest to see how the estimates from Gaussian Progress regression compare to estimates from other methods for large regression tasks. We report the performance of two methods. XGBoost (chen2016xgboost ()) is a recent tree based algorithm using the gradient boosting framework that is very scalable and gaining adoption among practitioners. Trees partition the predictor space and can account for interaction. This method uses boosting instead of bagging. Difference in accuracy estimates between XGBoost and the proposed method for a particular dataset could be attributed to the influence of boosting or effects of interaction among variables. In most cases the estimates from XGBoost were comparable to the estimates from the proposed method. As discussed in 6.4, it is possible to combine estimates from XGBoost with the estimates from the proposed method using stacking. The Generalized Additive Model (GAM) (friedman2001elements ()) is a scalable regression technique. As the name suggests, it fits an additive model in terms of smooth nonparametric functions (typically splines) of the predictor variables. The GAM estimate serves as an independent performance estimate from another nonparametric regression algorithm based on an additive model. In most cases the accuracy obtained from GAM’s were similar to those obtained from Gaussian Process regression.
6.2 Feature Relavance
As discussed in section 3, an additive model worked well for the datasets used in this study. Further, in all these datasets, the response depended only a small number of attributes. XGBoost can report feature importance. This summarized in Table 1 below. The response depends on a small subset of predictors in all the datasets reported in this study. GPy, (gpy2014 ()) the package used to implement the experiments reported in this work, implements the Automatic Relevance Determination (see (Rasmussen06gaussianprocesses, , Chapter 5)) feature. Features that are not relevant are dropped from the model.
Dataset  Features  Important Features  

1  Airline Delay  11  1 
2  Ailerons  40  3 
3  Power Plant  4  2 
4  Delta Elevators  6  2 
5  California Housing  8  3 
6  House Hold Power  12 
4 
6.3 Accuracy
As discussed in section 4, the choice of a GP method to use for a regression task depends on the problem context. The algorithm presented in this work does not require the hyperparameter values to be specified. In terms of the inputs to the algorithm, the algorithm presented in this work is similar to the Sparse Gaussian Process and the Stochastic Variational Gaussian Process. These algorithms only require the dataset, kernel and the subset size as input. We report the performance of the Sparse Gaussian Process, Stochastic Variational Gaussian Process and the proposed algorithm for each dataset. These algorithms are good choices to begin the knowledge discovery process in large datasets. The results of applying Algorithm 1 to all the datasets listed in section 5.1 are shown in Table 2.
Dataset  BM1  BM2  POE  SVGP  SPGP  XGBoost  GAM  SD 

Ailerons  0.000214  0.000208  0.000220  0.000220  0.000220  0.000220  0.000200  0.000410 
Delta Elevators  0.001510  0.001450  0.001547  0.001460  0.001460  0.001432  0.001444  0.002370 
CCPP  4.32  4.24  4.27  5.19  4.10  3.73  4.11  17.07 
Cal Housing  0.293  0.294  0.293  NA  NA  0.240  0.281  0.569 
Airline Delay  8.75  8.75  8.74  8.85  10.94  8.74  9.45  31.31 
HPC  2.039  1.992  2.120  NA  NA  1.63  2.18  3.24 
Sinc Function  0.0267  0.0200  0.0132  0.0627  0.0170  NA  NA  0.0634 
For each dataset, Table 2 includes the following columns:

POE: This is the result of applying the Product of Experts Algorithm (Equation 11) instead of bagging. The subset size used is the same as that for Bagging (method 1).

SVGP: This is the result of applying the SVIGP algorithm.

SPGP: This is the result of applying the Sparse GP algorithm.

XGBoost: This is the result obtained from the XGBoost algorithm.

GAM: This is the result obtained from GAM.

SD: This is the standard deviation of the response. If we use a model that simply predicts the mean response for each test point, then the standard deviation represents the error associated with such a model. For regression models to be useful, they must perform better than this model.
For all datasets the split between the training and the test sets was 70%  30%. The reported accuracies are on the test set. The subset sizes required for each dataset is captured by the parameter. The values associated with Table 2 are presented in Table 3. The values for the bagging methods and POE correspond the most accurate estimates we could obtain for the datasets on a laptop with 16 GB of RAM. The SVGP and the SPGP columns of Table 3, represent the at which the best performance was obtained for the Stochastic Variational GP and the Sparse GP methods. For the airline dataset, the best performance for Stochastic Variational GP was obtained at a smaller subset size than the one used for the proposed method. For the airline dataset, Sparse GP showed no improvement in performance for greater than
Dataset  BM1  BM2  POE  SPGP  SVGP  

Ailerons  0.30  0.30  0.30  0.30  0.30  
2  Delta Elevators  0.30  0.30  0.30  0.30  0.30 
3  CCPP  0.6  0.6  0.6  0.6  0.6 
4  Cal Housing  0.625  0.625  0.625  NA  NA 
5  Airline Delay  0.56  0.49  0.0.56  0.33  0.523 
6  HPC  0.43  0.38  0.43  NA  NA 
7  Sinc Function  0.50  0.50  0.50  NA  NA 
The kernels used for the datasets are shown in Table 4.
Dataset  Kernel  

1  Ailerons  sum of linear and RBF 
2  Delta Elevators  sum of linear and RBF 
3  CCPP  sum of linear and RBF 
4  California Housing  sum of PeriodicMatern32, Linear, RBF and a product kernel of Linear and RBF on the medianIncome attribute 
5  Airline Delay  sum of BF, linear and White Noise 
6  HPC  sum of Bias, Cosine, RBF, Linear and Brownian Motion 
7  Sinc Function  RBF 
The experiments reported in this work used gpy2014 () for implementation. The Sparse GP implementation in this package is based on the Variational Sparse GP method (see titsias2009variational ()). The Stochastic Variational GP implementation is the one used in hensman2013gaussian (). Both these methods use Stochastic Gradient Descent for hyperparameter learning. Kernel support for these methods is also limited to simple kernels and kernels like the Brownian Motion Kernel or the Periodic Matern kernel are not supported. For this reason, we do not report accuracy for SVIGP or the Sparse GP for datasets that required these complex kernels (Household Power Consumption and California Housing ). This also highlights the fact that Algorithm 1 may be a good candidate for datasets with the desired characteristics described earlier but requiring a complex kernel to model the underlying function.
An analysis of Table 2 shows that the proposed method performs as well as XGBoost or GAM’s in most cases. In some cases, XGBoost does marginally better. In these cases we explored if it is possible to construct a better estimator using both these estimators. This is discussed next.
6.4 Combining Estimators
Combining estimators is not a new idea wolpert1992stacked (). More recently, lloyd2014gefcom2012 () has combined gradient boosted tree models with Gaussian Process models. This prompted us to explore combining estimators for the data sets where the Gaussian Process model produced a slightly lower accuracy. There were three datasets where the performance of the Gaussian Process model was slightly lower than the XGBoost model. These were the Combined Cycle Power Plant, California Housing and the Household Power Consumption dataset. To combine estimators, the output of the XGBoost models and the GP models were used as the inputs to a classifier (the stacking classifier). The response from the classifier was the best model for a given set of XGBoost and GP model responses. A Knearest neighbor classifier was used for this purpose. The best value of K was determined through experimentation. Given a test point, the estimates from the XGBoost model and the GP model can be obtained. The classifier then predicts best model for this set of estimates. This model is then used to provide the estimate for the test point. This procedure improved the accuracy for the California Housing and Combined Cycle Power Plant datasets. A Knearest neighbor regression performed better than the Knearest neighbor classifier for the Household Power Consumption dataset. The results are shown in Table 5
Dataset  RMSE  

1  California Housing  0.211 
2  CCPP  2.943 
3  HPC  1.610 
A comparison of Table 5 with Table 2 shows that combining estimators yields solutions that are more accurate. The idea of combining estimators can be refined in many ways. We have not included the Sparse GP and Stochastic Variational GP in this solution. The choice of the classifier or regression model to combine estimates from the component models is another modeling decision. The intent here is to show that it is possible to combine the GP model developed using the method presented in this work with other regression models to produce solutions that are better than the individual solutions. An optimal choice of the estimators and the method used for stacking is beyond the scope of this work.
7 Conclusion
There are many methods to scale Gaussian Process regression to large datasets. The appropriate choice depends on the problem context. The proposed method is appropriate for large datasets with a small set of important features for which the kernel characteristics are unknown. Kernel choices can be determined through exploratory data analysis duvenaud2014automatic (). Kernel hyperparameters can be learned from the data. The data for the experiments reported in this work came from diverse application areas and in a wide range of sizes (few thousand to two million). In these datasets, the target variable depended on a small set of variables and an additive model matched the performance of models that permit interaction like XGBoost. This suggests that datasets with these characteristics are not uncommon. Results from Minimax theory for nonparametric regression indicate that additive models that depend on a small set of predictors have an attractive rate of convergence. This suggests that we can develop adequate regression models with a small subset of samples from the dataset. This was consistent with results observed in the experiments conducted as part of this study. The Stochastic Variational Gaussian Process and the Sparse Gaussian Process are also good candidates for problems with these characteristics. The results of this study show that the proposed algorithm can match or exceed the performance of the Sparse Gaussian Process or the Stochastic Variational Gaussian Process. The results of this study also show that Gaussian Processes can be as effective as ensemble methods like Gradient Boosted Trees on large datasets. Gaussian Processes are based on a probabilistic frame work and can provide uncertainty estimates directly as compared to other tools for large regression tasks like XGBoost. This could very important for some applications. For example, an analyst may be interested in the probability of a particular amount of delay given information for a particular flight. The regression function is also interpretable in the case of Gaussian Process models in contrast to methods like gradient boosted trees. Therefore Gaussian Process models can be good explanatory models as well as good predictive models. An important feature of this algorithm is the simplicity of implementation. In Internet applications that process continuous streams of data, frequent model development and deployment is needed. An algorithm that is simple but effective may fit these applications well. Finally it should be noted that it is possible to combine this algorithm with other algorithms like Gradient Boosted Trees using model stacking to achieve performance gains.
8 References
Footnotes
 journal: Journal of LaTeX Templates
 timestamp is counted as a feature. Variables derived from timestamp are not included in this count.
References
 C. E. Rasmussen, Gaussian processes for machine learning, MIT Press, 2006.

Z. Ghahramani,
A tutorial on
gaussian processes (or why i don’t use svms), mLSS Workshop talk by Zoubin
Ghahramani on Gaussian Processes [Accessed: 2016 07 19] (2011).
URL http://mlss2011.comp.nus.edu.sg/uploads/Site/lect1gp.pdf  C. J. Stone, Additive regression and other nonparametric models, The annals of Statistics (1985) 689–705.
 M. K. Titsias, Variational learning of inducing variables in sparse gaussian processes., in: AISTATS, Vol. 12, 2009, pp. 567–574.
 J. Hensman, N. Fusi, N. D. Lawrence, Gaussian processes for big data, in: Conference on Uncertainty in Artificial Intellegence, auai.org, 2013, pp. 282–290.
 T. Chen, C. Guestrin, Xgboost: A scalable tree boosting system, arXiv preprint arXiv:1603.02754.
 J. Friedman, T. Hastie, R. Tibshirani, The Elements of Statistical Learning, Vol. 1, Springer series in statistics Springer, Berlin, 2001.
 D. H. Wolpert, Stacked generalization, Neural networks 5 (2) (1992) 241–259.
 J. R. Lloyd, Gefcom2012 hierarchical load forecasting: Gradient boosting machines and gaussian processes, International Journal of Forecasting 30 (2) (2014) 369–374.
 D. Duvenaud, Automatic model construction with gaussian processes, Ph.D. thesis, University of Cambridge (2014).
 C. J. Stone, Optimal global rates of convergence for nonparametric regression, The annals of statistics (1982) 1040–1053.
 E. Snelson, Z. Ghahramani, Sparse gaussian processes using pseudoinputs, in: Advances in neural information processing systems, 2005, pp. 1257–1264.
 C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning), The MIT Press, 2005.
 C. M. Bishop, Pattern recognition, Vol. 128, 2006.
 L. Breiman, Random forests, Machine learning 45 (1) (2001) 5–32.
 Y. Cao, D. J. Fleet, Generalized product of experts for automatic and principled fusion of gaussian process predictions, arXiv preprint arXiv:1410.7827.
 A. B. Tsybakov, Introduction to nonparametric estimation. revised and extended from the 2004 french original. translated by vladimir zaiats (2009).
 L. Györfi, M. Kohler, A. Krzyzak, H. Walk, A distributionfree theory of nonparametric regression, Springer Science & Business Media, 2006.
 Y. Yang, S. T. Tokdar, et al., Minimaxoptimal nonparametric regression in high dimensions, The Annals of Statistics 43 (2) (2015) 652–674.
 Y. Yang, D. B. Dunson, et al., Bayesian manifold regression, The Annals of Statistics 44 (2) (2016) 876–905.
 A. Zaytsev, E. Burnaev, Minimax approach to variable fidelity data interpolation, in: Artificial Intelligence and Statistics, 2017, pp. 652–661.
 J. QuiñoneroCandela, C. E. Rasmussen, C. K. Williams, Approximation methods for gaussian process regression, Largescale kernel machines (2007) 203–224.
 A. G. Wilson, Covariance kernels for fast automatic pattern discovery and extrapolation with gaussian processes, Ph.D. thesis, University of Cambridge (2014).
 A. Wilson, R. Adams, Gaussian process kernels for pattern discovery and extrapolation, in: Proceedings of the 30th International Conference on Machine Learning (ICML13), 2013, pp. 1067–1075.
 P. Drineas, M. W. Mahoney, On the nyström method for approximating a gram matrix for improved kernelbased learning, journal of machine learning research 6 (Dec) (2005) 2153–2175.
 A. Rahimi, B. Recht, Random features for largescale kernel machines, in: Advances in neural information processing systems, 2008, pp. 1177–1184.
 F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, Scikitlearn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
 GPy, GPy: A gaussian process framework in python, http://github.com/SheffieldML/GPy (2012–2014).
 N. Srivastava, G. E. Hinton, A. Krizhevsky, I. Sutskever, R. Salakhutdinov, Dropout: a simple way to prevent neural networks from overfitting., Journal of Machine Learning Research 15 (1) (2014) 1929–1958.
 M. Belyaev, E. Burnaev, Y. Kapushev, Computationally efficient algorithm for gaussian process regression in case of structured samples, Computational Mathematics and Mathematical Physics 56 (4) (2016) 499–513.
 V. Tresp, A bayesian committee machine, Neural Computation 12 (11) (2000) 2719–2741.
 M. P. Deisenroth, J. W. Ng, Distributed gaussian processes, in: International Conference on Machine Learning (ICML), Vol. 2, 2015, p. 5.
 S. Park, S. Choi, Hierarchical gaussian process regression., in: ACML, 2010, pp. 95–110.
 R. B. Gramacy, D. W. Apley, Local gaussian process approximation for large computer experiments, Journal of Computational and Graphical Statistics 24 (2) (2015) 561–578.
 R. B. Gramacy, laGP: Largescale spatial modeling via local approximate gaussian processes in R, Journal of Statistical Software 72 (1) (2016) 1–46. doi:10.18637/jss.v072.i01.

M. Lichman, UCI machine learning
repository (2016).
URL http://archive.ics.uci.edu/ml 
L. Targo.,
Large
regression datasets (2016).
URL http://www.transtats.bts.gov/DL_SelectFields.asp?Table_ID=236 
B. USDOT,
Rita
airline delay data download (2016).
URL http://www.transtats.bts.gov/DL_SelectFields.asp?Table_ID=236