The Theory and Practice of Highly Scalable Gaussian Process Regression with Nearest Neighbours

Robert Allison; Tomasz Maciazek; Anthony Stephenson

stat.MLApr 8, 2026

The Theory and Practice of Highly Scalable Gaussian Process Regression with Nearest Neighbours

Authors: Robert Allison, Tomasz Maciazek, Anthony Stephenson

Organizations: School of Mathematics University of Bristol Bristol, BS8 1UG, UK

Abstract

Gaussian process ( $GP$ ) regression is a widely used non-parametric modeling tool, but its cubic complexity in the training size limits its use on massive data sets. A practical remedy is to predict using only the nearest neighbours of each test point, as in Nearest Neighbour Gaussian Process ( $NNGP$ ) regression for geospatial problems and the related scalable $GPnn$ method for more general machine-learning applications. Despite their strong empirical performance, the large- $n$ theory of $NNGP/GPnn$ remains incomplete. We develop a theoretical framework for $NNGP$ and $GPnn$ regression. Under mild regularity assumptions, we derive almost sure pointwise limits for three key predictive criteria: mean squared error ( $MSE$ ), calibration coefficient ( $CAL$ ), and negative log-likelihood ( $NLL$ ). We then study the $L_2$ -risk, prove universal consistency, and show that the risk attains Stone's minimax rate $n^{-2α/(2p+d)}$ , where $α$ and $p$ capture regularity of the regression problem. We also prove uniform convergence of $MSE$ over compact hyper-parameter sets and show that its derivatives with respect to lengthscale, kernel scale, and noise variance vanish asymptotically, with explicit rates. This explains the observed robustness of $GPnn$ to hyper-parameter tuning. These results provide a rigorous statistical foundation for $NNGP/GPnn$ as a highly scalable and principled alternative to full $GP$ models.

Figures & tables

	Response Model	Predictive Mean	Predictive Variance
$GPnn$	$f(\boldsymbol{x})+\xi(\boldsymbol{x})$	$\Gamma\,{\boldsymbol{k}_{\mathcal{N}}^{*}}^{T}\,K_{\mathcal{N}}^{-1}\,\boldsymbol{y}_{\mathcal{N}}$	$\hat{\sigma}_{\xi}^{2}+\hat{\sigma}_{f}^{2}-{\boldsymbol{k}_{\mathcal{N}}^{}}^{T}\,K_{\mathcal{N}}^{-1}\boldsymbol{k}_{\mathcal{N}}^{}$
$NNGP$	$\boldsymbol{t}(\boldsymbol{x})^{T}.\boldsymbol{b}+w(\boldsymbol{x})+\xi(\boldsymbol{x})$	$\boldsymbol{t}_{}^{T}.\hat{\boldsymbol{b}}+\Gamma\,{{k}_{\mathcal{N}}^{}}^{T}\,K_{\mathcal{N}}^{-1}\left(\boldsymbol{y}_{\mathcal{N}}-T_{\mathcal{N}}.\hat{\boldsymbol{b}}\right)$

Table 1: Summary of response models, predictive mean and variance in

GPnn

and

NNGP

. See the main text for the explanation of the symbols. As explained in the main text, for

NNGP

these formulas apply only in the conditional Gaussian predictive distribution which fixes the hyperparameters. The predictive means are corrected by the coefficient

\Gamma

making them unbiased when

m

is fixed and in the limit

n\to\infty

, as opposed to standard formulas used in the literature.

Figure 1: The

GPnn

workflow, including the above-described cheap hyper-parameter estimation and the calibration step of the predictive variance (see the main text for explanation).

Panel of Figure 2 — Figure 2: Experiment results on a suite of UCI datasets. Optimal calibration performance is 1 (indicated by a black dashed line). Lower is better for NLL and RMSE. Y-axis truncated for readability for Calibration due to large values on “Ctslice”. “symlog” is logarithmic axis rescaling applied to the $y$ -axis on both positive and negative values (“symmetric”). Adapted from Allison et al. (2023) .

Table 5

Figure 5: Risk landscape (shifted) as a function of the hyperparameters and training set size.

NNGP

regression with

P_{\mathcal{X}}

uniform on

d_{\mathcal{X}}=2

disk and

\nu=1/2

. The parameter

\hat{\boldsymbol{b}}

is chosen as

\hat{\boldsymbol{b}}=(\hat{b},\hat{b})

. Note the extreme flatness of

\hat{\boldsymbol{b}}

-landscape.

$n$	$10^{4}$		$10^{5}$		$10^{6}$
	before	after	before	after	before	after
$\widehat{CAL}$	$4.90\pm 0.01$	$1.036\pm 0.002$	$9.84\pm 0.01$	$1.025\pm 0.001$	$20.17\pm 0.02$	$1.040\pm 0.001$
$\widehat{NLL}$	$1.494\pm 0.004$	$0.355\pm 0.006$	$3.563\pm 0.003$	$0.2989\pm 0.0006$	$8.34\pm 0.01$	$0.2755\pm 0.0003$

Table 6: Empirical

\widehat{CAL}

and

\widehat{NLL}

in the synthetic-data

NNGP

experiment before and after the post-hoc calibration. Results show excellent effectiveness of the calibration procedure, even in a strongly mismatched hyperparameter setting. The errors are obtained by calculating the STD over

4

independent realisations of the training set.

Appendix figures & tables1 asset

Supplementary material from the paper’s appendix.

Appendix

	Response Model	Predictive Mean	Predictive Variance
$GPnn$	$f(\boldsymbol{x})+\xi(\boldsymbol{x})$	$\Gamma\,{\boldsymbol{k}_{\mathcal{N}}^{*}}^{T}\,K_{\mathcal{N}}^{-1}\,\boldsymbol{y}_{\mathcal{N}}$	$\hat{\sigma}_{\xi}^{2}+\hat{\sigma}_{f}^{2}-{\boldsymbol{k}_{\mathcal{N}}^{}}^{T}\,K_{\mathcal{N}}^{-1}\boldsymbol{k}_{\mathcal{N}}^{}$
$NNGP$	$\boldsymbol{t}(\boldsymbol{x})^{T}.\boldsymbol{b}+w(\boldsymbol{x})+\xi(\boldsymbol{x})$	$\boldsymbol{t}_{}^{T}.\hat{\boldsymbol{b}}+\Gamma\,{{k}_{\mathcal{N}}^{}}^{T}\,K_{\mathcal{N}}^{-1}\left(\boldsymbol{y}_{\mathcal{N}}-T_{\mathcal{N}}.\hat{\boldsymbol{b}}\right)$

Appendix

Table 7: Summary of response models, predictive mean and variance in

GPnn

and

NNGP

. See the main text for the explanation of the symbols. The predictive means are corrected by the coefficient

\Gamma=\frac{\hat{\sigma}_{\xi}^{2}+m\hat{\sigma}_{f}^{2}}{m\hat{\sigma}_{f}^{2}}

making them unbiased when

m

is fixed and in the limit

n\to\infty

, as opposed to standard formulas used in the literature.

Explore similar work

Scalable Gaussian process inference via neural feature maps
May 11, 2026Anthony StephensonGaussian Processes Kernel Methods
School of Mathematics University of Bristol
From Shallow Bayesian Neural Networks to Gaussian Processes: General Convergence, Identifiability and Scalable Inference
Feb 26, 2026Gracielle Antunes de Araújo, Flávio B. GonçalvesBayesian Neural Networks Shallow Neural Networks
Universidade do Estado de Minas Gerais (UEMG), FaEnge – Faculdade de Engenharia. · Universidade Federal de Minas Gerais (UFMG), ICEx – Instituto de Ciências Exatas.
Information-Based Calibration of Uncertainty Quantification in Product-of-Experts Gaussian Process Models
Aug 29, 2026Yean Hoon Ong, Paolo Barucca, Wei Pan +1Uncertainty Quantification Gaussian Processes
University College London, London, UK · The University of Manchester, Manchester, UK

May 11, 2026J/K move · Enter open · S save

Scalable Gaussian process inference via neural feature maps

Anthony Stephenson

School of Mathematics University of Bristol

Gaussian Processes Kernel Methods Gaussian Process Regression

Open paper AskPDF arXiv

We present a theoretically grounded Gaussian process framework that leverages neural feature maps to construct expressive kernels. We show that the learned feature map can be interpreted as an optimal low-rank approximation to a Gram matrix derived from an implied RKHS, from which we establish consistency of the GP posterior. We further analyse the spectral properties of the induced kernels and introduce product feature-map kernels to address oversmoothing. This simple yet powerful approach enables fast, scalable, and accurate exact GP inference with minimal upfront work. The flexibility of kernel design supports seamless application to both regression and classification tasks across diverse data modalities, including tabular inputs and structured domains such as images. On benchmark datasets, this approach surpasses pre-existing methods in terms of accuracy and training and prediction efficiency.

Abstract

Explore similar work

Scalable Gaussian process inference via neural feature maps

From Shallow Bayesian Neural Networks to Gaussian Processes: General Convergence, Identifiability and Scalable Inference

Information-Based Calibration of Uncertainty Quantification in Product-of-Experts Gaussian Process Models