Leveraging Locality and Robustness to Achieve Massively Scalable Gaussian Process Regression
Authors: Robert F Allison
Anthony Stephenson
Samuel F
Edward Pyzer-Knapp
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian Processes
Bayesian Inference
Regression
Bayesian Nonparametrics
Kernel Methods

Abstract: The accurate predictions and principled uncertainty measures provided by GP regression incur $O(n^3)$ cost which is prohibitive for modern-day large-scale applications. This has motivated extensive work on computationally efficient approximations. We introduce a new perspective by exploring robustness properties and limiting behaviour of GP nearest-neighbour (GPnn) prediction. We demonstrate through theory and simulation that as the data-size $n$ increases, accuracy of estimated parameters and GP model assumptions become increasingly irrelevant to GPnn predictive accuracy. Consequently, it is sufficient to spend small amounts of work on parameter estimation in order to achieve high MSE accuracy, even in the presence of gross misspecification. In contrast, as $n \rightarrow \infty$, uncertainty calibration and NLL are shown to remain sensitive to just one parameter, the additive noise-variance; but we show that this source of inaccuracy can be corrected for, thereby achieving both well-calibrated uncertainty measures and accurate predictions at remarkably low computational cost. We exhibit a very simple GPnn regression algorithm with stand-out performance compared to other state-of-the-art GP approximations as measured on large UCI datasets. It operates at a small fraction of those other methods' training costs, for example on a basic laptop taking about 30 seconds to train on a dataset of size $n = 1.6 \times 10^6$.

Leveraging Locality and Robustness to Achieve Massively Scalable Gaussian Process Regression
Robert Allison
Anthony Stephenson
Samuel F Alan Turing Institute
Department of Mathematics
Department of Mathematics Bristol University
Bristol University
marfa@bristol.. ac.uk
Edward Pyzer-Knapp
IBM Research
Abstract
The accurate predictions and principled uncertainty measures provided by GP regression incur O (n²)) cost which prohibitive for modern-day large-scale appli- cations. This has motivated extensive work computationally efficient approxi- mations. We introduce a new perspective by exploring robustness properties and limiting behaviour of GP nearest neighbour (GPnn) prediction. We demonstrate through theory and simulation that as the data-size n increases, accuracy esti- mated parameters and GP model assumptions become increasingly irrelevant GPnn predictive accuracy. Consequently, it is sufficient to spend small amounts in the
work parameter estimation in order achieve high MSE accuracy, presence of gross misspecification. In contrast, as n - oo, uncertainty calibration and NLL are shown remain sensitive just one parameter, the additive noise- variance; but we show that this source of inaccuracy can be corrected for, thereby achieving both well-calibrated uncertainty measures and accurate predictions remarkably low computational cost. We exhibit very simple GPnn regression algorithm with stand-out performance compared to other state-of-the- proximations as measured on large UCI datasets. It operates at a small fraction of those other methods' training costs, for example on basic laptop taking about 30 seconds to train on dataset of size 1.6 x 10".
GP ap-
1
Introduction
We first briefly review the computational cost of exact GP regression and the motivation for this paper: Given n training samples X, y. where X E Raxd has feature vector E Rdd in its i'th row and y R". exact GP regression (36) makes use of n gram matrix K Kx.e constructed from
pre-specified positive definite covariance function o(s.)) Rdd Rd R together with length-scale, additive-noise variance and kernel-scale "hyperparameters" = (1, o of). In the training phase
estimates the hyperparameters, (i, o². 0f). are obtained by minimising the loss function
loss(f) () log p( y X 0) 1 (2" + log Kol + n log(2m))
Then for subsequent predictions the predictive distribution point x" E Id defined by y" X, _Nppp,. o²²
(2)
n¹ = k-
(3)
0² 07 k. T
K K. +
(4)
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

where K = Kg with components [5] (g the vector k' has components k ka ); and (o( x') [fc(=/1,a//1 + 6 o with a "normalised" covariance function c(:,)) such that c(æ; æ) = 1. The derivation of these steps is based the assumption that the underlying random field Gaussian, as the additive noise with variance o..
The single cost-O(n²) step of inverting K' is needed repeatedly to compute the loss. Sophisticated implementations reduce this toward O (n²) ((35), but even that cost is generally impractical for > 10 For a survey of numerous GP approximations and their reduced costs see [9].
Machine learning methods must tackle massive data problems handle many modern day appli- cations. Revolutionary developments in neural network methodologies achieve this, but Bayesian predictive methodologies, particular GP regression with its major advantages of robustness and uncertainty measures, are somewhat behind the curve. This motivates development of fast implemen- tations retaining the accuracy and princirie uncertainty of exact GPs.
2 Background and Paper Outline
A feature common to all mainstream GP approximations is that training and prediction processes make joint use the same underlying mathematical constructions. In the "subset-of-data" method the same subset of data is used both for parameter estimation and prediction. Similarly, in the various Bayesian committee methods (118) hyperparameters are estimated using a collection subsets of data and then the same subsets are used in combination to make predictions. In the variational ([14 31) and other inducing point methods parameters are estimated using a low rank approximation to the kernel gram matrix and then the same low-rank matrix approximation is used make predictions. Despite being almost universally adopted there is no obvious reason why constraining algorithms to the same constructions for estimation and prediction will help rather than hinder the end goal of high performance at low cost. Whilst there are some passing mentions of decoupling prediction and estimation in the literature e.g- 28
it has not been adopted as mainstream approach.
Our first observation that allowing parameter-estimation and prediction processes to become decoupled may provide the flexibility greatly improve cost-accuracy trade-off. As shown in Figure GP approximations first obtain point estimate of the kernel hyperparameterss Ô from training data and then feed 0 into predictive process. Our end-goal is only to obtain accurate and well-calibrated predictive distributions of y" at each target point x"; obtaining accurate parameter estimates is not goal in itself. It follows that the computational budget devoted to parameter estimation need only be sufficient to provide parameters capable of delivering accurate and well- calibrated predictions.
Training data (X, y)
Training data (X, 3)
Parameter estimation
Prediction
Distribution of y" X, at xt
Target point xx
Figure 1: Flowchart of the GP regression procedure. The dashed box indicates the usual approach of combining the parameter estimation and prediction tasks under one strategy.
This need not mean that Ô is an accurate estimate of the parameters which leads us to the second observation and main theoretical component of this paper: In section 5] theory and simulations reveal that under widely applicable circumstances, as increases the mean squared error (MSE) predictive accuracy obtained from GP nearest neighbour prediction becomes increasingly insensitive model misspecification, i.c. insensitive the wrong choice covariance function, the choice of i, o3 and o2 and even insensitive to departures from Gaussian model assumptions made for the underlying stochastic process and additive noise. Similarly the negative log likelihood (NI..) predietive-accuracy becomes insensitive to all of those factors apart from the variance of the additive noise of. In 6.1

---

we describe simple calibration step that corrects for the latter inaccuracy thereby achieving near optimal limiting NLL values in addition well-calibrated uncertainty measures whilst leaving the well-behaved MSE values completely unaltered. We apply these overall observations to construct highly efficient, accurate and well calibrated regression algorithm in section 6
Our key contributions:
Demonstration of GPnn robustness against model and parameter misspec- ification through theory and simulation ( derivation of explicit formulae for the limiting MSE, NLL and calibration performance of GPnn as - xc I.T translation of this theory into a new GP approximation framework with stand-ou performance relative to other ffate-ofthhe-ee mations (67.1) a simple generic method for re- re-calibrating uncertainty measures in GP regression with immediate applications to improving calibration of other GP approximations such as SVGP [6.]); achievement of massive scalability for GPs, for example a 100x speed-up over state-of-the-art methods on a 1.6x 10th training set whilst also improving upon their performance (7.1 demonstrating that provably best possible MSE, NLL and calibration performance can be closely approached on data that is grossly misspecified relative to GP model assumptions (7.2)
GP approxi-
3 Performance Measures, Weak and Strong Calibration
Along with many other GP publications we use mean squared error MSE (or its square root RMSE) and negative log likelihood (NLL). both computed from held-aside test data, to assess predictive performance. These are simply the mean values of e = (yg 2²)² and 1 = 0.5 (log o + (,₁ ";)² /op log 2T) respectively. However, we find those measures alone inadequate for determining how well calibrated a predictive distribution is. We define "weakly calibrated" prediction to mean that 1 {("" 1² /0)²//0²) = and accordingly use "calibration" to be measure of how well the average value of 2 = (vf nu "zt)//a"" over test-data agrees with 1. This choice of metric (also made use of in [16]) can be motivated as follows: For well-calibrated GP the expected squared deviation
of S from μ4," should match the corresponding predictive variance of² ( i.e. Ey- X.y {2"} 1. Hence, observing an average of 2" values close to 1 is consistent with necessary condition for effective calibration. In practice, we find that GP approximation methods can fall well short of this condition (e.g. see the LHS plot in Figure and TTall 3) for those results in tabular form) whilst this is not evident from their MSE and NLL values alone. A better measure of calibration "sstrong-calibration would have been to see how well percentiles of the predictive distributions agree with those observed in test data, e.g. see (24] but defer such a refinement to future work.
4 Prediction Method and Sources of Misspecification
4.1 GP Nearest Neighbour Prediction
We now describe what we mean by "GP nearest neighbour (GPnn) prediction". Assume that we are given parameters θ = (i, a², 03) obtained from the parameter estimation phase of Figure 1 Then to compute the estimated pointwise-distribtton of y" at x" indicated in Figure I we find the m nearest training-set neighbours N = N ((=*) x" and apply exactly the same GP prediction formulae as in ( (3) and (4) but with X Rnxd replaced by N Rmxd and y E R" replaced by yN R Note that in this setup conditioning on N(x*) ) is equivalent to conditioning on the full input matrix X. We obtain:
5" N(x'), YN N(NNN
(5)
1/ = :x RNN
N
UN
(6)
N 2 of N
K Ki
where we have used hatted notation in a generic manner to cover all the potential sources of mis- specification (4.2) that might arise when carry out these predictions. The niv, x² parameters substituted for o²2 when computing the performance measures described in section
Note: In this paper we use Fuclidean distance for nearest neighbour assignment but more generally could employ metric defined by the covariance function see A. For the covariance functions used in this paper these metrics are "cquivalent" because one is monotonic function of the other.
3

---

In our algorithmic implementations we replace exact m nearest neighbour algorithm with much more efficient approximate nearest neighbour algorithm as discussed in section 6 However for the purpose of the theoretical analysis of robustness in /ection this distinction can be ignored.
4.2 Sources of Misspecification
For the remainder of the paper extend our theory and notation to encompass several (possibly simultaneous) sources of misspecification: standard GP theory assumes that data comes from a latent Gaussian random field GRFFof(((//L./1)) specified by covariance function c(.-) and parameters 1. The construction the matrix Kg in section assumes data to have arisen from this GR.F with i.dd (0,, (²) additive noise. Henceforth, we limit covariance functions c(ar, æ²) to be stationary, i.e. to vary only with (æ x'). The forms of (possibly simultaneous) misspecifications be accounted for in the theoretical treatment 5.1 are: (a) parameter o wrongly specified as o (b) (normalised) covariance function -( wrongly specified as 6(.-), (c) parameters 1, misspecified as 1, (relevant only e(:,)) not misspecified), (d) true additive noise is not Gaussian and (e) the data generated by non Con-saasiinn weakly stationary random field WSRJF rather than GR.F.
5 GP nearest neighbour Limits and Robustness
In this section we investigate the behaviour of MSE, NLL and calibration for GPan prediction as oo, showing how all of these performance measures become increasingly less sensitive to hyperparameter accuracy, kernel choice and the above departures from the GP model assumptions.
5.1 Theory
The true generative model from which the data arises is y4 = f(x.) + with & Pe. f(x) Ws and yi f(;) ~ Pş where the variance of the distribution Pg is og. Neither the WSRF nor the additive noise distribution Pg need be Gaussian. The training values are i.i.d. The MSE, NLL and calibration statistics the test set derived according to the nearest neighbour GP prediction process (4.1) and subject to any/all forms of misspecification in 4.2 Additionally, assume that if and only if the m nearest neighbour converges to the test point under c, then it also converges under è (A Definition IOP.
Assumptions:
Given a size-n training set X and test point x" in the support A Definition 9) of the measure of let ( Ey.. }, (NLL(() Ey.y. [[ (CAL(()) Ey.y { J; where expectations are w.r.t. the true generative process for y and y² and the performance measures 1N,N,,,, 2 (section 3) are for the nearest neighbour prediction process. Note that these the expected (rather than mean) values of the performance measures described in /sectionn 3 and the dependence on is implicit in the construction of the nearest neighbour sets N = -Nn (xx") used for prediction. Then we have:
Result:
Theorem (GPan limits), As - oc, JMSE JCAL NLL (i.i..)) measure on x E X and x", and pointwise as functions of where:
MSE , ,NLL ae w.r.t the
(i) MSE (0) = = + + (
(ii) (CA((O)
oz + =
(iii) for
(NLL = flog (6²(1+ m
+ 2
+ log 2 00
=0 (0mmm)))
Setting Ô = 0 (and in particular 6² = o?)) in the above provides matched matched-parameter limiting results.
derive expressions for each of the expectations JASE SCAL, ,NLL since these only depend on the known marginal covariance matrices of the (misspecified) W.RFF. We then use results concerning asymptotic convergence of Fuclidean nearest neighbours, in combination with some standard linear algebra results and continuity properties, to
Proof sketch. It is quite straightfoward

---

NLL (1=0.5)
Calibration (/=0.5)
(/=1.0)
1.25
2.0
1.50
1.25
1.00
1.5
1.00
0.75
1c+05
1.0
0.75
0.50
1c+07
0.50
0.5
0.25
0.25
(,10)
0.00
0.0
0.00
0.0
0.5
1.0
0.0
0.5
1.0
2
1
Figure 2: Behaviour of performance metrics as functions of kernel hyperparameters for increasing training set sizes n. The black dashed line denotes the true parameter value; the red dashed line shows the limiting behaviour as - oo and the green dashed line shows the limiting behaviour when the hyperparameters are correct (the red and green dashed lines coincide for MSE). True shown in the title; additionally of 0.1, of 0.9, d = 20. When not varied, the assumed parameters are 02 = 0.2, 07 0.8, i Finally we generate the input data from the measure Px = N (0, Td).
obtain the stated limits. Note that the expression for for MMS above the right hand side must always exceed of since this is an absolute lower bound on MSE performance; likewise fCAL, is constrained to be non-negative. See full proef (A).
The MSE results of Theorem show that to within small factor (e.g. m 0.0025 when 400 as for all reported runs our algorithm) the best possible MSE will be achieved in the limit. The NLL results also tell us (by setting of o what the best possible limiting NLL value is, but only according to the possibly misspecified Gaussian model. The corrupting influence of an incorrect value of o2 on the limiting NLL value clearly evident from the expression for fNLL and the picture is similar for calibration.
Interpretation:
Remark 2. Theorem shows that isotropic (e.g. RBF and Matérn) kernels converge to the best possible MSE as n. + 0O even on data generated with independent lengthscales on each x coorlinate.
Note that referr to pointwise convergence whereas we believe uniform convergence results should also be obtainable, e.g. perhaps of the form (or similar):
Conjecture Ex,x i -
(MMSE(O) + MSE = a1 + +m- Om -2) uniformly as a function of {
This particular conjecture would hold, for example, if the Lh.s. were shown to be a continuous function of I reducing monotonically and pointwise to the limit with (by Dini's theorem). We also have initial results on rate of convergence Theorem which we defer to later publication once more fully extended.
5.2 Simulation of Limits and Robustness Scale
At first sight it seems infeasible to demonstrate the above robustness and limit properties empirically on GP data-sets of size 10yh above. One major obstacle being the generation of GP synthetic datasets at this size which is computationally prohibitive even allowing for the speedups described in 29). Fortunately major efficiencies, by adopting the approach described in Algorithm I
avoid the need for large-scale data-generation, in addition to achieving other
The simulation algorithm gains its efficiency by exploiting the locality of the GPnn prediction process at x whereby the predictive distribution only makes use of size-(m + 1) marginal distribution of the full distribution of (y, y/ over 1=++1, (By the definition of a Gaussian process this marginal is (low dimensional) multivariate Gaussian distribution from which samples can cheaply be generated). The following lemma is proved in Appendix B
5

---

Algorithm Simulation of GPnn Robustness and Limiting Behaviour
Input: n (training size), n² (test size), (number of nearest neighbours), (æ-dimension), c(:;))
(generative covariance function), 0 (generative kernel parameters), P(.) (x-distribution) Set-up phase A:
1. Produce the training æ-values and n" test -t-vuuues by sampling independently from P(.).
2. Find the n² size- m nearest neighbour sets (x;)) for test points x; (i [ ,.....,""))
3. Generate n GP samples y: Rm*+1 from (0, K where U N U {
4. Store (N(x)), (ya. f=1)) which will be used for predictions in Phase B together with
y (m) E R, the true value at x...
Robustness Evaluation phase:
For several choices of assumed covariance function o(:;)) .) and parameters 0: For i 1,..., n² do:
1. Compute predictive distribution (Np 2) at xt based on è, t and (N(z)), (mn))s=))
2. Update NLL., MSE and calibration statistics using "i, of² and yf
Output: NLL., MSE and calibration stats for range of covariance and parameter assumptions.
Lemma 4 (Algorithm validity). The MSE, NLL and calibration estimates returned by Algorithm 1 are equally valid those that would be obtained by applying the full GPnn predictive process (exactly as described in |sussection 4.7) T) to synthetic data sets of size n.
Figure shsows how the observed performance metries approach the limiting behaviour as n increases. In particular, from the RHS plot we see that as n increases, MSE becomes increasingly insensitive to departure of from the true value I = 1. This is consequence of (what appears to be uniform) convergence of MSE toward constant value of 0.1 (the best achievable MSE) as predicted by a practitioner selects particular kemnel family, the accuracy of the hyperparameter i becomes less and less critical to MSE predictive accuracy with increasing n, so that expenditure on estimating it accurately provides diminishing returns.
Theorem ( In practical terms:
The interpretation of the leftmost two plots is similar albeit somewhat more involved: The dotted red lines show the asymptotic dependence on the misspecification of the noise-variance as predicted
by Theorem
(ii) and (iii) i.e. plots of y 0.1 and y flog o 0.1 de
ac log 2 where 0.1 of
is the true noise value used to generate the synthetic GP data. We again see evidence of uniform convergence toward this limiting behaviour with increasing n. The green horizontal dotted lines show the limiting values of NLL and calibration (y [log of + 1 log 2 and y 1 respectively) that can be achieved the incorrect of value is replaced by the correct value o.. This underlines the importance of estimating this particular parameter more accurately in order obtain improved NLLL and uncertainty calibration for large n. Further plots from Algorithm I showing dependence each metric on all of the parameters, are given in Figure
A Highly Scalable GP Nearest Neighbour Regression Algorithm
6.1 Parameter Estimation
Parameter Estimation Phase 1 The first step of parameter estimation (Figure involves randomly selecting a small subset of the training data to obtain a first-pass estimate (i, o2, az). Small
subsets yield sub-optimal Ô values, yet shown in 7.1 these are capable of yielding strong MSE performance due to the robustness properties of section 5) We use the method in section 3.1 of 7 to estimate parameters from E, randomly partitioning E into w size-s subsets (ws = e) and using block diagonal approximation (with blocks of size x s to the full X e gram matrix. For Table we set e 3000, s = 300, m 10. For strong computational efficiency we set = |E| to a small constant value no matter the size of n. Thus as grows, an increasingly small portion of the data is used for this phase of parameter estimation and the associated cost does not increase with n. Note that other choices of cheap parameter estimation could be substituted here.

---

Parameter Estimation Phase 2 (calibration) As shown in section 5 NLL and calibration per-
formance derived from = (1, (2,02.07 03) remain very sensitive to inaccuracies in o An additional
"calibration step" is used to refine those parameters: We randomly select a size calibration set C (which is otherwise unused) from the training data and proceed according to Algorithm
Algorithm 2 Calibration of Predictive Distribution
Input: A size subset C' of (x","") pairs from the training data, parameters 0 = (1, 0l, 03).
1. For each (æf yf) E C' use the efficient GPnn predictive algorithm of with covariance function o(.,.) and parameters e = (i. 87), to obtain an estimate of the mean and variance, |1, o,² of the predictive distribution of y at a..
2. Compute a = Ei=
(07-1) - of
Output: Calibrated parameters O' (i, (1.a..6... a o.)).
Note that this process not only adjusts the noise variance estimate of but also the kernel scale parameter o7 In doing it simultaneously calibrates the predictive distribution and improves NLL performance whilst leaving unchanged the MSE performance obtained from the original parameter estimates ① oZ of)) The lemma below is straightforward to prove (Appendix C):
Lemma 5 (Calibration). The parameters O' output from Algorithm 2 produce GPnn predictions that (a) achieve perfect (weak) calibration on C, (b) minimise NLL on C' over all choices of a and (c) produce the same MSE as does on any choice of test set.
Remark Algorithm 2 can applied other GP methods, such as SVGP, to improve calibration.
Table T Jsees 1000; a simple refinement would be to select automatically (e.g. using bootstrap) with optional manual override. Where accurate uncertainty calibration is paramount practitioners could devote much larger CPU resources to this phase (which is also easily distributed); when no uncertainty measures are to be used this calibration step can be bypassed altogether.
6.2 Efficient Nearest Neighbour Prediction
In order to implement GPan prediction described in 4.1 we use the scikit-learn NearestNeighbors package ((22) This implements an efficient approximate nearest neighbour algorithm whereby one-time work is carried out to construct table (at O(dn log (see e.g. [9. (1] and counted within the total training times quoted) which subsequent calibration/test predictions then make use of. Query compute-coste are described in the associated documentation, c.g. Old n) for the Ball-tree algorithm which the default automated algorithm selection in SciKit-Learn should at least match. In contrast, exact KNN costs are listed in the documentation ((dn). As is evident in Table 2, Figure 3 and the quoted query and table-setup complexities, the nearest neighbour work increases with di-dinnnsoon d. Alternative nearest neighbour algorithms and/or comeension-reductie techniques to help address this are yet to be investigated. We set the number of nearest neighbours to be 400 for all usages in this paper having observed minimal sensitivity this parameter on independent synthetic datasets. Although a simple cross-validation procedure could be followed if tuning of is desired (at an increase in computational overhead), we wished to minimise such fine-tuning to emphasise the simplicity and robustness of the method we present, noting the strong performance we obtain despite this. At first glance, prediction complexity might appear restrictive, but some empirical tests on a laptop reveal comparable performance to SVGP prediction (Table 5).
7 Experimental Performance of GPnn Regression
7.1 Performance on Real World Datasets
Implementational Detail!! !. Comparisons are made between our method and the state-of-the-
art approaches of SVGP [1] and five distributed methods ([1) 2 7 and [18] following the
recommendation in 4 ) We have chosen not to include other perfooperfo approximations
htttps://github.comm/me http://aatttttttehsooooeeeeeeeee htt.://sooongmmmmmooooooeeemmsss experiments/
7

---

Table 1:
RMSE and NLL results (mean and standard deviation over 3 runs) for the best distributed method (w.rt. MSE). SVGP and our method.
NLL Distributed
RMSE Distributed
OURS
SVGP
OURS
SVGP
Dataset
Poletele 4.6c+03
19
0.0091 + 0.015
+0.214 0.019
-0.0667 1 0.017
0.241 1: 0.0033
0.195 i 0.0042
0.226 a 0,0059
Bike
1.c+04
13
0.977 + 0.0057
=
0.93
0.0043
0.634 # 0.004
0.624 + 0.0079
0.606 = 0.0033
Protein
3.6c+04
9
1.11 0.0051
1.1 0.0016
0.0059
0.733 0.0038
0.666 0.0014
0.688 0.0043
Ctslice
4.2c+04
378
-0.159 + 0.052
-1.26 0.01
0.016
0.237 # 0.012
0.132 + 0.00062
0,384 = 0.0064
Road3D
3.4e005
2
0.685 0.0041
0.371
0.004
0.608 + 0.018
0.478 + 0.0023
0.351 # 0.0014
0.443 0,008
Song
4,6c405
90
1.32
1.18
0,0045
1.24 12040000012
0.851 # 6.7c05
0.787 + 0.0045
0.834 0.0011
Houself
1.6e+166 8
-1.34 + 0.0013
-1.56 + 0.0066
-1.46 0.00146
0.0626 + 5.2e-05
0,0506 = 0.00072
0.11666 0,00011
(e.g. structured kernel interpolation (SKI) methods and their extensions ((37 39 ),, since, to our knowledge, none have supplanted these methods in the community general datasets. Parameters for our method are given in [6.16 and 6.2] SVGP used 1024 inducing points; the distributed methods all used randomly selected subsets of sizes as close as possible 625. The learning rate for the Adam optimiser was 0.01 for SVGP and 0. for our method and distributed methods. All runs in Table 1 used the the squared exponential ("RBF") covariance function. A "pre-whitening" process (E.1) was applied to x values for all methods and the y values normalised (using training data-derived means and sds) have mean zero and variance 1. More complete details are given in E SVGP was run on a single Tesla V100 GPU with 16GB memory; all distributed methods run on eight Intel Xeon Platinum 8000 CPUs sharing 32GB of memory. Our method was run on Macbook Pro with 2.4 GHz Intel core i5, See D a full explanation of our selection and pre- pre-processing of datasets which, apart from Protein, are taken from the UCI repository.
ubiquitous benchmarks on
Results Runs were made on three randomly selected 7 / 2/9 splits into training and test sets. Table 1 shows MSE and NLL results for our method alongside SVGP and distributed method (note:
training set size). The table shows only the best of the five distributed methods' results (w.r.t. MSE) but full results and details of all methods and all three performance measures are given in B Complete calibration results are also plotted in Figure 4] With the exception of the Bike dataset our method is found to outperform all methods simultaneously for both MSE and NLL, and calibration likewise bar a narrow second place on the Song dataset. (Talle 2 and Figure 3] show that this is achieved whilst undereutting the training costs of the other methods, an effect that is very pronounced for large training sets (c.g. approximately 100x faster than the other methods at n 1.6 x 10" on House Electric). Figure 3/shows that a significant portion of time involves calibration; this can be parallelised (or eliminated if uncertainty is not required). Note also that larger timings observed for higher dimensional datasets are due to slower performance of the approximate nearest neighbour algorithm (6.2) in that regime, both for nn table construction and calibration. As discussed in 6.2 future improvements may reduce this effect. It is very interesting that course-of-dimensoonn has not impacted on the method's MSE, NLL or calibration competitiveness at large d. This was despite the fact that PCA analysis of the training x values showed no concentration within a low dimensional space (as to be expected given the prewhitening that has been applied (subsection E.I)
Conjecture 7. Robustness "curse-of of dimensionality" is at least partially explained by the increase in the intrinsic data-length-scade by a factor of order dd that must arise in order for GP methods to be effective.
The heuristic reasoning behind this conjecture is as follows: Unless length scale increases with d the kernel gram matrix will exhibit an abundance of exceptionally small off- off-diagonal entries and hence be unable able to gain significant predictive power. A vd increase serves to counterbalance this effect and seems consistent with length scales recovered from real data in practice.
7.2 Performance on Massive Synthetic Datasets
We generated size 5 10" datasets using the 15-variable deterministic Oakley and O' H'Haaan function 20 30) with i.i.d. variance-o additive noise sampled from zero-mean Laplacian distribution (with much wider tails than (00 of)). This function has 5 inputs contributing significantly to output
8

---

Table Corresponding recorded training times (with mean and standard deviation from runs) associated to the metrics Table I i.e. recorded the same time and with the time 10 given for the "distributed" method relating to the best per- forming model terms of MSF. Mean times are rounded to 3 s.f. and standard deviation to 2.
Model
SVGP
OURS
OURS
>20
102
Train time/s Distributed
=20
OURS
SVGP
Dataset
Poletele 4.6e+03
19
17.1 = 0.66
28.8 + 0.22
11.9 = 0.081
10a
Bike
1.4++04
13
43.5 0.64
0.12
32.3 + 0.15
104
10
10c
Protein
3.6e+04
9
98.9 1.7
27.7 + 0.19
81.1
Ctslice
4.2e+04
378
86.9 1.7
7.6
98.2
Road3D
3.4e+05
1200.0 + 110.0
27.9 1.3
760.0 8.0
Figure 3: Training times (s) for each
Song
4.6e+05
90
1050.0 + 110.0
138.0 =
1080.0
14.0
model with "high" dimensional datasets highlighted. f: without calibration.
HouseE
1.6e+06
8
3110..0 + 250.0
32.0 #: 0.34
3720.0 # 17.0
cymlog(Callbration
symi (Relative extreme)
RMSE
BCM
GPOE
GRBCM
0.6-
RBCM
SVGP
0.4-
uiee
Cenice
Housestecctic
Pnleriia
Protiin
RoadSD
Song
Bike
Cholie
Poleteie
Portein
ReadhD
Song
Bãc
Ceniic
Housieseecnn
Prosiin
No Boaern
Song
Figure 4: Experiment results on a suite of UCT datasets. Optimal calibration performance is (indicated by a black dashed line). Lower is better for NLL and RMSE. Y-axis truncated for readability for Calibration due to very large values on the "Ctslice" dataset. NLL is rescaled relative the most extreme model performance. "symlog" refers logarithmic axis rescaling applied the y-axis on both positive and negative values ("symmetric").
variance, with smaller effect, and 5 with almost no effect. These properties are poorly matched by the isotropic covariance functions being applied, resulting in gross misspecification of the assumed GRFF model and the additive noise. Figure shows performance achieved with both the squared exponential ("RBF") covariance function and the exponential (Matérn 1/2) covariance function. It is very interesting to note the improvement in convergence rate achieved by the exponential kernel. (see Remark for potential explanation of why isotropic covariance functions are so effective at large n).
Remark 8. We checked to see whether this strong exponential kernel performance extended to UCI datasets. Surprisingly, given that it is not recommended for GP regression (e.g. (36) page 85), it produced best RMSE performance across the board when compared with RBF and Matérn 3/2 kernels Table 4), with Road3D RMSE reducing from 0.351 to 0.098. Exponential- performance was everywhere best apart from Ctslice and calibration was also better in most cases.
NLL
8 Discussion
The basic "subset-of-date approximation (B3) also achieves training efficiency by using a small portion of training data and can achieve surprisingly goods results (B3), [35] section 5.1, [3). But it typically would need much greater proportion of training data than we are using for large n due to its failure to leverage the power of large training sets for prediction; this explains why it is not consistently competitive with other methods. Passing references to the decoupling of prediction and estimation have been made, e.g. [8] I 3) but not shown to consistently powerful as have found, nor justified in terms of robustness theory explored as mainstream
Related work:

---

Calibration
MSE
NLL
1.20
0.50
0.15
1.10
1.00
0.10
0.90
0.05
-0.11
0
2
0
2
4
0
2
1e7
n
1e7
1e7
ExD, of 0.05 RBF, of 0.05
Exp, of 0.1
RBF, of 0.11
Figure 5: Behaviour of performance metrics on the Oakley and O' O'Hagan function derivee dataset (20) (0) as function of data size, The horizontal green line indicates the limiting behaviour if the predictive and generative were to match. Shaded regions indicate 95% confidence intervals for the fitted curves.
approach. Various works (primarily from the geospatial community) make use of nearest neighbour (NN) techniques for GPs (e.g. [34 likelihoods for parameter estimation, whilst Stein (28) adapts this work to REML using more distant points as well as NNs, again for parameter estimation purposes. In contrast, our focus is on using NNs for prediction, and whilst we have found passing references to its explicit use for this purpose (e.g. (27)) we have found little or no discussion of effectiveness in comparison other methods on large datasets and no detailed accompanying analysis of its robustness properties and how these can achieve very high efficiency at scale. (5] gives a construction of hierarchical fully Bayesian model (" ("NNGP") derived using collections of NN sets. This approach adds considerable computational overhead and code-come colplexxy,,,,y e.g. use of Gibbs sampling. Whilst fully Bayesian treatment of hyperparameters is explored in the ML literature, e.g. [17]. it has not been adopted by the ML community for use scale due to its high computational cost relative to empirical Bayes methods (177) section 5). Bearing this in mind, we consider the extension given in S for improved scalability (Algorithm 5 "conjugate NNGP') to be more relevant. In this hybrid method some hyperparameters recovered 'empirical-Be Bayes' point-estimates via grid-based search and the remainder treated in This
28 6 8
32
38) Vecchia ((44) uses NNs approximate
a fully Bayesian fashion using a conjugate prior with some choice of nyper--yperparametoo results in Student-t predictive distribution, rather than Gaussian, but with equal first moment to ours and variance differing only by a ttyper-hyperparanee that recalibration step would render redundant (see I.... Recent work ((32 38) ) extends local geospatial GP methods into sparse variational ML. applications. 'VNNGP' is shown to be competitive with other methods in (38) despite adding further approximations into the model. We note that when using all observations as the inducing points their predictive mechanism matches ours, up to choice of parameter estimation. We believe these pre-existing works, which differ significantly in approach and perspective, complement our own, which emphasizes the benefits of decoupling parameter estimation from prediction, the robustness properties that can be achieved at large scales, the efficacy of simple recalibration step and can be run at high scale with a simple algorithm on an off-the-shelf laptop.
dependent) multiplicative factor, an effect
Limitations and Future Research: Our results exhibit leap in speed and performance for GP regression at scale, but there remains more to be done to fully explain and extend performance (as evidenced by our remarks and conjectures). This is particularly so for high dimensional problems where (a) a faster nearest neighbour algorithm would have particularly big pay-off and (b) there need to explain why "curse-of-dimensionalll competitiveness (see(Conjecture 7). Extensions of theory to broader aspects of GP robustness, rates of convergence and "strong calibration" (section 3) are current areas of some the authors' ongoing work.
appears not have damaged the method's
Acknowledgments
We would like to thank IBM Research and EPSRC for supplying iCase funding for Anthony Stephen- son and the UK National Cyber Security Centre for contributing toward Robert Allison's 's funding. We also thank the anonymous reviewers of this paper for their constructive comments and suggestions.
10

---

References
[1] F. Bachoc, N. Durrande, D. Rullière, and C. Chevalier. Properties and Comparison of Some Kriging Sub-model Aggregation Methods. Mathematical Geosciences, 2022.
[2] Y. Cao and D. Fleet. Generalized product of experts for automatic and principled fusion of gaussian process predictions. arXiv preprint arXiv:1410.7880 2014.
[3] K. Chalupka, C. K. Williams, and Murray. A framework for evaluating approximation methods for gaussian process regression. Journal of Machine Learning Research, 14:333-350, 2013.
[4] S. Cohen, R. Mbuvha, T. Marwala, and M. Deisenroth. Healing products of Gaussian process experts. In H. D. III and A. Singh, editors. Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1068-2077. PMLR, 2020-13.
[5] A. Datta, S. Banerjee, A. Finley, and A. E. Gelfand. Hierarchical Nearest-Neighbor Gaus- sian Process Models for Large Geostatistical Datasets. Journal of the American Statistical 2016.
Association, 11 111(514):800-88
[6] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand. On nearest-neighbor Gaussian process models for massive spatial data. Wiley Interdisciplinary Reviews: Computational Statisties, 15):1122177, 2016.
[7] M. Deisenroth and W. Ng. Distributed gaussian processes. In International Conference on Machine Learning, pages 1481-1490. PMLR, 2015,
[8] A. O. Finley, Datta, B. D. Cook, D. Morton, H. Andersen, and Banerjee. Efficient algorithms for bayesian nearest neighbor gaussian processes. Journal of Computational and Graphical Statistics, 2812):01-414,, 2019, PMID: 31543693.
[9] J. Friedman, J. Bentley, and R. Finkel. An Algorithm for Finding Best Matches in Logarithmic Expected Time. ACM Transactions on Mathematical Softwar (TOMS), 3:209-226, 1977.
[10] J. Gardner, G. Pleiss, K. Q. Weinherger, D. Bindel, and A. G. Wilson Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. Advances in neural information processing systems, 31, 2018.
[11] R. Gardner, G. Pleiss, R. Wu, K. Q. Weinberger, and A. G. Wilson. Product kernel interpolation for scalable gaussian processes. International Conference on Artificial Intelligence and Statistics, AISTATS 2018, 84:140771416, 2018
[12] L. Györfi, M. Kohler, A. Krzyžak, and H. Walk. A Distribution- Theory of Nonparametric Regression. Springer Series in Statistics, 2010.
[13] K. Hayashi, M. Imaizumi, and Yoshida. On random subsampling of gaussian process regression: A graphon-based analysis. In International Conference on Artificial Intelligence and Statistics, pages 2055-2065. PMLR, 2020.
[14] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. arXiv preprint arXiv:1309.688 2013.
[15] G. F. Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 147111177--1800, 2002.
[16] M. Jankowiak, G. Pleiss, and J. Gardner. Parametric Gaussian process regressors. In H. D. III and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 1 119 of Proceedings of Machine Learning Research, pages 4702-4712, PMLR, 13-18 Jul 2020.
[17) V. Lalchand and C. E. Rasmussen. Approximate Inference for Fully Bayesian Gaussian Process Regression, Apr. 2020. 1tXiv.1912.13444 [cs, stat].

---

[18] H. Liu, J. Cai, Y. Wang, and Y. S. Ong. Generalized robust bayesian committee machine for large-scale gaussian process regression. In International Conference on Machine Learning. pages 3131-3140. PMLR, 2018.
[19] H. Liu, Y.-S. Ong, X. Shen, and Cai. When gaussian process meets big data: A review of scalable gps. IEEE transactions on neural networks and learning systems, 31(11):4405 -4423, 2020.
[20] J. E. Oakley and A. O'Hagan. Probabilistic sensitivity analysis of complex models: A bayesian approach. Journal of the Royal Statistical Society. Series B (Statisfical Methodology), 66(3):751- 769, 2004.
[21] S. Omohundro. Five Balltree Construction Algorithms. Technical report, International Computer Science Institute, 1989.
[22] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-288 2830, 2011.
[23] B. Schölkopf. The Kernel Trick for Distances. In T. Leen, T. Dietterich, and V. Tresp, editors, Advances Neural Information Processing Systems, volume 13. MIT Press, 2000.
[24] H. Song, T. Diethe, M. Kull, and P. Flach. Distribution calibration for regression. 36th International Conference on Machine Learning. ICML 2019, 1019--Jue::0347-1035556 2019.
[25] M. L. Stein. The screening effect in kriging, Annals of Statistics, 30(1):298-323, 2002.
[26] M. L. Stein. 2010 Rietz lecture when does the screening effect hold? Annals of Statistics, 39(6):2795 28199, 2011. arXiv: 120v.1
[27] M. L. Stein. Limitations on low rank approximations for covariance matrices of spatial data. Spatial Statistics, 8:1-19, 2014. Spatial Statistics Miami.
[28] M. L. Stein, Z. Chi, and L. J. Welty. Approximating likelihoods for large spatial data sets. Journal of the Royal Statistical Society. Series B: Statistical Methodology, 66(2):275-22) -266, 2004.
[29] A. Stephenson, R. Allison, and E. Pyzer-Knapp. Provably reliable large-scale sampling from gaussian processes. arXiv preprint arXiv:2211.08055: 2022.
[30] S. Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved May 5, 2023, from htttp://www.sful wtw.ss/suu.aasurrrrraaaa ano.
[31] M. Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
[32] G.-L. Tran, D. Milios, Michiardi, and M. Filipponee Sparse within sparse gaussian processes using neighbor information. In M. Meila and T. Zhang, editors, Proceedings of the 38th Inter- national Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 10369-103778.. PMLR, 2021-18.
[33] V. Tresp. A bayesian committee machine. Neural computation, 12(11):2719-2744 2000.
[34] A. V. Vecchia. Estimation and Model Identification for Continuous Spatial Processes. Journal of the Royal Statistical Society. Series B (Methodologicall, 50(2):997-31,, 1988. Publisher: [Royal Statistical Society, Wiley].
[35] K. Wang, G. Pleiss, J. Gardner, S. Tyree, K. Q. Weinberger, and A. G. Wilson. Exact gaussian processes on a million data points. Advances in Neural Information Processing Systems, 32, 2019.
[36] C. Williams and C. E. Rasmussen, Gaussian processes for machine learning, volume MIT press Cambridge, MA, 2006.
12

---

[37] A. G. Wilson and H. Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). 32nd International Conference on Machine Learning, ICML 2015, 3:1775-1784, 2015.
[38] L. Wu, G. Pleiss, and J. Cunningham. Variational Nearest Neighbor Gaussian Process. Proceed- ings of the 39th International Conference on Machine Learning, 2022.
[39] M. Yadav, D. Sheldon, and C. Musco. Faster kernel interpolation for gaussian processes. In A. Banerjee and K. Fukumizu, editors, Proceedings of The 24th International Conference Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 2971-2979. PMLR, 13-15 Apr 2021.
13

---

A Theoretical GPnn Results
A.1 Preliminary results
be the kernel-induced distance function over Rd ([33 We define X(j.n) (a") as the jth nearest neighbour random variable to a test point x under p. which the realised jth nearest (x)) ø") and
Let p(z, x') = -oo//-ccx/l,/////
we abbreviate to X(j) when the context is clear, and (z") E (x""
neighbour of the test point x: from training set X. From this we define
Eij (x)) X6)).
Definition 9 (Support). Let Px be the probability measure of and S the closed ball of radius > 0 under the metric centred at æ. Then we define support (Px)rtt(ppp))))))))))>) = {x Px(S2,) 0.0 > 0}.
Definition 10 (Weakly-faithfull We define a pair of metrics p(::),, )( be weakly-faithful w.E.t. each other if the following condition holds: The we nearest neighbour under á converges to the test point as + do if and only if the thth nearest neighbour under p converges to the test point in the limit.
Assumptions
sid (A1) x Px and x" E support(Px) under the generative metric defined by c(.,).
(A2) c(; )))) are stationary kernels whose induced distance functions are weakly faithful metrics (Definition 10)
(A3) = f(x:;) + with & (( ( and yı f(m: Pe and [ = 0, - o.
Note: Assumption (A2) is not overly restrictive and encompasses commonly used kernels such as all those mentioned in this paper.
Lemma 11. e; + 0 and eij + 0 aasnn + oo a.e. with respect to the measure over x E Rdd for
i, < m, - 0 and under (A1-2)
Proof. Lemma 6.1 of [12] states that [x(m,n) (2) - æ0|| n=200 0 with probability one (with respect
Px). Their proof can be generalised immediately state that ()) m=000 0 by using our definition of support. / that directly invokes the metric Hence e; -> 0 for all i S m (since xx is in support( (Px)). Since is metric it satisfies the triangle inequality; hence p(X()) X()) S
p(x()): x²) + p(x()). x²) R-000, for all i, S m.
Lemma 12. For an m-GPnn under the assumptions (AliicoonsSAAA (1-3)
lim k K= kiy = oî -o?m +0m
n+000
Proof. From Lemma 11
we have that limn-roc k(x()) (x=), x²)" = limm (07 - <,) of and
lim (x=)) (6) (x")) = lim, -00n (0) - Cij) = 07. As result, ky + 031 and
K = lim KN of l + 1111
(8)
n=000
14

---

Now using Sherman Morrison and the continuity of matrix inverse and matrix matrix-matrix products:
A A bc A-1
(A bc"))===-=" A-1
(9)
1 c
11
1K²0)-1-002 10)11111)-11
0f o + of
2
(10)
mof
T 1k(K0
m
1 0 mo?
+ oe
m
1-mot
1 - mot
1 1
0mh
mo}
fom
m²o}
1
30
+ 0(m
(11)
0
mo
4
Thus,
lim n-000
k K kiy (K>)) = o3 - omm -1 + 0m
(12)
Lemma 13 (WSRJ expectations). Under (A3) Ey.y {yv"} k' and By (vy)) } K.
Proof. By assumption the covariance properties of y and the independence and zero-mean of the additive noise, Fy [v:y}} = k( æ.)) Extending this to the joint distribution over y, y is straightforward and gives the results stated.
Lemma 13 is subsequently assumed to be in use throughout A.2
A.2 Limit proof's
In the following statements only misspecificationn of type (d) and/or (e) (subsection 4.2 is considered be at work.
Lemma 14 (MSE limit). Under the assumptions (AI-3) for fixed m <, oc, the predictive GPan given in subsection 4.1 converges pointwise in the sense of MSE wrt P Px-a.e as
lim p (MSE()) n
= 1 + m m-) 0 (m-
Proof. This follows from Lemma 12 by expanding the definition of MSE:
MSE (0) = lim Ey.y" {lu- -HN1?
lim f 2=:000
n n00
lim
Ey 22²2) + Ey uuN" 23
Ey. (KW KN "y''
= of + o
lim Fyy (1ix 2) n+30
= o2 m 1 O (m -)
Since Ey (ny Ey ( K yNyR kg) Ey.- (yNU" } ky, even under a WSRF generative model Lemma 13)
T-T K k and by assumption
Corollary 15 (NLL limit).
lim (NLL 1 #=000
1
log (o log(o 1mmmm111 1))
log 2n (mm
2
2
15

---

Proof. The proof follows straightforwardly from Lemma 12 and because oN = oî + 0g k Ky kiv
2 B y. (1v) = Ey log oN + (v- - -N3)
+ log 2nT
oF2²
log N
2
+ 1 + log 2T
lim 2 Ey ('y) log j + oz (o omm - + O0 2m)= + 1 + log 2T
n-+oc
log (0² + m mm1m ) O m 2) + 1 log 2T
= log oî + m 1 + 1 + log 2T - (m" 22
A.2.1 Full misspecificatiion
For the remainder of A.2 we assume that the full range possible misspecifications ((a)-ee) outlined in infsubsection 4.2 are in action. We refer to this case as fully-misspeccfied" and introduce the notation jr, Nf be understood to mean the predictive mean and variance under these misspecifications. Lemma 16 (Fully misspecified MSE limit). For a fully misspecified model, asymptotically
lim fa ( =1 (1(++mm +m= 3m
provided the misspecified kernel distance metric is weakly faithful in the sense that the mh nearest neighbour converges under both the true and misspecified metrics @Definition 10p.
Proof.
Ey Ey. 1(0" - 1xx) ly]
Eyy
Eyy
2y AN + (AN)) y 1
Ey {oi²
+ 2
2p 2p + (1xx)}
NN
k
K ky-
K k
KN t
(b)
(0)
(d)
We can use standard results to state that (a) + (b) of + oz. expand it in terms of m-1:
Then we define f
87 mm+ and
89
o o
my
+ (m-
mo²
m²8
In manner similar Lemma 12 we use this result to compute:
lim (c) = 03 2(1 - ;11 116}
0700
m(1
m³)
6²
00
02
oz
+ (m)))
02
87
mot
0
of
0 00m0mm
m
f
16

---

and
0, -1 (o +0
lim
(d)
R-000
6
1"
[o]] + o3 202 + 3 of 203Am11" + 103 1
+ mlmm
moo
2m1 + m2
mmoo)( + ml0
mf)²
g
0
0
2 oh
mldg
2
mmo
+ -00mm)))
m20
m²a
of 02 of
of
+ {m
where we have used the expansion of 1 - mř given earlier. Putting these results together gives
lim fn
MMS (0) = lim ((a) + (b) 2(c) + (d)
20000
n-+00
o of
o
o² o
= = + oz 2 0)
+
o
2 -2 = 0 3)
of
= oe + m ) + O( (2-
Lemma 17 (Calibration limit under full misspecificaaion).
o²
lim (
+ =0(m-2).
0²
Proof. We use continuity to write
lim Ey.y (y*
13
( (im MSE (0)
lim
12-3000
02²
2
000
(n
By direct application of Lemma 12 3
n 076 + m (m and thus
CCL 0 = o²
lim
f
+ =m
n-000
Corollary 18 (NLL limit under full misspecification).
1 0²
1
lim (LL(() (Õ) log (o² m -)
log 2 1 0m
+
2 0²
Proof. We start with
2/ NLL (0) = Ey.r log 22? + (y* - 1Nx)
+ log 2T
032
For the second term we use Lemma 17 that have
0²
lim 2f NLL = log o² + m + R+000
+ log 2n 0mm (2)
We construct the proof using all of the intermediate results given above. In particular item followss from Lemma 16 item (ii) from Lemma 17 and item (iii) from Corollary 18
Proof of Theorem /
17

---

B Validity of Algorithm 1
(Proof of Lemma 4
Proof of Lemma 4 identical.
We prove the result for MSE, the proofs for NLL and calibration being essentially
The proof involves showing that Algorithm is exactly equivalent to Algorithm lb al- ternative scheme which self-evidently provides valid estimates of MSE on size-n synthetic datasets:
Algorithm 1b Simulation of GPnn Robustness and Limiting Behaviour (Expensive)
1. Generate size-n set of i.i.d. training æ-values X constant henceforth).
{ with x. Px (X is held
2. Generate n i.i.d. test points x. ~ Px-
3. Generate n" separate independent synthetic GP samples y: R(++1) corresponding each of the values (2;, X).
4. From within X take the m nearest-neighhours Naa)) to the test point x;.
5. Corresponding to each test point o; evaluate the function eq' = fi(y.) = (yt uxx (3())² Where, as defined in Equation 6 py k R y with N = N (æ,) and where y E R is formed from the components of y. corresponding to (2;)).
6. Compute the average + " El e,' obtain the MSE statistic.
Algorithm 1b cleclerly provides a valid evaluation of the MSE arising from the full GPan prediction process described in section 3) (albeit at prohibitive expense for large n). Hence it is sufficient to prove that Algorithm iis equivalent to AAlgorithm 1b We do so by applying two minor alterations to Algorithm lb in succession whose combined effect is to convert it to Algorrihhm 1. We also show that throughout this process equivalence is maintained with Algorithh 1btthus completing the proof:
Let yi = (vf y', y where (y, y E (m++) with y corresponding to the x-values in N N(x)). ). Since each function fe(y)) in line 5 above only depends on the components (yf' yf) of the (n + 1)-long vector yi we can equally write ex' = f. v/) and hence truncate each of the vectors y: to (vi, y/) before evaluating ex" in and this clearly leaves the output of Algorithm 1B completty unchanged.
Change 1:
Change 2: Note (by a standard property of multivariate normal distributions) that each of the truncated vectors (uf y/) are i.i.d. from the multivariate normal distribution whose covariance Σ is the kernel gram matrix corresponding to just (æ,, N (Nx))). Hence it is legitimate to change the (inefficient) way of sampling (u," v/) (via the large sample y) direct sampling from (0 2) whilst still maintaining equivalence with Algorithm 1b
In this way we arrive at Algorithm
C Parameter Calibration (Proof of Lemma 5)
Proof of Lemma 5 (a) Replacing parameters Ô - (i, 03) with 0' - (i, ad? adz) changes all of
the ot² values to a and therefore changes the calibration value on from 1 c E
(v. (m)=
o7
a/a 1. (b) The NLL on C' arising from parameters a6, ad?) is 1. flog [oo? (vt 2 + log 2 which, on taking first and second derivatives w.r.t. a, is found to be uniquely minimised by Di= o;
(c) It is easily shown that replacing parameters (og 07) by koō;,ko?) (for any k > 0) in the formula for 2² [Equation 3 and Equation 6) does not
18

---

alter n². Hence the value of MSE = /: Σi= (1) 14)2 on any size-n" test set is unchanged when parameters O' are used in place of Ó.
D Real world datasets
We consider variety of datasets from the standard UCI machine learning repository These datasets are commonly used in the GP literature (see [10] for instance) and are, in principle, easily available online. In practice, encountered some difficulties: the dataset documentation is often limited; the dataset names commonly used in other published papers do not always match the UCI database naming and important details about data pre-processing, which features to use etc, often omitted. There are numerous attempts on GitHub and elsewhere at cataloguing these datasets along with any pre- pre-processing. however we had limited success using them, with many appearing unmaintained. Our focus in this work is on testing our methods on a variety of real world datasets and in way that is, as far as possible, consistent with other papers. We therefore rejected datasets about which there ambiguity over the correct features to use, or even which column to regress on or for which outlier rejection required but undocumented elsewhere.
Referring to the datasets used in 10 we were able to locate the following:
SearPreedctt.onlSDD.ttt.....
Bike
http:///.......ee.oo..//////aaad Bike-ShariingDDaastt..ijj
Poletele
Keggundirected
machine-
.data)
CTSlice http:://heiiiii....ttaa/a//aaaat
Road3d http:://ccciiii.cc...cc/./////at
Protein
http:://cthaaa.aaaaaoaoaoaaooaoa 00265/CASP csv)
regression.. tz)
HouseFllectric
(HouseE)
We were unable to find any documentation on the Kegg datasets to indicate which of the columns should be used as the independent variable (the regressor) and neither is this mentioned in any literature of which we are aware. Initial runs of standard exact GP training and prediction produced RMSEs much higher that reported in [0]. Combining these two observations, we chose to exclude both Kegg datasets. Likewise we faced problems with Buzz. An analysis of the values revealed small proportion of extremely large outliers that we found could unduly distort performance results (e.g. depending on whether these outliers appeared in the test set for some of the random splits). With the lack of documentation we were unable to identify an outlier rejection scheme that we were confident would be consistent with results quoted in other papers. For this reason we have excluded Buzz.
The choice of (x, value that we applied for each of the used datasets is as follows:
https://archhiv--etaauu wcw.uci.si....o.c edu
accessed April 2023.
19

---

Song. The first column is all remaining columns are
Bike. We use hour. csv. The value is cnt. dteday (the date) is transformed to just be the integer representation of the day. instant just index so is dropped. registered and casual are dropped as registered + casual = cnt.
Poletele. The y value total total_UPDRS. _PDRS. The columns subject# and test_ time are not relevant to the problem so dropped.
CTSlice. value is the final column. The first column is dropped as it is just an index. We additionally drop six columns which are constant over the majority of the dataset, namely columns 59, 69, 179, 189, 279 and 351.
Road3d. value is the final column. The first column is dropped as is just an index.
Protein, This dataset was processed as per htttps///githut com/ hughsal imbeni/ whereafter we used our own random (seeded) train/test split.
boyesian_benchhag
HouseElectric. y value is the column labelled "Global active power". rescaled by 1000//0 and with "Sub metering 1,2,3" columns subtracted. We convert the date column into day-of- year/365 and the time column into time of day in minutes. Further, we remove any rows with null entries.
We note that although we are using standard set of real-world datasets, not always clear exactly how others in the field have carried out their own preprocessing, limiting the ability to make direct comparisons to other results reported in the literature.
E
Additional Implementation Details
E.1
Pre-whitening of Data
For all datasets covered in subsecttionn 7.1 following "whitening" preprocessing step is adopted: Let y be the vector of all regressor values in the training dataset only, and X the matrix of all regressands in the training dataset only, where each row of X is a feature. Let Myy og be the sample mean and variance of y respectively in the training dataset, then the whitened y values used both the training and test set are simply oyy (y My)) Let μx Σx be the sample mean and covariance matrix of X respectively Let Σx M MMT then the whitened values in both the training and test data are 1. M- (x: xxx), ) where d is the feature dimension of X. Note: the performance metrics v7
given in subsection 7.1 are expressed in terms of the whitened values rather than the y values in their original form. This appears to common practice the literature and has no bearing on the comparative performance of the different methods within this paper.
E.2 Test-Set Batching
To prevent excessive memory consumption, we perform all predictions for the distributed and variational methods in batches of 1000 points at time. Where this is not possible (e.g. for especially large datasets), we use smaller batches of 500 or 250 points, as appropriate.
E.3 Additional Implementation Details for SVGP
We use the sparse variational inducing point approach of (14). following the implementation provided by GPyTorch, which in particular uses a Choleksy decomposition to parameterise the covariance matrix of the variational prior. We broadly follow the SVGP implementation exam- ple provided by https: //do:s///ccs.. spytorch. ai/en/eeeeeeeeaa htt/e//stable//////e/e/// www.s//eeeaapeee////ee/ Variatiooallll.........aaa and__ Approximate GPs SVGP Regression_ CUDA, html In particular, we follow their example in us- ing the Adam optimiser to train our model over 100 epochs with a minibatch size of 1024 and learning rate of 0.01. We opt to use 1024 inducing points. All experiments under this method are run on SageMaker mll.p3.2xlare instance, consisting of single Tesla V100 GPU with 16GB of memory.
20

---

E.4 Additional Implementation Details for Distributed methods
A good introduction to distributed methods for Gaussian process inference is 77]. Here we run the product-of-experts (PoE) [[3], generalised product-of-experts (gPoE) [2] Bayesian committee machine (BCM) [33. robust Bayesian committee machine (rBCM) 7 and generalised robust Bayesian committee machine (GrBCM) [8] following the recommendation in 4 to aggregate in f-space. There are three components to any distributed method: the hyperparameter inference, the partitioner and the aggregator. Hyperparameter estimation is the same for all the methods: we use the method in section 3.1 of 77.. randomly partitioning the entire training set into subsets of size 625 (or as close as possible with equal-sized experts given that in general n is not a multiple of 625). A block diagonal approximation (with n2/225 blocks) is then used approximate to the full n gram kernel matrix. To recover hyperparameters with this we use Gaussian Process models with a zero prior mean and scaled square- expooeeexppoential kernel. Training is conducted using the Adam optimiser with a learning rate of 0.1 over 100 optimiser iterations. Once the hyperparameters are trained, we distributed prediction mechanism to evaluate performance against the test-set. The 625-sized partitioned blocks are referred to as "experts" and the shared hyperparameter values are distributed to each expert and held fixed thereafter. In the aggregator, distributed prediction phase, each expert produces an individual predictive distribution and these are then aggregated to final predictive mean and variance for each of our test points. GRBCM prediction is little more complex than this makes use of an additional "communications" expert as explained in [8 aggregating in f-space as recommended in We provide timing statistics for training these models.
We use our own GPyTorch-bare based implementation of distributed GP approximations. All exact GP calculations performed using GPy Torch using the default settings (so 20 Lanczos iterations throughout and a CG tolerance of for hyperparameter inference, and 10-3 for posterior predictions). For all of our experiments, we utilise AWS 13. lx instance (consisting of 8 Intel Skylake Processors and 32 GB of RAM).
E,5 Reproducibility
All code used to generate tables and figures in this document can be found at bestps://githut com/ htt/aanttsteppensoono /ttp//gnnnnngge peppn-experimentie
F
Related Work
F1
NNGP
An hierarchical Bayesian approach to nearest neighbour GPs is derived in [5) who construct a full stochastic process allowing end-to-end probabilistic approach they term 'NNGP' derived from a 'parent' GP using collections of nearest neighbour sets forming a 'reference' set. This work is Bayes/aulll--aaaeeeeeaa
modified in various ways in [8] including by adaption to a hybrid approach to improve scalability. This is described in Algorithm 5 which presents an MCMC-free approach ('conjugate NNGP') where some parameters are estimated via K-fold cross-validation and the remainder are given conjugate priors to allow exact posterior inference. Under this model, the marginal predictive distribution is Student-t with mean and variance expressions that match that given by GPnn up to scaling factor on the variance. Here we will go into more detail to make this connection explicit. First of all, note that they use an alternative parameterisation to us, specifying an inverse lengthscale (Φ), ratio (a) of noise variance (r²) to kernelscale (σ²) and vector of coefficients for their linear mean function (3, which we can neglect since we focus exclusively on mean-zero GPs).
21

---

In Algorithm 5 they go to give the following expressions for the predictive mean and variance:
2 = (s, (s, k))
solve (M (s)
= y(s) - dot(r(s), g(k)) + dot(w, (v(N(s, (oo(
u s(s) dot (s,, w)
Un = dot(u, u, gemv u)) + + 2)
Var(y(s)) = 6y (4)
We can see that the predictive mean me matches our own in the mean-zen setting. In particular, in this case the only non-7ero term is dot(w, y (s, k)]) which can be re- re-written in our notation as k K yv- Note that this expression is, for fixed o, independent of choice of kernelscale. This can be seen using the notation given in section k K yN = of (ojON cf (C'x +al)-yy.. yy.
In addition to this, the predictive variance is equivalent to ours up to multiplicative scaling 8y (k)) (R) 1)]; i.e. Var(y(s)) 2 which be seen as follows: In the mean-zero have Mp = 1+ - (ootww, 2) which can again be re-written in our notation as 1+0ğ/oj -cy Ca.N (where Ca.N CN + al). Hence, if we rescale this by factor of of we exactly obtain the GPnn
predictive variance, so that Var(y(s)) bog (k)ve/ claimed. Since GPnn
(a) (k) - 1) by (1)) additional recalibration step rescale all predictive variances by single shared multiplicative factor, the factor S between the two methods can be made effectively redundant: i.e. if our recalibration step (Algorithm 2p were applied their method, the methods would become equivalent with respect to RMSE and veak-calibration (setting aside differences in parameter estimation and their associated costs).
1)] as
We reemphasise that whilst GPnn and Conjugate NNGP pointwise predictions can be related as above, the latter is given in the form of t-distribution which, even with matched first two moments, will have a different shape to that of a Normal distribution. As such, the (Gaussian) NLL performance measure is longer valid choice to assess this model. RMSE and weak-calibration remain distribution-agnosticc however.
22

---

G Further simulation results
0.8
1.00
0.6
0.6
0.75
MSE
0.4
0.4
0.50
0.2
0.2
0.25
0.0
0.0
1.8.00
2
0.0
0.5
2
2.0
1.5
1.0
1.5
1000
1.0
10000
1.0
100000
1000000
0.5
10000000
0.5
(.10)
(.10)
0.0
0.0
0.0
2
0.0
0.5
1.0
2
Calibration
2
1.5
1.0
0.0
0.0
0.5
1.0
0
2
67
6
7
Figure 6: Behaviour of performance metrics as functions of kernel hyperparameters for increasing training set sizes The black dashed line denotes the true parameter value; the red dashed line shows the limiting behaviour as > co and the green dashed line shows the limiting behaviour when the hyperparameters are correct. Simulations run with d 20, 0.5, o 0.1, oj = 0.9. Assumed
parameters when constant: oz = 0.2, 03 0.8, i = 0.5.
23

---

H
Further Results on UCI Datasets
H.1 Results for all distributed methods
Table Results for all methods on all metries.
Calibration
NLL
RMSE
Dataset
Model
BCM
1.02 # 0.02
1.0 0.0065
0.66 = 0.0043
GPOE
0.873 + 0.012
1.03 0.0069
0.664 + 0.0054
GRBCM
0.893 0.014
0,977 = 0.0057
0.634 + 0.004
Bike
1.4e+04
13
OURS
0.974 0.087
0.953 = 0.013
0.624 0.0079
POE
1.03 1.022
1.01 # 0.0083
0.664 0.0054
RBCM
1.01
1.0 0.0065
0.659 = 0.0043
SVGP
0.898 0.011
0.93 0.0043
0.606 + 0.0033
BCM
5.04 0.28
1.43 + 0.13
0.311 = 0.0052
GPOE
0.435 0.013
0.422 0.0015
0.347 0.0027
GRBCM
0.11
-0.159 = 0.052
0.237 +
Ctslice
4.2e+04
378
OURS
1.04 + 0.085
1.26 0.01
0.132 + 0.00062
POE
6.39 0,27
0.347 + 0.0027
RBCM
4.26
0.987 0.11
0.28 0.0048
SVGP
0.865 0.026
0.467 = 0.016
0.0064
BCM
1.27 = 0.0046
-1.33 0.0009
0.0634
GPOE
0.908 0.0065
0.006
0.0638
GRBCM
1.25 0.011
0.0099
0.00220
Houseelectric 1.6e+06
8
OURS
1.08 = 0.21
-1.56
0.0506
POE
1.28 0.006
0.0018
7.76-5
RBCM
1.24 0.0054
0.0013
2 5.2e-05
SVGP
0.911
1.46 0.0066
0,0566 - 0.00011
BCM
1.07 = 0.029
0.00035 + 0.019
0,243 + 0.0048
GPOE
0,917 + 0.02
0,246 + 0.0038
GRBCM
0.872 + 0.024
0,0091
0,241 0.0033
Poletele
4.6e+03
19
OURS
1.03 0.073
-0.214 0.019
0.195 + 0.0042
POE
1.1
0.00772 0.016
0,246 0.0038
RBCM
1.08 + 0.029
0.00309 = 0.018
0.243 + 0.0048
SVGP
0.862 + 0.035
-0.0667 0.017
0.226 + 0.0059
BCM
0400000997
1.14 0.003
0.0022
GPOE
0.925 + 0.007
1.15 0.0035
0.763 + 0.0024
GRBCM
0.95 + 0.012
1.11 0.0051
0.733 + 0.0038
Protein
3.6e+04
OURS
0.991 0.029
1,01 0.0016
0.666 0.0014
POE
1.07 + 0.0088
1.15 0.0033
0.763 0.0024
RBCM
1.03 = 0.0096
1.13 0.003
0.752 0.0022
SVGP
0.908 + 0.016
1.05 # 0.0059
0.688 + 0.0043
BCM
1.01 0.017
0.753 # 0.007
0,514 + 0.0035
GPOE
0.756 + 0.012
0.819 0.0054
0.529 + 0.0037
GRBCM
0.873 0.011
0.685 0.0041
0.478 0.0023
Road3D
3.4c+05
OURS
+ 0.041
0.371 0.004
0.351 + 0.0014
POE
1.07 + 0.019
0.783 # 0.0076
0.529 # 0.0037
RBCM
0.016
0.735 # 0.0066
0.505 # 0.0034
SVGP
0.9 0.00094
0.608 = 0.018
0.443 + 0.008
BCM
1.56 0.0063
1.32 0.0012
0.851 + 6.7c-05
GPOE
0.926 0.00049
1.27 3,4e-55
0.864 + 7.5c-05
GRBCM
1.61 0.11
1.66 0.058
0.961 + 0.035
Song
4.6e+05
90
OURS
0.99 + 0.037
1.18 0.0045
0.787 + 0.0045
POE
1.61 + 0.0067
1.34 0.0013
0.864 + 7.5e-05
RBCM
1.56 0.0062
1.31 +30001
0.851 + 6.4e-05
SVGP
0.991 0.02
1.24 0.0012
0.834 0.0011

---

H.2 Performance of different kernels
Table 4: Results on the UCI datasets using different kernel choices for our method and demonstrating the apparent superiority of the exponential kernel in these cases.
Calibration
Distributed
Ours (Exp)
Ours (Matérn)
Ours (RBF)
SVGP
Poletele
4.6e+03
0.872 +: 0.024
0.994 = 0.15
0.971 = 0.13
1.03 0.073
0.862 #: 0.035
Bike
1.4+004
0.893 + 0.014
0,988 0.098
0.971 0.086
0,974 0.087
0.898 0.011
Protein
3.6c+04
9
0.95 = 0.012
0.995 0.038
0.993 = 0.031
0,991 = 0.029
0.908 # 0.016
Cislice
378
0,082
0.085
0.865
Road3D
3,4c+05
2
0.873 + 0.011
1.09 + 0.065
1.0 0.054
0.991 0.041
0.9 0.00094
Song
4.6e+05
90
1.56 0.0063
0,995 0.033
0,994 0.035
0.037
0.991 : 0.02
Houseelectrie
1.6c+06
1.24 a 0.0054
1.11 0.29
1.08 0.27
1.08 0.21
0.911 0.038
RMSE Distributed
Ours (Exp)
Ours (Matérn)
Ours (RBF)
SVGP
Polltelle
4.6c+03
19
0.241 # 0.0033
0.169 = 0.0076
0.17 0.0076
0.195 = 0.0042
0.226 + 0.0059
Bike
0.565 0.0036
0.6 0.0044
0.624 0.0079
0.606 0.0033
Protein
3.6c+04
0.733 +: 0.0038
0,58 0.0068
0,629 0.004
0,666 0.0014
0.688 a: 0.0043
Cislice
Road3D
3,4c+05
2
0.478 + 0.0023
0,0976 + 0.013
0.27 + 0.01
0.351 0.0014
0.443 + 0.008
Song
4.6e+05
90
0.851 + 6. 75
0,776 = 0.004
0,778 = 0.0045
0,787 0.0045
0.834 + 0.0011
Houseelectric
1.6e406
0.0626 5.2c-05
0.045 0.00025
0.0485 0.0004
0.0506 0.00072
0.0566 0.0001
Distributed
Ours (Exp)
Ours (Matén)
Ours (RBF)
SVGP
Dataset
13
Poletele
4.6c+03
19
0.0091 =: 0.015
-0.397 0.028
0,346 0.032
.0211 =: 0.019
-0.066 = 0.017
0.0057
0.0077
Protzin
3,6c+04
9
1.11 = 0,0051
0.853 0.013
0.95 0.0061
1.01 + 0.0016
1.05 0.005
Cislice
4.2e+04
378
+0.159 1 0.052
-1.05 1 0.027
-1.31 + 0.017
-1.26
0.467 1 0.016
RoudDD
3.4e+05
2
0.685 0.0041
-0,931 0.14
0,109 0.039
0.371 0.004
0.608 0.018
Song
4.6e+05
90
1.32 + 0.0012
1.16 0.0046
1.17 0.0051
1.18 0.0045
1.24 0.0012
Houseelectric
0.0066
-1.46
H.3
Prediction times
Table 5: Prediction times (in seconds) for GPnn and SVGP with 400 nearest-n and 1024 inducing points respectively, over a small range of dataset sizes and dimensions.
d
GPnn
SVGP
4.2e4
378
0.06
0.02
3.6e4
0.02
0.02
1.le5
50
0.03
0.06
1.le5
10
0.02
0.06
25

---

I Overall Computational Expenditure
Our distributed and variational method experiments were conducted using cloud computing resources. Experiments using our own method have been carried out on author's 's laptop. SVGP experiments were run using SageMaker virtual machine on a single Nvidia Tesla V100 GPU with 16GB memory. Distributed method experiments were run using eight Intel Xeon Platinum 8000 CPU cores (13, (13.2xlarge EC2 instances).
Below we will attempt to give reasonable indications of the amount of computational work expended obtain the results in this paper, though note that development and research stages that did not directly contribute to the runs in the paper. As such, the costs presented are representative of the costs of replicating our paper, repeating the research from scratch. Instead of reporting costs in dollars, we will report approximate computing hours for each instance type. The reader can then estimate their costs using the current instance costs in the region of their choice, or under other cloud providers or even using on- nn-rremise compute.
neglecting the work expended in the
Dataset
Billed hours (1 GPU, 3 runs)
Billed hours (8 CPUs, 3 runs of' 5 methods)
bike
0.027
0.222
ctslice
0.082
0.546
houseelectric
3,713
256.776
poletele
0.010
0.079
protein
0.068
0,680
road3d
0.635
31.674
song
0.904
12.237
This gives a total of around 5.4 hours compute time on a GPU VM and 302.2 hours on CPU VM.
26


Reviews:
[{"review_id": "ZIQ2SpMiSy", "reviewer": "Reviewer_jFQi", "summary": " This work at one level is about proposing a scalable GP approach called GPnn where the prediction step uses only the nearest neighbours of a test point in order to form the predictive distribution. This approach is well-motivated as the authors argue that there isn't a strong mathematical reason to couple the training step / estimation of hyperparameters and the prediction / generalisation over unseen inputs. At the end of the day, we care about the latter (prediction performance) and not really the former. At another level, there is a theoretical component where they show that predictive distributions obtained through GPnn are robust to a wide array of model misspecification, for instance, wrong kernel choice or Gaussian noise in the large data limit. ", "questions": "\n- What about a sensitivity analysis to the choice of M? or are GPnn predictions insensitive to this in the large data limit. \n- Can M be learnt instead of set? \n- Perhaps some areas of the input space require more M for good predictive performance than areas where the function isn't changing much.\n- It would be interesting to see the performance on non-stationary functions.  \n- There isn't sufficient insight on the selection of nearest neighbours  exact or approximate - at first blush one would dismiss the idea because of the added compute entailed in finding neighbours per test data point. What is the complexity of this as a function of N and d precisely?", "limitations": "There is some discussion of limitations.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and clear. There is always a demand for new GP approximations which scale to millions of data points as SVGP has a known issue with over-estimating the aleatoric uncertainty. The algorithms clearly describe the simulation and evaluation step. There is a good summary of other distributed methodologies and clear delineation of how this work is different. Some components from this paper can be applied to other methods in order to refine estimates and calibration (Remark 5).\n\n  ", "weaknesses": "Figure 2 is hard to read - there is a lot going on and ideally should have been plotted with N on the x-axis. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "0. We agree that the figure is not the most straightforward to interpret, but since the intention was to convey the dependence of the performance metrics on the estimated parameters, for varying dataset sizes, reproducing the figures with n along the x-axis would lose much of this information (please see the pdf attachment for our best attempt). For example the third plot in the paper shows the  tendency, with increasing n, for \"MSE verses $l$\" to flatten toward the limiting horizontal line  (i.e. for MSE to become insensitive to $l$ ). For the final version of the paper we will modify our caption and text to clarify, and may include an additional figure with n along the x-axis, although we suspect the information this conveys will be largely covered by Fig. 5 (which will now reference in 5.2).\n\n1. Theorem 1 in the paper indicates the degree of sensitivity of prediction performance to m in the large n limit, e.g. for MSE, noise variance scaled by 1/m. For finite n it might be interesting to do such an analysis. We chose the value of m in our experiments without much experimentation beyond a brief comparison on otherwise unused synthetic data. We found minimal impact on performance. For even larger datasets it may be possible to reduce m and decrease training and test time at little cost to performance, but we did not investigate this further since we wished to emphasis the minimal amount of tuning our method requires.\n\n2. Yes, for example by cross-validation, as is commonly done in the literature (on kNN). This would of course add additional computational cost, but could be done efficiently by starting with a maximum number, and iteratively evaluating performance on decreasing subsets. \n\n3. This is an interesting but non-trivial suggestion which we would like to explore in future refinements to the method.\n\n4. Our current theoretical analysis relies on the stationarity of both the generative and predictive model, although we suspect this could be relaxed. We appreciate the potential interest in non-stationary functions, but do note that results using the Oakley and O’Hagan function (non-stationary) are included in 7.2 and provide initial encouragement. \n\n5. The literature on wide ranging approximate kNN methods is extensive and it is difficult to find a general complexity for the procedure. We used the SciKit-Learn implementation, whose costs are described in the associated documentation, e.g. O(dlogn) query compute-cost for the Ball-tree algorithm which the default automated algorithm selection in SciKit-Learn should at least match. In contrast, exact kNN costs O(dn) which is why we have chosen approximate kNN in preference. In practice the predictive timings are fast enough for most applications, taking on a laptop, for example, approximately 0.06s per prediction on the high-dimensional (d=378 and therefore relatively extreme) Ctslice dataset. We will add both example timings and the N,d complexities to the paper, which should also reassure readers that the cost of kNN search is not detrimental to the test-times of GPnn in practice (see also our response to reviewer N8D3 bullet point 4 and table 1 from the pdf)."}, {"review_id": "5E4mJmWp6d", "reviewer": "Reviewer_MqBs", "summary": "This paper starts by proposing use of a simple nearest neighbors scheme for prediction, where instead of using the training set, one uses the $m$ nearest neighbors of the training set to the test set. They then, in Theorem 1, analyze the asymptotic expected MSE, calibration and negative log-likelihood of this prediction approach under a given estimator of the parameters. They find that there is robustness to poor estimation except in noise variance. They show empirically that their simulations match theory in algorithm 1. They then propose a scalable GP regression algorithm: noting that when using their nearest neighbor prediction approach, the MSE, calibration and negative log-likelihood are *somewhat* robust to wrong parameters (except noise variance), they use multiple approximation steps, including using a subset of the training data and using a structured covariance matrix. In order to improve noise variance estimation, they add a calibration step. They show that it outperforms several baselines on real world datasets.", "questions": "Why did you use this particular, somewhat strange, paper structure?\n\nWhy is your calibration definition what it is?\n\nI'm a little confused: you said you replace $X$, but eqn. 5 still has it. Is this a typo?\n\nIs SVGP really state of the art? What about SKI [1] and its extensions?\n\n[1] Wilson, Andrew, and Hannes Nickisch. \"Kernel interpolation for scalable structured Gaussian processes (KISS-GP).\" International conference on machine learning. PMLR, 2015.", "limitations": "Yes.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "This is a very interesting paper: the resulting technique is very simple but pops out in a very non-obvious way. Further, it leads to a powerful consequence: you don't need the greatest parameter estimates for training if you do nearest neighbors prediction, as long as you add a calibration step to improve the estimate of noise variance. Despite my slight reservations about the baselines and the strange paper structure, I strongly recommend acceptance.", "weaknesses": "This paper is written in a very non-standard way and completely ignores the standard ML paper format. It completely ignores the standard introduction, replacing it by what is often the 2nd or third section, which is some technical background. It then has the 2nd section, which is *similar* to what often comes at the end of an intro, but is somehow different. The related work section is very short. In general everything is somewhat terse. I don't completely hate it, but it's jarring to read this, and I'm not sure what the justification for doing it this way is, particularly since the authors seem aware of literature and the techniques used in the literature so presumably are aware of how standard ML papers are written to flow a certain way.\n\nThe main weakness is that I'm not sure how close to state of the art the baselines are.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "1. Thank you for your comments on the structure. It was not a conscious decision to write substantially outside the stylistic norms of the community and on reflection perhaps we should have tried to conform more strongly. We do not wish to drastically change the structure of the paper at this stage in case that would lead to a requirement for other reviewers (who did not object) to reassess. We hope that you do not think this unreasonable and will definitely take your comments on board for future work. We are planning to elaborate and extend the related work section.\n\n2. We define calibration this way for a few reasons:\n\t- For a well-specified GP, the $E_Y[MSE]$ matches the predictive variance, i.e. the uncertainty in mean prediction (or the magnitude of the residuals) is reflected in the variance predicted by the model.\n\t- It provides a convenient numerical baseline for what \"well-calibrated\" means, when the value is 1.\n\t- It allows us to carry out our simple recalibration procedure (Algorithm 2), to improve both our measure of calibration and NLL, whilst leaving MSE unchanged.\n\t- Although a \"weak\" measure of calibration, obtaining a  value near 1 is a necessary condition for effective calibration, and marked departures from 1 were detected for some of the methods in figure 4 and table 3.\n\t\tAdditionally, we have subsequently found that [1] uses the same definition and we will reference that in the paper.\n[1] Jankowiak, Pleiss, and Gardner, ‘Parametric Gaussian Process Regressors’.\n\n3. This is not a typo, but we appreciate your pointing this notational ambiguity out so that we can rectify this. In a sense, under this model specifying $X$ and $N(x^*)$ are equivalent, since $N(x^*)$ is directly derived from the training set $X$.\n\n4. Given the potential for impact that the findings in this paper imply, we were keen to ensure that our benchmarks included widely used methods, so that an interested user might easily determine the potential for these findings to inform their work.  We are not aware of a method which has replaced SVGP in this position, in terms of generality of uptake and usage. We were aware of the paper on SKI that you cited, but were under the impression it had limited applicability to datasets with more than a handful of dimensions. We would like to thank you for exposing us to the various extensions however, and would aim to include them as additional baselines to compare against in future publications.  We do not believe that the omission of this comparison harms the conclusions which we draw in this paper, however."}, {"review_id": "LV9XW4BEQf", "reviewer": "Reviewer_N8D3", "summary": "The paper proposes a change in perspective in how Gaussian process regression can be leveraged by conditioning predictive distribution on only the neighboring data points. The paper argues against the common practice of using one singular set of data points for model hyperparameter estimation and prediction, and presents a framework for running scalable approximation of GP regression using locality information. ", "questions": "None. ", "limitations": "I have listed the limitations in the \"weakness\" section. ", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "I believe that the overall idea of leveraging locality in Gaussian process models has a lot of potential: while I am skeptical as the arguments in the paper seem unconvincing, I think there is merit in good empirical performance for such algorithms. However, I believe that a comparison with existing methods should be done in a more proper manner. ", "weaknesses": "While using KNN in combination of GP regression might seem like a straightforward idea to massively speed up the inference of GP models, I believe that many key aspects of GPs are overlooked, constituting a major weakness from this paper. I list some of the specific examples below. \n- Combining KNN with GP is overall not a novel idea, the practice is not exactly new. A cursory search reveals that many research papers outside the machine learning community are already using it as a proxy for the less scalable fullscale GP regression. Some examples relevant to machine learning include: 1. Wu L, Pleiss G, Cunningham JP. Variational nearest neighbor Gaussian process. In: Proceedings of the 39th International Conference on Machine Learning, PMLR; 2022; 2. Chen H, Zheng L, Kontar RA, Raskutti G. Gaussian Process Parameter Estimation Using Mini-batch Stochastic Gradient Descent: Convergence Guarantees and Empirical Benefits. Journal of Machine Learning Research. 2022;23(227):1–59. \n- The paper focuses on _pointwise_ uncertainty prediction, while GP posteriors provide not only a pointwise uncertainty, but _pairwise_ covariance as well. The covariance kernel is such a crucial element that defines the hypothesis space of GP models that given a configuration of pointwise uncertainty, there exists countless GPs that takes identical pointwise variance. \n- While KNN is guaranteed to speed up massively in lower dimensions, it suffers from curse of dimensionality itself as it becomes harder to find nearest neighbors. \n- While it is not mandatory for parameter estimation and prediction to be conditioned on the same set of data, SVGP does not incur much computational expense at prediction level. Conventional GP regression and SVGP only need one matrix inversion operation to predict the uncertainty level on an arbitrary number of test points, while GPnn needs a matrix inversion operation for every separate data point -- essentially it becomes a balancing act of whether optimizing an SVGP or apply KNN to every test point is more expensive. \n- Theorem 1 discusses the asymptotic behavior as $n\\rightarrow\\infty$. While the result of insensitivity w.r.t. kernel parameter carries some value, the _rate_ of convergence matters in a finite data regime. For example, Ackermann function and inverse Ackermann function both converge to infinity but at drastically different rates. \n- While the common squared exponential kernel is designed to look at neighborhood information, many more expressive kernels can have longer-range connections. I think the paper tacitly admits that GPnn can only be used for interpolation tasks, but it is somewhat unfair to state that kernel hyperparameter does not matter as it is always a misspecified model. \n- The ease at obtaining easy uncertainty information given neighborhood information is conditioned on conjugate likelihood, with a Gaussian observational model. As some degree of variational inference is required for non-conjugate likelihood, GPnn struggles to apply beyond regression tasks. \n- Simulating synthetic data using Algorithm 11 seems flawed. Algorithm 1 generates test points as if every 2 test points are un-correlated, which naturally gives GPnn an advantage. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "1. We briefly mention use of NN in other GP literature in sec. 8 \"related work\" but appreciate your feedback pointing to the need to expand that further, including how our use differs from theirs. We agree that  the mere act of including kNN in some form within a paper cannot be regarded as novel. However, the novelty is in how we use it, the resulting strength of results and the substantial computational efficiency savings made. Novel aspects include:  robustness in the limit to kernel choice and hyperparameters, substantial reductions in training costs, justification of detachment of prediction process from training process, improvements in calibration over other methods, algorithmic simplicity. These innovations are more general than use of NN alone.\n\n2. Thank you for this perspective. We have focused on pointwise prediction and, in light of your feedback, will clarify that in the paper. Whilst accepting pointwise prediction does not cover the entire spectrum of GP applications it is an exceptionally important area in its own right (and the most heavily used in practice) for which nearly all research papers quote their performance results, e.g. [2,4,7,11,12,13,22,24,25]. The importance of making improvements in this area is also reflected in the comments of all reviewers of this paper. Selling points of GPs in the pointwise context are provision of accuracy and well principled uncertainty measures. Our method is shown to generally outperform the other methods in both respects at a small fraction of their training cost.\n\n3. Prior to our analysis, we too had thought that curse of dimensionality would impact badly on comparative performance for high-dimensional datasets. Surprisingly, we found mse, nll and calibration still beat other methods at large d (e.g. Ctslice with d=378). This led to conjecture 6 and related follow-on work in progress. Computational cost does rise with d, as covered in 6.2, table 2 and figure 3, but even for the most extreme case (d=378) GPnn outperforms other methods on training time and is comparable at test time.  Given that the precise implementation of the approximate kNN algorithm is beyond the scope of this paper, and should not negatively affect the locality arguments on which this work is structured, we believe that advances in this area will only improve the performance beyond its already strong baseline.\n\n4. Whilst true that SVGP prediction is fast, GPnn prediction is comparable (see table 1 in the pdf) and training time much faster than that of SVGP, meaning for a given computational budget covering both train and test time GPnn is very competitive, especially at large n. On test point timings, table 1 gives exemplar prediction times in seconds, all obtained on a laptop (with 400 NNs and 1024 inducing points as in the paper). We will add similar results to the paper to reassure readers on this point.\n\n5. We agree that the rate of convergence is of practical value and we have preliminary results for this but wish to delay publication of them (see 5.1). Nevertheless we believe the importance of sharing existing results with the community outweighs delaying publication for completeness. This view appears to be supported by the other reviewers. \n\n6. GPnn _is_ designed for interpolation tasks on which we demonstrate excellent performance, even in the presence of somewhat severe misspecification, at modest expense. It is ongoing work to generalize our existing results and determine the dependence of the rate on the kernel properties.\nPractitioners generally fix a choice of kernel(s) and then optimize hyperparameters in relation to that choice. A key point that our paper makes is to caution against excessive effort in hyperparameter fine-tuning. Our empirical and theoretical results demonstrate that the predictive measures commonly used by the community are robust to non-optimal parameter specification for pre-picked kernels, which runs somewhat counter to conventional belief.\n\n7. As covered under bullet point 2, improving GP regression is of major importance in its own right so we did not view a failure to address non-regression applications, e.g. classification, to be a weakness of the paper. In fact, we note that the use of GP for regression tasks is more common than for classification, and thus we would argue that our findings, although focussed in the regression domain, have wide reaching and significant impact to the community. Having said that, we take on board your comments on the potential advantages of variational methods beyond regression whilst also remaining open minded as to whether GPnn could potentially play a role in that domain. \n\n8. \n- Re \"Algorithm 11 seems flawed\": Thank you for pointing out the need to justify this. We are confident the simulation is both valid and appropriate for what we are aiming to do and hope the following brief explanation helps: Firstly, we are focussing on point prediction and lack of correlation between test points does not impact that. Secondly, the validity of the estimates is provable: The gist of the argument is to show for fixed size-$n$ training set $X$ that $\\\\{(x_i^*,N(x_i^*),y_i^*,y_i)\\\\}\\_{i=1}^{n^*}$ is a set of iid vector-RVs (here the $y_i$ are m-dimensional). Then since $e_i^*$ is a deterministic function of $(x_i^*,N(x_i^*),y_i^*,y_i)$, $\\\\{e_i\\\\}\\_{i=1}^{n^*}$ is also a set of iid RVs. Finally $e^* =  E(e_i)$ so that $e^* = E[\\frac{1}{n^*} \\sum_{i=1}^{n^*} e_i]$ as required. The same argument holds for $l^*$ and $z^*$. \n- Re \"gives GPnn an advantage\": We do not believe this to be the case because (a) the estimated mean mse, nll and calibration are valid for (point-estimate) GPnn and (b) we are in any case only using this simulation to experimentally demonstrate convergence toward the theoretical limits of theorem 1, not to compare GPnn performance with the other methods."}, {"review_id": "gTgurXz8th", "reviewer": "Reviewer_3iYQ", "summary": "The paper proposes a Gaussian process regression technique that only uses nearest neighbors for prediction. The question of choosing hyperparameters is addressed and asymptotic behavior is analyzed. The method shows substantial speed up on UCI datasets while also providing improved predictive performance.", "questions": "Due to time constraints and NeurIPS not accommodating my request for a reduced reviewing workload, I did not have the opportunity to thoroughly review the proofs or examine all the details.", "limitations": "The limitations are addressed.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The empirical evaluation shows very good results. Even so, the method is theoretically analyzed and discussed with interesting insights and hypotheses. The paper is well-written.", "weaknesses": "A comparison with other nearest neighbor based techniques might have been helpful (e.g. ones cited in the paper). Also, it might make sense to investigate more into the geostatistical work on the subject and cite something from the area. I know that nearest neighbor based predictions are quite often considered to be standard in the area.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "1. Thank you for what time you had available and your comments. We take your point regarding references to other work on similar topics and intend to expand our section on \"related work\" to reflect this, including contributions from the geospatial community. \n2. We acknowledge that a comparison to other nearest-neighbour based methods might have been interesting, but we neglected to do so on the basis that the thrust of this paper was to emphasise, for example, the following points:\n- The decoupling of parameter estimation and prediction and the associated advantages to computational cost this can incur\n- The asymptotic insensitivity of GPnn prediction to model misspecification making use of nearest-neighbours to do so, rather than the focus being primarily on NNs itself.\n3. Please also see our response to reviewer MqBs bullet point 4 concerning our choice of \"state of the art\" comparisons."}]

Author rebuttal:
We thank all of the reviewers for their useful and constructive feedback and for the time and effort that they have voluntarily set aside for this task.

Meta-review: This paper investigate nearest neighbor approximations to GP posteriors and provides a theoretical asymptotic analysis. The proposed method offers many advantages over SVGP (a scalable GP approximation used extensively by the ML community), which the authors demonstrate experimentally. Moreover, despite its asymptotic limitations, the theoretical analysis is sounds and provides interesting insights. Based on these strengths, I recommend this paper is accepted into NeurIPS.

However, there is a glaring weakness off this paper, which I would strongly encourage the authors to address in the revisions: the treatment of prior NN-approximation work from the GP literature. While the authors nod to a few recent works, they do not include any discussion of them, and how these approaches differ from what is proposed. Most notably, the authors have also missed the pioneering work of Vecchia (1988) and many of the follow-up works - a method which has much popularity in the geostatistics community. The theoretical analysis of the authors is an important contribution; however, I fear that this analysis might not actually be too useful since the NN approximation analyzed by the authors differs significantly from these well-established/highly used approximations.

I would encourage the authors to:
- Include more relevant (and popular) NN-based GP approximations, on top of the more recent ones in the related work section
- Provide a substantive description of these methods and how the proposed method differs. (Would the same theoretical analysis apply to these other methods?)
- Include an experimental comparison of the proposed method against other NN-based methods.

Vecchia, A. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society, Series B, 50(2):297–312.
Decision: Accept (poster)