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

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.