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

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.

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__ Feb 20 at 14:11

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.