# gaussian process regression for large datasets

I've been learning about Gaussian process regression from online videos and lecture notes, my understanding of it is that if we have a dataset with $n$ points then we assume the data is sampled from an $n$-dimensional multivariate Gaussian. So my question is in the case where $n$ is 10's of millions does Gaussian process regression still work? Will the kernel matrix not be huge rendering the process completely inefficient? If so are there techniques in place to deal with this, like sampling from the data set repeatedly many times? What are some good methods for dealing with such cases?

• Why do you want to use Gaussian process and not something that is destined for dealing with large data?
– Tim
Commented Feb 2, 2018 at 17:08

in the case where 𝑛 is 10's of millions does Gaussian process regression still work?

Not in the standard sense of constructing and inverting a large matrix. You have two options: 1)choose a different model or 2) make an approximation.

1) Some GP-based models can be scaled to very large data sets, such as the Bayesian committee machine linked in the answer above. I find this approach rather unsatisfactory though: there are good reasons for choosing a GP model, and if we are to switch to a more computable model we might not retain the properties of the original model. The predictive variances of the BCM depend strongly on the data split, for example.

2) The 'classical' approach to approximation in GPs is to approximate the kernel matrix. There's a good review of these sorts of methods here: http://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf. In fact, we can usually see these matrix approximations as model approximations, and lump them in with the Bayesian committee machine: they're changes to the model and it can be hard to understand when those changes might be pathological. Here's a super review: https://papers.nips.cc/paper/6477-understanding-probabilistic-sparse-gaussian-process-approximations.pdf

The way I advocate for making approximations for large GPs is to avoid approximating the kernel matrix or the model, and approximate the posterior distribution using variational inference. A lot of the computations look like a 'low rank' matrix approximation, but there is one very desirable property: the more computation you use (the more "ranks") the close the approximation is to the true posterior, as measured by the KL divergence.

These articles are a good starting point: http://proceedings.mlr.press/v5/titsias09a/titsias09a.pdf https://arxiv.org/pdf/1309.6835

I wrote a longer article on the same argument here: https://www.prowler.io/blog/sparse-gps-approximate-the-posterior-not-the-model

In practice, the variational approximation works really well in a lot of cases. I've used it extensively in real applications. And more recently there's been some excellent theory to back up why it should work (https://arxiv.org/abs/1903.03571).

A final plug: variational inference in GPs is implemented in gpflow (https://github.com/GPflow/GPflow)

• Just for the record if someone comes across this and doesn’t know, James would be one of the authorities in the field along with Mike Osborne, Neill Lawrence and so on - kind of cool he’s answering Q&A on SE
– j__
Commented Feb 20, 2020 at 14:11

There are a wide range of approaches to scale GPs to large datasets, for example:

Low Rank Approaches: these endeavoring to create a low rank approximation to the covariance matrix. Most famously perhaps is Nystroms method which projects the data onto a subset of points. Building on from that FITC and PITC were developed which use pseudo points rather than points observed. These are included in for example the GPy python library. Other approaches include random Fourier features.

H-matrices: these use hierarchical structuring of the covariance matrix and apply low rank approximations to each structures submatrix. This is less commonly implemented in popular libraries.

Kronecker Methods: these use Kronecker products of covariance matrices in order to speed up the computational over head bottleneck.

Bayesian Committee Machines: This involves splitting your data into subsets and modeling each one with a GP. Then you can combine the predictions using the optimal Bayesian combination of the outputs. This is quite easy to implement yourself and is fast but kind of breaks your kernel is you care about that - Mark Deisenroth’s paper should be easy enough to follow here.

Usually, what you can do is to train Gaussian Processes on subsamples of your dataset (bagging). Bagging is implemented in sk learn and can be used easily. See per example the documentation.

Calling $n$ the number of observations, $n_{bags}$ the number of bags you use and $n_{p}$ the number of points per bag, this allow to change training time from a $O(n^3)$ to a $O(n_{bags}n_{p}^3)$. Therefore, with small bags but using all the data, you can achieve a much lower training time. Unfortunately, this often reduces the performance of the model.

Apart from bagging techniques, there is some active research about making the Gaussian Process Regressions scalable. The article Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP) proposes to reduce the training time to a $O(n)$ and comes with a matlab code.

You may wish to take a look at the spatial statistics literature where different classes of massively scalable Gaussian processes have witnessed significant development with the number of space or space-time coordinates (inputs into the Gaussian processes) into several millions. The idea is not to rely solely upon computations but build different classes of Gaussian processes that can scale up to tens of millions of locations. Most of these methods are agnostic to any specific algorithm and can be adapted to full inference using MCMC or faster approximate Bayesian algorithms such as Variational Bayes or INLA.

Some recent pointers are available here:

Sparsity-inducing Gaussian processes including Nearest Neighbor Gaussian Processes (NNGPs) based upon an approximation due to Vecchia (arranged chronologically):

Vecchia's original paper https://doi.org/10.1111%2Fj.2517-6161.1988.tb01729.x

While the paper was published in 1988, interest in this approach was rather lukewarm until the following NNGP paper was published in 2016:

For dynamic NNGP with space-time inputs (where neighbors of inputs are estimated rather than fixed): https://doi.org/10.1214/16-AOAS931

https://doi.org/10.1080/10618600.2018.1537924 (focusing on efficient Bayesian algorithms for NNGP models)

https://doi.org/10.1214%2F19-STS755 (accompanying R package: https://cran.r-project.org/web/packages/GPvecchia/index.html)

Meshed Gaussian processes https://doi.org/10.1080/01621459.2020.1833889 (accompanying R package: https://cran.r-project.org/web/packages/meshed/index.html)

Apart from this, there is ExaGeoStat here:

https://cemse.kaust.edu.sa/stsds/exageostat

Here is some recent work on Variational NNGPs:

https://proceedings.mlr.press/v162/wu22h/wu22h.pdf

In terms of pure scalability, these approaches have been applied to data sets with number of inputs exceeding millions.

Review articles include:

https://doi.org/10.1214/17-BA1056R (only Bayesian)