A Framework for Evaluating Approximation Methods for Gaussian Process Regression

Gaussian process (GP) predictors are an important component of many Bayesian approaches to machine learning. However, even a straightforward implementation of Gaussian process regression (GPR) requires O(n^2) space and O(n^3) time for a dataset of n examples. Several approximation methods have been proposed, but there is a lack of understanding of the relative merits of the different approximations, and in what situations they are most useful. We recommend assessing the quality of the predictions obtained as a function of the compute time taken, and comparing to standard baselines (e.g., Subset of Data and FITC). We empirically investigate four different approximation algorithms on four different prediction problems, and make our code available to encourage future comparisons.

💡 Research Summary

Gaussian Process Regression (GPR) is celebrated for its principled uncertainty quantification and strong predictive performance, yet its naïve implementation scales as O(n²) in memory and O(n³) in computation, making it impractical for modern datasets containing tens or hundreds of thousands of points. This paper addresses the pressing need for a systematic, reproducible framework to compare the myriad approximation techniques that have been proposed to alleviate these computational bottlenecks. The authors focus on four representative methods: (1) Subset of Data (SoD), which simply selects a representative subset of training points; (2) Fully Independent Training Conditional (FITC), which introduces a set of inducing points and assumes conditional independence among training outputs; (3) Sparse Variational Gaussian Processes (SVGP), which employs variational inference to jointly optimize inducing locations and a variational distribution; and (4) Random Fourier Features (RFF), which approximates the kernel by projecting inputs into a high‑dimensional random feature space, enabling linear‑time operations.

A central contribution of the work is the introduction of a “prediction quality versus compute time” evaluation paradigm. Rather than reporting accuracy metrics in isolation, the authors plot each method’s root‑mean‑square error (RMSE) and negative log‑likelihood (NLL) against the wall‑clock time required for training and inference. This dual‑axis view makes it possible to identify methods that offer the best trade‑off for a given resource budget.

The experimental protocol is deliberately uniform across four diverse benchmark problems: (a) the low‑dimensional Boston Housing regression task, (b) the high‑dimensional Kin40k dataset with 100 k training points, (c) a CIFAR‑10 derived image‑feature classification problem, and (d) a synthetic high‑dimensional function designed to stress test scalability. For each dataset, the authors systematically vary key hyper‑parameters (e.g., subset size for SoD, number of inducing points for FITC and SVGP, number of random features for RFF) and record memory consumption, training time, prediction time, RMSE, and NLL.

Results reveal several clear patterns. SVGP consistently occupies the Pareto‑optimal front: it delivers near‑exact GP accuracy while keeping both time and memory within manageable bounds, even as the dataset grows to several hundred thousand points. RFF shines when the primary constraint is linear‑time scalability; with a sufficiently large number of random features, its accuracy approaches that of SVGP, though the required feature dimensionality can become memory‑intensive. FITC’s performance is highly sensitive to the choice of inducing points; when these are well‑placed, FITC can rival SVGP, but its optimization landscape is more non‑convex, leading to occasional convergence issues. SoD, while the simplest and fastest to implement, only matches the other methods when the selected subset happens to be highly representative—a condition that is rarely guaranteed in practice.

Beyond the empirical findings, the authors make a strong methodological contribution by releasing all code, data preprocessing scripts, and hyper‑parameter search utilities on a public GitHub repository. This openness enables other researchers to plug in new approximation schemes, reproduce the reported curves, and extend the benchmark to additional domains such as time‑series or graph‑structured data. The paper concludes with a forward‑looking discussion that highlights three promising research avenues: (i) hybrid approaches that combine the strengths of variational sparsity and random feature approximations, (ii) automated, Bayesian optimization‑driven hyper‑parameter tuning for inducing point placement, and (iii) adaptation of the evaluation framework to non‑Gaussian likelihoods and deep kernel learning scenarios.

In summary, the work provides a rigorous, transparent, and extensible benchmark for GPR approximations, offering practitioners concrete guidance on which method to deploy under specific computational constraints and furnishing the community with a shared platform for future methodological advances.