# 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 Feb 2 '18 at 17:08

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.

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.