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 L2-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(x)+ξ(x)
ΓkN∗TKN−1yN
σ^ξ2+σ^f2−kN∗TKN−1kN∗
NNGP
t(x)T.b+w(x)+ξ(x)
t∗T.b^+ΓkN∗TKN−1(yN−TN.b^)
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 Γ making them unbiased when m is fixed and in the limit n→∞ , 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).
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) .
Figure 4Table 5
Figure 5: Risk landscape (shifted) as a function of the hyperparameters and training set size. NNGP regression with PX uniform on dX=2 disk and ν=1/2 . The parameter b^ is chosen as b^=(b^,b^) . Note the extreme flatness of b^ -landscape.
Figure 7
n
104
105
106
before
after
before
after
before
after
CAL
4.90±0.01
1.036±0.002
9.84±0.01
1.025±0.001
20.17±0.02
1.040±0.001
NLL
1.494±0.004
0.355±0.006
3.563±0.003
0.2989±0.0006
8.34±0.01
0.2755±0.0003
Table 6: Empirical CAL and 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(x)+ξ(x)
ΓkN∗TKN−1yN
σ^ξ2+σ^f2−kN∗TKN−1kN∗
NNGP
t(x)T.b+w(x)+ξ(x)
t∗T.b^+ΓkN∗TKN−1(yN−TN.b^)
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 Γ=mσ^f2σ^ξ2+mσ^f2 making them unbiased when m is fixed and in the limit n→∞ , as opposed to standard formulas used in the literature.
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.
In this work, we study scaling limits of shallow Bayesian neural networks (BNNs) via their connection to Gaussian processes (GPs), with an emphasis on statistical modeling, identifiability, and scalable inference. We first establish a general convergence result from BNNs to GPs by relaxing assumptions used in prior formulations, and we compare alternative parameterizations of the limiting GP model. Building on this theory, we propose a new covariance function defined as a convex mixture of components induced by four widely used activation functions, and we characterize key properties including positive definiteness and both strict and practical identifiability under different input designs. For computation, we develop a scalable maximum a posterior (MAP) training and prediction procedure using a Nyström approximation, and we show how the Nyström rank and anchor selection control the cost-accuracy trade-off. Experiments on controlled simulations and real-world tabular datasets demonstrate stable hyperparameter estimates and competitive predictive performance at realistic computational cost.
Gracielle Antunes de Araújo, Flávio B. Gonçalves
Universidade do Estado de Minas Gerais (UEMG), FaEnge – Faculdade de Engenharia. · Universidade Federal de Minas Gerais (UFMG), ICEx – Instituto de Ciências Exatas.
Gaussian process (GP) regression with a single global GP (GP-glo) incurs cubic computational cost, limiting scalability to large datasets. Product-of-experts GP models (GP-pro), which combine local GP models to capture global correlations, alleviate this computational burden. However, training local experts on disjoint data subsets can lead to overestimated posterior variances. We propose GP-pro-c, a product-of-experts GP model that calibrates these variances using an information-based method. The method exploits the monotonicity and submodularity of information gain in GPs to define a calibration ratio that reduces the posterior variance of individual local GP models. We evaluate GP-pro-c using negative log-likelihood (NLL), root mean squared error (RMSE), and expected normalised calibration error (ENCE). Experiments on four synthetic functions and six regression datasets show that GP-pro-c achieves average reductions of 2.3% in NLL and 12.0% in ENCE compared with the uncalibrated GP-pro model. The proposed method mitigates posterior variance overestimation while maintaining predictive accuracy and reducing computational complexity. GP-pro-c provides a promising approach for uncertainty estimation in scalable GP models and may serve as a useful surrogate model for Bayesian optimisation with high-dimensional and large-scale data.
Yean Hoon Ong, Paolo Barucca, Wei Pan +1
University College London, London, UK · The University of Manchester, Manchester, UK