Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent
Authors: Jihao Andreas Lin
Javier Antoran
Shreyas Padhy
David Janz
José Miguel Hernández-Lobato
Alexander Terenin
Venue: NeurIPS 2023 oral (2023)
Keywords: Gaussian processes
scalable learning
posterior sampling
Bayesian optimization
TL;DR: We sample from GP posteriors using SGD and develop a spectral characterization for why it works, even in cases of non-convergence.

Abstract: Gaussian processes are a powerful framework for quantifying uncertainty and for sequential decision-making but are limited by the requirement of solving linear systems. In general, this has a cubic cost in dataset size and is sensitive to conditioning. We explore stochastic gradient algorithms as a computationally efficient method of approximately solving these linear systems: we develop low-variance optimization objectives for sampling from the posterior and extend these to inducing points. Counterintuitively, stochastic gradient descent often produces accurate predictions, even in cases where it does not converge quickly to the optimum. We explain this through a spectral characterization of the implicit bias from non-convergence. We show that stochastic gradient descent produces predictive distributions close to the true posterior both in regions with sufficient data coverage, and in regions sufficiently far away from the data. Experimentally, stochastic gradient descent achieves state-of-the-art performance on sufficiently large-scale or ill-conditioned regression tasks. Its uncertainty estimates match the performance of significantly more expensive baselines on a large-scale Bayesian~optimization~task.

Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent
Jihao Andreas Lin Javier Antorán"!
Shreyas Padhy" Alexander Terenin 1.4
David Janz
José Miguel Hernández-Lobato"
University of Cambridge University of Alberta
Max Planck Institute for Intelligent Systems
*Cornell University
Abstract
Gaussian processes are powerful framework for quantifying uncertainty and for sequential decision-making but are limited by the requirement of solving linear sys- tems. In general, this has a cubic cost dataset size and is sensitive conditioning. We explore stochastic gradient algorithms as a computationally efficient method of approximately solving these linear systems: we develop low-variance optimization objectives for sampling from the posterior and extend these to inducing points. Counterintuitively, stochastic gradient descent often produces accurate predictions, even in cases where it does not converge quickly to the optimum. We explain this through a spectral characterization of the implicit bias from mon-convergence. We show that stochastic gradient descent produces predictive distributions close to the true posterior both regions with sufficient data coverage, and in regions sufficiently far away from the data. Experimentally, stochastic gradient descent achieves state-of-the-ar performance on sufficiently large-scale or ill conddiiiioned regression tasks. Its uncertainty estimates match the performance of significantly more expensive baselines on a large-scale Bayesian optimization task.
Introduction
Gaussian processes (GPs) provide comprehensive framework for learning unknown functions in an uncertainty-awe aware manner. This often makes Gaussian processes the model choice for sequential decision-making achieving Kate-offthhe-rr performance in tasks such optimizing molecules in computational chemistry settings 21 and automated hyperparameter tuning 39 24
The main limitations of Gaussian processes is that their computational cost is cubic in the training dataset size. Significant research efforts have been directed at addressing this limitation resulting in two key classes of scalable inference methods: (i) inducing point methods [33 23 which approximate
19 which approximate the computation needed obtain the GP posterior. Note that in structured settings, such geospatial lcarning in low dimensions, specialized techniques are available (49 48 [5]. Throughout this work, we focus on the generic setting, where scalability limitations are of yet unresolved.
GP posterior, and (ii) conjugate gradient methods [20
In recent years, stochastic gradient descent (SGD) emerged the leading technique for training machine learning models at scale (37, 42 in both deep learning and related settings like kernel methods [14) and Bayesian modeling 29 While the principles behind the effectiveness of SGD are not yet fully understood, empirically, SGD often leads to good predictive performance- when it does not fully converge. The latter is the default regime in deep learning, and has motivated researchers to study implicit biases and related properties of SGD
52
"Equal contribution, order chosen randomly. Code available at: HTTPS://GITHUB COw.///C/////CCCCCCCCC MLG/SGD-OO
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

In the context of GPs, SGD is commonly used to learn kernel hyperparameters by optimizing the marginal likelihood [8]
L
19
12 IT or closely related variational objectives [43 23 In
this work, we explore applying SGD to the complementary problem of approximating GP posterior samples given fixed kernel hyperparameters.
one of his seminal books on statistical learning theory, Vladimir Vapnik [45] famously said: "When solving a given problem, try to avoid solving a more general problem as an intermediate step." Motivated by this viewpoint, as well as the aforementioned property of good performance often not requiring full convergence when using SGD, we ask: Do the linear systems arising in GP computations necessarily need to be solved to a small error tolerance? not, can SGD help accelerate these computations?
We answer the latter question affirmatively. with specific contributions as follows: (i) We develop a scheme for drawing GP posterior samples by applying SGD to quadratic problem. In particular,
we re-cast the pathwise conditioning technique of Wilson et al. (50, 51] as an optimization problem to which we apply the low-variance SGD sampling estimator of Antorán et al. [3]. We extend the proposed method inducing point Gaussian processes. (ii) We characterize the implicit bias in <GD-approximated GP posteriors showing that despite optimization not fully converging, these match the true posterior in regions both near and far away from the data. (iii) Experimentally. we show that SGD produces strong cesupas-compared to variational and conjugate gradient methods on both large-scale and condiy-conditoone regression tasks, and on a parallel Thompson sampling task, where error bar calibration is paramount.
2 Optimization and Pathwise Conditioning in Gaussian Processes
random function f: X + R over some set X is a Gaussian process 35| if. for every finite set
of points x: E X.N, f(x=) is multivariate Gaussian. A Gaussian process GP k) is uniquely
determined by a mean function p(-) = E-)(()) and a covariance kernel k(,, ()) = [C)()())
We denote by K the kernel matrix [k(I2,, (r;)) J...=.....N We consider the Bayesian model y f() + where E ~ N(0, Σ) gives the likelihood, f f~ GP k) is the prior, and H, are the training data. The posterior of this model is fl y (P k(l)).
with
(-) = (()) (Kaz + Σ) 1
(.) = ((,)) K(-ja (Kaz + )) ->>"Ka(>)> (I)
Using pathwise conditioning 50
51
can also write the posterior directly the random function
(/ v-) - ʃ(-) + K()) K
Σ) (y f(x) +) ~~ N(O Σ) ʃ ~ G). (2)
We assume throughout that Σ is diagonal. Due to the matrix inverses present, the cost to directly compute each of the above expressions is ((N3).
2.1 Random Fourier Features and Efficient Sampling
Our techniques will rely on random Fourier features 34 40
Let X Rd and let k be stationary-
that is, k(æ,, x²) = k(x - x²). Assume, in this section only and without loss of generality, that k(x,x)) 1. Random Fourier features are sets of L, random functions + : X RR whose
- (2 for odd, and Φr(-) L-1/2 (2r :)) for É even. By taking We to be samples distributed according p. the nor- malized spectral measure of the kernel, for any x, x' E X we have
components, indexed by é < L, are (e(-)
k(2,x') = E.w..... u ( (x)) ((x'))
which we will use to construct unbiased estimators of kernel matrices in the sequel. Fourier features can also be used to efficiently sample functions from the GP prior. Naively, evaluating the random function f at N. points requires computing a matrix square root of K at O(N3) computational cost. However, given a precomputed set of Fourier features &, the Monte Carlo estimator
f() (() = 0T (()
Σ (.
0 N(O,I)
approximately samples from GP k) at ON.) cost. The approximation error in 4 decays as goes to infinity. Following Wilson et al. (50) 51), our proposed algorithms will approximately sample from the posterior by replacing if with in the pathwise conditioning formula 2

---

Stochastic Gradient Descent
Conjugate Gradients
Variational Inference
5.0
Infill asymmtotics
2.5
0.0
-.5
5.0
dorgen asymntotics
0-
-2
3.5
-7
-3.5
3.5
7
-3.5
3.5
Exact GP
Approximations
Figure 1: Comparison of SGD, CG [46] and SVGP (23 for GP inference with squared exponential
kernel on 10k datapoints from sin (2x + cos(5x) with observation noise N(0, 0.5). We draw 2000 function samples with all methods by running them for 10 minutes on an RTX 2070 GPU. Infill asymptotics considers z; ~ N(0, 1). A large number of points near zero result in very ill-conditioned kernel matrix, preventing CG from converging. SGD converges in all of input space except at the edges of the data. SVGP can summarise the data with only 20 inducing points. Note that CG converges to the exact solution if one uses more compute, but produces significant errors stopped too early, occurs under the given compute budget. Large domain asymptotics considers data regular grid with fixed spacing. This problem is better conditioned, allowing SGD and CG to recover the exact solution. However, 1024 inducing points not enough for SVGP summarize the data.
2.2 Optimization-baset Learning in Gaussian Processes
Following Matthews et al. 31 and Antorán et al. (3). both a Gaussian process' posterior mean and posterior samples can be expressed as solutions quadratic optimization problems. Letting 2 = (Ka + E) in (T.. we can express the GP posterior mean in terms of the quadratic problem
N
(yi U)
(-) = K (()w''
Σ (
O' arg min BERN
+ ||
(5)
En
We say that k((x4,)) :) are the canonical basis functions, are the representer weights || I = - K the regularizer. The respective optimization problem for obtaining posterior
38 and that
samples is similar, but involves a stochastic objective, which will be given and analyzed in Section 3.2
Conjugate gradients (CG) are the most widely used algorithm to solve quadratic problems, both
in the context of GPs 22 19 4 and more generally [33,[0] 10]. CG reduces GP posterior inference to series of matrix-vecte vector products, each of O(N2) cost. Given the system A -b,, the number of matrix-vector products needed to guarantee convergence of CG to within a tolerance of e [41. is
of (econd( (conddA)
cond(A) Amaze ( (A)
(6)
Auiin (A)
where Amaxx (A) and Amii (A) are the maximum and minimum eigenvalues of A. CG performs well in many GP use cases, for instance Gardner et al. [19 and Wang et al. [46). In general, however, the condition number cond(K + E) need not be bounded, and conjugate gradients may fail to converge quickly [41]. Nonetheless, by exploiting the quadratic structure of the objective, CG obtains substantially better convergence rates than gradient descent 8 52 This presents a pessimistic outlook for gradient descent and related methods such as SGD which do not take explicit advantage of the quadratic nature of the objective. Nonetheless, the success of SGD in large-scale machine learning (9] motivates us to try it anyway.
3

---

Gaussian Process Predictions via Stochastic Gradient Descent
We now develop and analyze techniques for drawing samples from GP posteriors using SGD. This is done by rewriting the pathwise conditioning formula 22 in terms of two stochastic optimization problems. As a preview of what this will produce, we showcase SGD's performance on pair of toy problems, designed to capture complementary computational difficulties, in Figure
3.1 A Stochastic Objective for Computing the Posterior Mean
We begin by deriving stochastic objective for the posterior mean. The optimization problem ( with optimal representer weight solution v", requires O(N²) operations to compute both its square eπrоr and regularizer terms exactly. The square error loss term is amenable to minibatching, which gives an unbiased estimate in (NV) operations. Assuming that k is stationary, following Section 2.1 we can stochastically estimate the regularizer with Fourier features using the identity (2:) Combining both estimators gives our SGD objective
|| K..
E Bw 8 $ )
N D
K 0) Σ
L
Σ [=1
(o" )e(x))))
t=1
where D is the minibatch size and L the number of Fourier features. This regularizer estimate is unbiased even when drawing a single Fourier feature per step: the number of features controls the variance. Equation (7) presents (NN complexity, in contrast with the O(N2) complexity of one CG step. We discuss sublinear inducing point techniques further on, but first turn to sampling.
3.2 Computing Posterior Samples
We now frame GP posterior samples in a manner amenable to SGD computation similarly to (7) First, we re-write the pathwise conditioning expression given in (2) as
(f (.))))
f(++
K
(() (K (Kaz
(f)=)
E g~~ N(O,E)
f f~ k). (8)
prior
postetior mean
uncertainty reductioo term
We approximate the prior function sample using the Fourier feature approach of (4) and the posterior mean, defined in Section 2.2 with the minimizer v of (7) obtained by SGD. Each posterior sample's uncertainty reduction term is parametrized by a set of representer weights given by a linear solve against noisy prior sample evaluated at the observed inputs, namely a (Kaz + Σ) (f(x)) +e). We construct an optimization objective targeting a sample's optimal representer weights as
(f(xz) + - K a))
f(x) ~ N(O, Kmn) N(O,E) E).
a arg min
||a:| K
(9)
QERN
Ei
Applying minibatch estimation to this objective results in high gradient variance, since the presence of e, makes the targets noisy. To avoid this while targeting the same objective, we modify ( as
(f(xs) - K a)2
f(x) N(O, K N(O,
a arg min
+ |la -
(10)
Σg
GERN
moving the noise term into the regularizer. This modification preserves the optimal representer weights since objective (10) equals (9) up to constant: a proof is given in Appendix D This
generalizes the variance reduction technique of Antorán et al. (3] to the GP setting. Figure
illustrates minibatch gradient variance for these objectives. Applying the minibatch and random feature estimators of 7, obtain per-step of O(NS), for S the number of posterior samples.
3.3 Inducing Points
So far, our sampling objectives have presented linear cost in dataset size. In the large-scale setting, algorithms with costs independent of the dataset size are often preferable. For GPs, this can be achieved through inducing point posteriors 43 23 to which now extend SGD sampling.

---

Gradient Variance
RMSE
NLL
Total Runtime
6k
Loss
0.75
-0.6
10h
4k
Loss 2
0.50
-0.8
2k
0.25
1.0
0.00
-12
0.1h
10
10
10a
10
10
10a
10
10
10a
10a
10
Iterations
Time in sec
Time in sec
Standard SGD
Inducing SGD
Inducing SGD
Inducing SGD
1844352
099 206
M 431 489
M 218 782
Figure 2: Left: gradient variance throughout optimization for a single-sample minibatch estimator (D of (9). labeled Loss and the proposed sampling objective (TO), labeled Loss 2, on the ELEVATORS dataset (N = 16k). Middle plots: test RMSE and negative log-likelihood (NLL) obtained by SGD and its inducing point variants, for decreasing numbers of inducing points, given in the rightmost plot, function of time on an A100 GPU, on the HOUSEELECTRIC dataset (N x 2M).
Let E XM be set of M E N inducing points, Applying pathwise conditioning to the Kullback- Leibler-optimal inducing point approximation of Titsias [43) gives the expression
(r(=) (-) (() " (-) (()= Ky ( K >)- ())(((((((x)) +e)
(11)
N(O, 2)
f ~ G))
(()) - K fz).
Following Wild et al.
(7 Theorem 5, the optimal inducing point mean / can therefore be written
(yi U,,
1 K 0)
(()))))))
arg min VERM
+ |
(12)
and we can again parameterize the uncertainty reduction term as K (()=a' with
a arg min (xx) + 022 + ||a||
(()) N(O, K= K Kzn) E N=0, E).
(13)
Σx
OERM
A full derivation is given in Appendix C Exact implementation of (13 is precluded by the need to sample from a Gaussian with covariance K2 Kgz- However, identify this matrix Nyström approximation to KE Thus, can approximate (13) by replacing f(x) with f. which can be sampled with Fourier features (4). The approximation error is small when M is large and the inducing points close enough to the data. Our experiments in Appendix [[show it to be negligible in practice. With this, we apply the stochastic estimator (7) to the inducing point sampling objective.
The inducing point objectives differ from those presented in previous sections in that there are ((M) and not (NN learnable parameters, and we may choose the value of M and locations freely. The cost of inducing point representer weight updates is thus O(SM), where S is the number of samples. This contrasts with the (M*)) update cost of stochastic gradient variational Gaussian processes 23 Figure 2 shows that SGD with M x 100k inducing points matches the performance of regular SGD on HOUSEELECTRIC (N Rs 2M), but an order of magnitude faster.
3.4 The Implicit Bias and Posterior Geometry of Stochastic Gradient Descent
We have detailed an SGD-based scheme for obtaining approximate samples from posterior Gaussian process. Despite SGD's s significantly lower cost per-iteration than CG. its convergence to the true optima, shown in Figure 3 is much slower with respect to both Euclidean representer weight space, and the reproducing kernel Hilbert space (RKHS) induced by the kernel. Despite this, the predictions obtained by SGD are very close to those the exact GP, and effectively achieve same test RMSE, Moreover, Figure shows the SGD posterior width close to the data, and which revert smoothly to the prior far away from the data. Empirically, differences between the SGD and exact posteriors concentrate at the borders of data-dense regions.
ID toy task exhibits error bars of the correct
5

---

RMSE to Test Data
RMSE to Exact GP
Enclidean Error
RKHS Error
3.6
0.9
Low Opt.
noise nolse
3.5
0.6
2.4
0.6
0.0
0.0
1.2
0.3
160
5
0.0
0.0
30
10
10
10
10a
10a
10
10-$10"
10"
10 10a
10
'10"
10'
10
10a
10
1001100 10' 102 10a
Time in sec
Time in sec
Time in sec
Time in sec
SGD (optimal noise)
SGD (low noise)
CG (optimal noise)
CG (low noise)
Figure 3: Convergence of GP posterior mean with SGD and CG as function of time (on A 100 GPU) on the ELEEVTOOSS dataset (N R 16k) while setting the noise scale to (i) maximize exact GP marginal likelihood and (ii) to 10-3, labeled low noise. We plot, in left-to-right order, test RMSE, RMSE to the exact GP mean at the test inputs, representer weight error |/ v" ll2, and RKHS error || - MSGD H In the latter two plots, the low-noise setting is shown on the bottom.
We now argue the behavior seen in Figure general feature of SGD: one can expect it obtain good performance even in situations where it does not converge to the exact solution. Consider posterior function samples pathwise form, namely (f (.) (())) where ~ GP k) is a prior function sample and are the learnable representer weights. We characterize the behavior of SGD-computed approximate posteriors by splitting the input space X into 3 regions, which call the prior, interpolation, and extrapolation regions. This is done as follows.
(I) The Prior Region. This corresponds to points sufficiently distant from the observed data. Here, for kernels that decay over space, the canonical basis functions k(z,,)) -) go to zero. Thus, both the true posterior and any approximations formulated pathwise revert to the prior. More precisely, let X = =dd
_ let k satisfy lime--oe k(x', c. x) = O for all x' and x in X, and let (f y)) be given by (f v) f()) + (),u, with IR Then, by passing the limit through the sum, for any fixed
it follows immediately that lim (ʃ (c x) = flc æ). Therefore, SGD cannot incur error in regions which are sufficiently far away from the data. This effect is depicted in Figure 4
(II) The Interpolation Region. This includes points close to the training data. We characterize this region through subspaces of the RKHS, where we show SGD incurs small error.
Let K az = UAUT be the eigendecomposition of the kernel matrix. We index the eigenvalues A diag À1. ..... AN) in descending order. Define the spectral basis functions eigenvector- weighed linear combinations of canonical basis functions
N
Us k(I,.)) u (0)) = Σ
(14)
These functions which also appear in kernel principal component analysis are orthonormal with respect to the RKHS inner product. To characterize them further, we lift the Courant-Fischer characterization of eigenvalues and eigenvectors to the RKHS H₃ induced by k. obtaining the expression
(0 arg max NeHt
Σ u(xx)) || (u, 20)) H. = 0, <i
(15)
This means in particular that the top spectral basis function, u (.)) is a function of fixed RKHS norm that is, of fixed degree of smoothness, as defined by the kernel k-which takes maximal values at the observations -n. Thus, (1) will be large near clusters of observations. The same will be true for the subsequent spectral basis functions, which also take maximal values at the observations, but are constrained to be RKHS-orthogonal proof of the preceding expression is given in Appendix E.3] Figure confirms that the top spectral basis functions are indeed centered on the observed data.
previous spectral basis functions. A
Empirically, SGD matches the true posterior in this region. We formalize this observation by showing that SGD converges quickly in the directions spanned by spectral basis functions with large eigenvalues. Let projuco (-) be the orthogonal projection onto the subspace spanned by (9)

---

SGD Posterior
Wasserstein-2 to Exact GP
W2
2.5
-5.0
-10
-5
10
Kernel Matrix Eigenvalues
Spectral Basis Functions
10c
0.15
MA M
16
0.07
10a
41
45
W
10-4 44
-0.07
10-9
-0.15
10
10a 10
-10
10
Index
Figure 4: SGD error and spectral basis functions. Top-left: SGD (blue) and exact GP (black, dashed) fit to a N = 10k, d = 1 toy regression dataset. Top-right: 2- 2-Wasserstein distance (W2) between both processes' marginals. The W2 values are low near the data (interpolation region) and far away from the training data. The error concentrates at the edges of the data (extrapolation region). Bottom: The low-index spectral basis functions lie on the interpolation region, where the W2 error is low, while functions of index 10 and larger lie on the extrapolation region where the error is large.
Proposition 1. Let 8' = for o² 0. Let PSGD be the predictive obtained by Polyak-averaged SGI) after steps, starting from an initial set of representer weights equal to zero, and using sufficiently small learning rale of 0- < the gradient is E-sub-Gaussion Then, with probability 1 ō, we have for 1,... that
o² ('( Assume the stochastic estimate of
1
1 n +
2 t
projum Pfl - projuc PSGD
I
G.
log
(16)
V
6
The proof, as well as an additional pointwise convergence bound and a variant that handles projections onto general subspaces spanned by basis functions, are provided in Appendix E. In general, expect G to be at most (>f bllll 20) with high probability.
The result extends immediately from the posterior mean to posterior samples. As consequence, SGD converges to the posterior GP quickly in the data-dense region, namely where the spectral basis functions corresponding large eigenvalues are located. Since convergence speed the span of each basis function is independent of the magnitude of the other basis functions' eigenvalues, SGD perform well even when the kernel matrix is ill-conditioned. This is shown in Figure
(III) The Extrapolation Region. This can be found by elimination from the input space of the prior and interpolation regions, in both of which SGD incurs low error. Consider the spectral basis functions (.) with small cigenvalues. By orthogonality of (.."),...!! be large near the observations while retaining prescribed norm. Their mass is therefore placed away from the observations. SGD converges slowly in this region, resulting in large error in its solution in both Euclidean and RKHS sense, seen in Figure 3 Fortunately, due the lack of data in the extrapolation region, the excess test error incurred due to SGD nonconvergence may be low, resulting benign nonconvergence (52). Similar phenomena have been observed in the inverse problems literature, where this is called derative regularization (22 27). Figure 4 shows the Wasserstein distance to the exact GP predictions is high in this region, as when initialized zero SGD tends return small representer weights, thereby reverting the prior.
(!,,,, such functions cannot
7

---

pol
keggdirected
3droad
houseelectric
N 15 000
N 48 827
NV
434 874
N
2 049 280
NLL RMSE
1.2
2.4
1.2
1.6
0.6
1.2
0.6
0.8
0.0
0.0
0.0
0.8
1.1
0.7
10
-0.3
0.0
0.0
-1.4
-0.7
2
SGD
10
10
10a
10a
10
10a
10a
10a
10
102
10 10*
10a
10'
10
10
10"
CG
Time in sec
Time in sec
Time sec
Time in sec
Figure 5:
Test RMSE and NLL as a function of compute time on TPUv2 core for CG and SGD.
Experiments
We now turn to empirical evaluation of SGD GPs, focusing on their predictive and decision-making properties. We compare SGD GPs with the two most popular scalable Gaussian process techniques: preconditioned conjugate gradient (CG) optimization [9] 46) and sparse stochastic variational inference (SVGP) 43 23 For CG, we pivoted Cholesky preconditioner of size 100, except in cases where this slows down convergence, where we instead report results without preconditioning. We employ the GPJax (32) SVGP implementation and use M = 4096 inducing points for all datasets, initializing their locations with the K-means algorithm. Full experimental details are in Appendix
4.1 Regression Baselines
We first compare SGD-based predictions with baselines in terms of predictive performance, scaling of computational cost with problem size, and robustness to the ill-conditioning of linear systems. Following Wang et al. [46]. we consider 9 datasets from the UCI repository (16) ranging in size from N 15k to N 2M datapoints and dimensionality from d = 3 to d = 90. We report mean and standard deviation over five 90%-train 10%-test splits for the small and medium datasets, and three splits for the largest dataset. For all methods, we use a Matérn-3/2 kernel with a fixed set of hyperparameters obtained using maximum marginal likelihood (2). as described in Appendix A.
We run SGD for 100k steps, with a fixed batch size of 512 for both the mean function and samples. For CG, we run maximum of 1000 steps for datasets with N S 500k, and a tolerance of 0.01. On the four largest datasets, CG's per-step cost is too large to run 1000 steps: instead, we run 100 steps, which takes roughly 9 hours per function sample on a TPUv2 core. For SVGP, we learn the variational parameters for M 4096 inducing points by maximizing the ELBO with Adam until convergence. For all methods, we estimate predictive variances for log- likelihood computations from 64 function samples drawn using pathwise conditioning (2).
Predictive performance and scalability with the number of inputs. Our complete set of results provided in Table including test RMSE, test negative log-like likelihood (NLL) and compute time needed to obtain the predictive mean single core of TPUv2 device. Drawing multiple samples the small setting (N < 20k),
requires repeating this computation, which we perform in parallel. taking 100k steps of SGD presents compute cost comparable to running CG to tolerance, which usually takes around 500-800 steps. Here, CG converges to the exact solution, while SGD tends to present increased test еитог due non-convergence. In the large setting (N 100k), where neither method converges within the provided compute budget, SGD achieves better RMSE and NLL and results. SVGP always converges faster than CG and SGD, but it only performs best on the Buzz dataset, which can likely be summarised well by M = 4096 inducing points.
From Figure 5 we see that SGD makes the vast majority of its progress in prediction space in its first few iterations, improving roughly monotonically with the number of steps. Thus, early stopping after 100k iterations incurs only moderate errors. In contrast, CG's ss initial steps actually increase test error, resulting very poor performance stopped too early. This interacts poorly with the number of CG steps needed, and the per-step cost, which generally grow with increased amounts of data 41).
8

---

Table 1: Regression task mean and std-err for GP predictive mean RMSE, low-noise RMSE (t). TPUv2 node hours used to obtain the predictive mean, and negative-log log-liielihood (NLL) computed with variances estimated from 64 function samples. SVGP is omitted for the low noise setting, where it fails to run. Metrics are reported for the datasets normalized to zero mean and unit variance.
Dataset
POL
ILEVATORS
BIKB 17379
PROTEIN
KBGGDIR
3DROAD 434874
SONG
BUZZ
HOUSCELBC
N
15000
16599
45730
48827
515345
583250
2049280
RMSE
0.00
0.38 0.00
0.11 0.00
4: 0.00
0.12
0.00
0.00
0.80 0.00
0.42:0.001
0.09 0.00
0.00
0.00
0.50 0.00
0.08
0.08:0.00 0.00
0.05
0.85
1.41 0.08
0.87 = 0.14
SVGP
0.00
0.37 0.00
0.07 0.00
0.57 - 0.00
0.09
0 0.00
0.49
0.09
0.11
0.33 + 0.00
0.11 = 0.01
RMSE
SGD
0.13 0.00
0.38
0.00
0.11 = 0.00
0.51 +: 0.00
0.12 0.00
0.11 0.00
0.80
0.00
0.42 + 0.01
0.09 0.00
0.06
0.68
0.09
0.05
0.23
9.79 1.06
0.02
0.02 0.83
1:
0.93 0.19
SVGP
Minutes
SGD
0.01
3,51 0.01
5.70 0.02
7.10 +: 0.01
15.2 : 0.02
27.6 11.4
220 - 14.5
347 61.5
162 : 54.3
CG
2.18
0.32
1,72 -: 0.60
281
0.22
9.07
+: 1.68
12.5 :: 1.99
85.2 = 36.0
195
:: 2.31
351
48.3
157 0.41
SVGP
21.2.1..2
21.3
0.112
20.5
0.02
1: 0.04
0.05
24.7
0.05
0.05
20.0
0.03
SGD
.0.70 = 0.02
0.47
0.00
-0.48
0.08
0.64
1: 0.01
-0.62 +: 0.07
-0.60 +: 0.00
1.21 1: 0.00
0.83
+: 0.07
-1,09 J: 0.04
NLL
CG
0.01
0.38 = 0.00
-2.62
0.06
0.62
+: 0.01
.0.92
0.10
0.10
0.45
1.36
2.38
0.08
2.07 0.58
SVGP
.0.71 =: 0.01
0.43,1 0.00
-1.27 1: 0.02
0.86
1: 0.01
-0.70
1: 0.02
0.67
=: 0.02
1.22 0.00
0.25
1: 0.04
-0.89 1: 0.10
Robustness to kernel matrix ill-conditioning. GP models are known to be sensitive to kernel matrix conditioning. We explore how this affects the algorithms under consideration by fixing the noise variance low value of x² = 10 6 and running them on our regression datasets. Table | shows the performance of CG severely degrades on all datasets and, for SVGP, optimization diverges for all datasets. SGD's results remain assentiall-unnchapee the smallest kernel matrix eigenvalues substantially and these do not affect convergence for the top spectral basis functions. This mirrors results presented previously for the ELEVATORS dataset in Figure 3
This because the noise only changes
Regression with large numbers of inducing points. We demonstrate the inducing point variant of our method, presented in Section [.3] on HOUSEELECTRIC, our largest dataset (N = 2M). We select varying numbers of inducing points with a K K-neareet-neighhoo algorithm, described in Appendix [A.2] Figure shows shows the time required for 100k SGD steps scales roughly linearly with inducing points. It takes 68m for full SGD and 50m, 25m, and 17m for M 1099k, 728k, and 218k, respectively. Performance in terms of RMSE and NLL degrades less than 10% even when using 218k points.
4.2 Large-scale Parallel Thompson Sampling
A fundamental goal of scalable Gaussian processes is to produce uncertainty estimates useful for sequential decision making. Motivated by problems in large-scale recommender systems, where both the initial dataset and the total number of users queried are simultancously large 36 17]. we benchmark SGD on large-scale Bayesian optimization task.
We draw a target function from a GP prior g GP(0, k)) and optimize it on X = [0, 1]d using parallel Thompson sampling (25). That is, we choose Fazw = arg max x (f (-) for a set of posterior function samples drawn in parallel. We compute these samples using pathwise conditioning with each respective scalable GP method. For each function sample maximum, we evaluate !/ew = g nenww) + e with 9 N(0, 10-6) and add the pair (Inew; !new) ) to the training data. We use an acquisition batch size of 1000 samples, and maximize them with a multi-stare gradient descent-basee approach described in Appendix A.3 We set the search space dimensionality to d = 8, the largest considered by Wilson et al. [50), and initialize all methods with a dataset of 50k observations sampled uniformly at random from X. To eliminate model misspecification confounding, we use Matérn-3/2 kernel and consider length scales of (0.1, 0.2, O.3. 0.4, 0.5) for both the target function and our models. For each length scale, we repeat the experiment for 10 seeds.
In large-scale Bayesian optimization, training and posterior function optimization costs can become significant, and predictions may be needed on demand. For this reason, we consider two variants of our experiment with different levels of compute. In the small-compute setting, SGD run for 15k steps, SVGP is given M 1024 inducing points and 20k steps to fit the variational parameters,

---

Small Compute Budget
Large Compute Budget
6.0
Max value
6.5
4.5
4.0
60k
70k
80k
2h
3h
70k
80k
15h
Acquisitions
Time elapsed
Acquisitions
Time elapsed
Random Search
SGD
CG
SVGP
Figure 6: Maximum function values (mean and std. err.) obtained by Thompson sampling with our approximate inference methods function of acquisition steps and of compute time on an A100 GPU. The latter includes time spent drawing function samples and finding their maxima. All methods share a starting dataset of 50k points and we take 30 Thompson steps, acquiring 1000 points in each.
and CG is run for 10 steps. In the large-compute setting, all methods for 5 times as many optimization steps.
Mean results and standard errors, across length scales and seeds, are presented in Figure 6 We plot the maximum function value achieved by each method. In the small-compute setting, the ~1.5 (A 100 GPU) hours required for 30 Thompson steps with SVGP and SGD are dominated by the algorithm used to maximise the models' posterior samples. In contrast, CG takes roughly twice the time, requiring almost 3h of wall-clock time. Despite this, SGD makes the largest progress per acquisition step, finding a target function value that improves upon the initial training set maximum twice as much as the other inference methods. VI and CG perform comparably, with the latter providing slightly more progress both per acquisition step and unit of time. Despite their limited compute budget, all methods outperform random search by significant margin. In the large- ange-compute setting, all inference methods achieve a similar maximum target value by the end of the experiment. CG and SGD make similar progress per acquisition step but SGD is faster per unit of time. SVGP is slightly slower than CG per both. In summary, our results suggest that SGD can be an appealing uncertainty quantification technique for large-scale GP-based sequential decision making.
5
Conclusion
In this work, we explored using stochastic gradient algorithms to approximately compute Gaussian process posterior means and function samples. We derived optimization objectives with linear and sublinear cost for both. We showed that SGD can produce accurate predictions- eveenn in cases when it does not converge to an optimum. We developed a spectral characterization of the effects of non-convergence, showing that it manifests itself mainly through error in an extrapolation region located away-but not too far away-from the observations. We benchmarked SGD on regression tasks of various scales, achieving state-of-the-art performance for sufficiently large or ill-conditioned settings. On a Thompson sampling benchmark, where well-calibrated uncertainty is paramount, SGD matches the performance of more expensive baselines at fraction of the computational cost.
Acknowledgments
We are grateful to David R. Burt for suggesting a few tricks which helped us complete the development of our theory and thank Bruno Mlodozeniec for identifying and correcting small mistake in the proof of Appendix E. JAL and SP were supported by the University of Cambridge Harding Distinguished Postgraduate Scholars Programme. JA acknowledges support from Microsoft Research, through its PhD Scholarship Programme, and from the EPSRC. JMHL acknowledges support from Turing AI Fellowship under grant EP/V023756/L. AT was supported by Cornell University, jointly via the Center for Data Science for Enterprise and Society, the College of Engineering, and the Ann S. Bowers College of Computing and Information Science. This work has been performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (http::/wwww. wwwpppp.p.......m///o/../ funded by an EPSRC Tier-2 capital grant. This work was also supported with Cloud TPUs from Google's TPU Research Cloud (TRC).
10

---

References
[1] J. Antorán, R. Barbano, Leuschner, M. Hemáández-Lobato, and B. Jin. Uncertainty
Estimation for Computed Tomography with Linearised Deep Image Prior. arXiv:2203.00479,
2022. Cited on page 14
[2] J. Antorán, Janz, J. U. Allingham, E. A. Daxberger, R. Barbano, E. T. Nalisnick, and
J. M. Hernández-Lobato. Adapting the Linearised Laplace Model Evidence for Modern Deep
Learning. International Conference on Machine Learning. 2022. Cited on page [
[3] J. Antorán, S. Padhy, R. Barbano, E. T. Nalisnick, D. Janz, and M. Hermández-Lobato. Sampling-based inference for large linear models, with application to linearised Laplace. In International Conference on Learning Representations, 2023. Cited on pages 14 20
[4] A. Artemev, D. R. Burt, and M. van der Wilk. Tighter Bounds on the Log Marginal Likelihood of Gaussian Process Regression Using Conjugate Gradients. In International Conference on Machine Learning, 2021. Cited on pages 3
[5] R. Barbano, J. Leuschner. J. Antorán, B. Jin, and J. M. Hemández-Lobato. Bayesian Experi- mental Design for Computed Tomography with the Linearised Deep Image Prior. In ICML Workshop on Adaptive Experimental Design and Active Learning in the Real World, 2022.
Cited on page 14
[6] M. Belkin, D. Hsu, S. Ma, and S. Mandal. Reconciling modern machine-learinng practice and the classical bias- variance trade-off. Proceedings of the National Academy of Sciences, 2019.
Cited on [
[7]
M. Bishop. Pattern Recognition and Machine Learning. Springer, 2007. Cited on page Blanchard and N. Krämer. Optimal Learning Rates for Kernel Conjugate Gradient Regres- sion. In Advances in Neural Information Processing Systems, 2010. Cited on page 3
[8]
[9] L. Bottou. Large-scale Machine Learning with Stochastic Gradient Descent. In International
Conference on Computational Statistics, 2010. Cited on page
[10] S. P. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2014.
Cited page
[11] H. Chen, I.. Zheng, R. Al Kontar, and G. Raskutti. Gaussian Process Parameter Estimation Us- ing Mini-batch Stochastic Gradient Descent: Convergence Guarantees and Empirical Benefits. Journal of Machine Learning Research, 2022. Cited page
[12] H. Chen, L. Zheng, R. Al Kontar, and G. Raskutti. Stochastic gradient descent in correlated settings: A study on Gaussian processes. Advances in Neural Information Processing Systems, 2020. Cited on page
[13] K. Cutajar, M. Osborne, J. Cunningham, and M. Filippone. Preconditioning Kernel Matrices. In International Conference on Machine Learning, pages 2529-2538, 2016. Cited on page [14] B. Dai, B. Xie, N. He, Y. Liang, A. Raj, M.-F. F. Balcan, and L. Song. Scalable Kernel Methods via Doubly Stochastic Gradients. Advances Neural Information Processing Systems, 2014. Cited [
[15] E. A. Daxberger, E. T. Nalisnick, J. U. Allingham, J. Antorán, and J. M. Hermández-Lobato. Bayesian Deep Learning via Subnetwork Inference. In International Conference on Machine Learning, 2021. Cited on page L
[16] D. Dua and C. Graff. UCI Machine Learning Repository, 2017. Cited on page 8
[17] M. Elahi, F. Ricci, and N. Rubens. A survey of active learning in collaborative filtering recommender systems. Computer Science Review, 2016. Cited on page |
[18] M. Filippone and R. Engler. Enabling scalable stochastic gradient-based inference for Gaussian
processes by employing the Unbiased LInear System SolvEr (ULISSE). In International
Conference on Machine Learning, 2015. Cited on page 2
[19] J. R. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. GPyTorch: Blackbox Matrix-matrix Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems, 2018. Cited on pages [
[20] M. N. Gibbs and D. J. C. MacKay. Efficient Implementation of Gaussian Processes for Interpolation. Technical report, University Cambridge, 1996. Cited on pages [
[21] R. Komez-Bombarelli, J. N. Wei, D. Duvenaud, J. M. Hernández-Lobato, B. Sánchez- T. D. Hirzel, R. P. Adams, and A. Aspuru- Guzik. Automatic Chemical Design Using a DrivnD Continuous Representation of
Lengeling, D. Sheberla, J.
Molecules. American Chemical Society Central Science, 2018. Cited on page [
11

---

[22] P. C. Hansen. Rank-deficiee and discrete ill-posed problems: numerical aspects of linear inversion. Society for Industrial and Applied Mathematics, 1998. Cited on page
[23] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian Processes for Big Data. In Uncertainty in Artificial Intelligence, 2013. Cited on p 8
[24] J. M. Hemmándzz-Lobtoo, M. W. Hoffman, and 7. Ghahramani, Predictive Entropy Search for Efficient Global Optimization of Black- -box Functions. In Advances in Neural Information
Processing Systems, 2014. Cited on page !
[25] J. M. Hemándezz-obatoo Requeima, E. O. Pyzer-Knapp, and A. Aspuru-Guzik. Parallel and Distributed Thompson Sampling for Large-scale Accelerated Exploration of Chemical Space. In International Conference on Machine Leurning, 2017. Cited on page
[26] R. A. Horn and C. R. Johnson, Matrix Analysis. Cambridge University Press, 2012. Cited on
page26 26
[27] B. Jin and Ž. Kereta, On the convergence of stochastic gradient descent for linear inverse problems in Banach spaces. SIAM Journal on Imaging Sciences, 2023. Cited on page 7
[28]
S. Lang. Real and Functional Analysis. Springer, 2012. Cited on page 27
[29] S. Mandt, M. D. Hoffman, and D. M. Blei. Stochastic Gradient Descent as Approximate Bayesian Inference. Journal of Machine Learning Research, 2017. Cited on page
[30] A. G. d. G. Matthews, J. Hensman, Turner, and Z. Ghahramani. On Sparse Variational Methods and the Kullback-Leibbe
Divergence between Stochastic Processes. In Artificial Intelligence and Statistics, Cadiz, Spain, 2016. Cited on page [8]
[3]] A. G. d. G. Matthews, J. Hron, R. E. Turner, and Z. Ghahramani. Sample-then-opimmize poste- rior sampling for Bayesian linear models. In NeurlPS Workshop Advances in Approximate Bayesian Inference, 2017. Cited on page
[32] T. Pinder and D. Dodd. GPJax: A Gaussian Process Framework in JAX. Journal of Open Source Software, 2022. Cited on page g
[33] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: the Art of Scientific Computing. Cambridge University Press, 2007. Cited on page 3
[34] A. Rahimi and B. Recht. Random Features for Large-scale Kernel Machines. In Advances in Neural Information Processing Systems, 2008. Cited on page 2
[35] C. E. Rasmussen and C. K. Williams. Gaussian Processes for Machine Learning. MIT Press,
2006. Cited on page 2
[36] N. Rubens, M. Elahi, M. Sugiyama, and D. Kaplan. Active Learning Recommender Systems. Recommender Systems Handbook, 2015. Cited on page
[37] S. Ruder. An overview of gradient descent optimization algorithms. 169i.01609.04777 2016. Cited on page [
[38] B. Schölkopf, R. Herbrich, and J. Smola. A Generalized Representer Theorem. In Compu- tational Learning Theory, 2001. Cited on p
[39] J. Snoek, H. Larochelle, and R. P. Adams. Practical Bayesian Optimization of Machine Learning Algorithms. In Advances in Neural Information Processing Systems, 2012. Cited on page[[]
[40] D. J. Sutherland and J. Schneider. On the error of random Fourier features. In Uncertainty in Artificial Intelligence, 2015, Cited on page
[41] A. Terenin, R. Burt, A. Artemev, S. Flaxman, M. van der Wilk, C. E. Rasmussen, and H. Ge. Numerically Stable Sparse Gaussian Processes via Minimum Separation using Cover Trees. Journal of Machine Learning Research, 2023. Cited on pages g
[42] Y. Tian, Y. Zhang, and H. Zhang. Recent Advances in Stochastic Gradient Descent in Deep Learning. Mathematics, 2023. Cited on P
[43] M. K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In
Artificial Intelligence and Statistics, 2009. Cited on pages D
18
[44] M. K. Titsias. Variational model selection for sparse Gaussian process regression. Technical
report, University of Manchester, 2009. Cited on page [8]
[45] V. Vapnik. The Nature of Statistical Learning. Springer, 1995. Cited on page 2
[46] K. A. Wang, G. Pleiss, J. R. Gardner, S. Tyree, K. Q. Weinberger, and A. G. Wilson. Exact
Gaussian Processes on a Million Data Points. In Advances in Neural Information Processing
Systems, 2019, Cited on pages g
[4
12

---

[47] V. Wild, M. Kanagawa, and D. Sejdinovic. Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes. arXiv:2106.0122,, 2023. Cited on pages S
18
[48] W. J. Wilkinson, P. E. Chang, M. R. Andersen, and A. Solin. State Space Expectation Prop- agation: Efficient Inference Schemes for Temporal Gaussian Processes. International Conference on Machine Learning, 2020. Cited on page
[49] G. Wilson and H. Nickisch. Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP). In International Conference on Machine Learning, 2015. Cited on page |I [50] J. T. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Efficiently Sampling Functions from Gaussian Process Posteriors. In International Conference on Machine
Learning. 2020. Cited on pages 2 9
[51] J. Wilson, V. Borovitskiy, A. Terenin, P. Mostowsky, and M. P. Deisenroth. Pathwise Conditioning of Gaussian Processes. Journal of Machine Learning Research, 2021. Cited on page
[52]
Zou, Wu, V. Braverman, Q. Gu, and S. M. Kakade. Benign Overfitting of Constant- Stepsize SGD for Linear Regression. In Conference on Learning Theory, 2021. Cited on pages]
3
13

---

A Experimental Setup
We employ the CG implementation of the JAX SciPy module and follow Wang al. [46) in using a pivoted Cholesky preconditioner of size 100. Our preconditioner implementation resembles the implementation of the TensorFlow Probability library. For small subset of datasets, we find the preconditioner to lead to slower convergence, and we report the results for conjugate gradients without preconditioning instead.
In addition to the experiment hyperparameters described in Section for all methods, we use 2000 random Fourier features draw each prior function used for computing posterior function samples via pathwise conditioning in
For SGD, each step, draw 100 new Fourier features to estimate the regularizer term. In all SGD experiments, we use a Nesterov momentum value 0.9 and Polyak averaging, following Antorán et al. 3]. For all regression experiments we use learning rate of 0.5 to estimate the mean function representer weights, and a learning rate of 0.1 to draw samples. For Thompson sampling, we use a learning rate of 0.3 for the mean and 0.0003 for the samples. In both settings, we perform gradient clipping using optax.
norm with max _orom__ommaooo set to 0.1.
Code available at:
A.1 Hyperparameter Selection for Regression
We prior mean function and the Matérn-3/2 kernel, and share hyperparameters across all methods, including baselines. For each dataset, we choose a homoscedastic Gaussian noise variance, kernel variance and separate length scale per input. For datasets with less than 50k observations, we tune these hyperparameters to maximize the exact GP marginal likelihood. The cubic cost of this procedure makes intractable a larger scale: instead, for datasets with more than 50k observations, we obtain hyperparameters using the following procedure:
1. From the training data, select a centroid data point uniformly at random.
2. Select subset of 10k data points with the smallest Euclidean distance the centroid.
3. Find hyperparameters by maximizing the exact GP marginal likelihood using this subset.
4. Using 10 different centroids, repeat the preceding steps and average the hyperparameters.
This approach avoids aliasing bias 515] due to data subsampling and tractable for large datasets.
A.2 Inducing Point Selection
For the inducing point SGD experiment shown in Figure 2, we choose inducing points as a subset of the training points. Due to the large number of datapoints in dataset, found k-means to converge very slowly. Instead, we develop an ad-hoc point elimination algorithm based on k- nearest-neighbors (KNN). We use the KNN implementation ANNOY which first run on the HIUSEELEECRICC dataset with num_ trees set to 50. We then iterate through our dataset, retrieving 100 nearest neighbors for each original point. Of these 100 nearest neighbors, only the ones closer to the original point than some length scale parameter, in terms of l2 distance, are selected. If the number of points selected larger than one, we eliminate the original point from our dataset and we also eliminate other points which are within the length scale neighborhood of the original and selected points simultaneously. The selected points are kept. We vary the number of points eliminated by our algorithm through the modification of the length scale parameter.
A.3 Maximizing Posterior Functions for Thompson Sampling
We maximize random acquisition function sampled from the GP posterior in three step process:
1. Evaluate the posterior function sample at a large number of nearby input locations. We find nearby locations using a combination of exploration and exploitation. For exploration, we sample locations uniformly at random from [0. 1]d, For exploitation, we first subsample the training data with probabilities proportional to the observed objective function values and then add Gaussian noise nearby N(0, o micurty), ). where neabby 1/2 and is the kernel
14

---

Objective function
Acquisition functions
Observed values
Ncarby locations
Top ncarby locations
-1
Acquisition maximizers
Optimization trace
-2
0.0
0.1
0.2
0.3
0.1
0.5
0.6
0.7
0.8
0.9
1.0
Figure 7: Illustration of single Thompson sampling acquisition maximization step on a 1D problem.
length scale. Throughout experiments, we find 10% of nearby locations using the uniform exploration strategy and 90% using the exploitation strategy.
2. Select the nearby locations which have the highest acquisition function values. To find the top nearby locations, we first try 50k nearby locations and then identify the location with the highest acquisition function value. We repeat this process 30 times, finding and evaluating total of 1.5m nearby locations and obtaining total of 30 top nearby locations.
3. Maximize the acquisition function with gradient- gradient-based optimization, using the top nearby locations as initialization. After optimization, the best location becomes Knew the maximizer at which the true objective function will be evaluated in the next acquisition step. Initializing the top nearby locations, we perform 100 steps of Adam on the sampled acquisition function, with learning rate of 0.001 to find the maximizer.
In every Thompson step, we perform this three step process in parallel for 1000 random acquisition functions sampled from the GP posterior, resulting in a total 1000 Znew, which will be added to the training data and evaluated at the objective function. Note that, although share the initial nearby locations between sampled acquisition functions, each acquisition function will, in general, produce distinct top nearby locations and maximizers. Figure 7 illustrates single Thompson step 1D) problem using acquisition functions, nearby locations and top nearby locations.
B Additional Experimental Results
Figure X provides optimization traces for SGD, CG and their low noise variants on all datasets with under 50k points. For these, we can compute the exact GP predictions, allowing us to evaluate error with respect to the exact GP in prediction space, representer weight space and in the RKHS. We observe the same trends reported in the main text. SGD decreases its prediction error monotonically, while CG does not. However, both take a similar amount of time to converge on smaller datasets and CG obtains lower error the end. While SGD makes negligible progress representer weight space, CG finds the optimal representer weights. CG's time to convergence is greatly increased in the low noise setting, but SGD practically unaffected.
Table 2 provides quantitative results for inducing point SGD on the HOUSEELECTRIC dataset. SGD's time to convergence is shown scale roughly linearly in the number of (inducing) points observed. However, for this dataset, keeping only 10% of observations and thus obtaining 10x faster
Table 2: Time to convergence and predictive performance for all methods under consideration, including inducing point SGD, on the HOUSEELECTRIC dataset.
Model
Inducing SGD
Standard SGD
CG 844 352
SVGP 499
M
218 782
431 489
1 1099 206
844 352
RMSE
0.08 = 0.00
0.08 0.00
0.08 + 0.01
0.09 0.00
0.11 0.11
Minutes
16.7 0.54
49.6 10.2 24.6 5.440
54.3
157
2.00
NLL
-1.10 0.05
-1.11 + 0.04
-1.13 0.04
-1.09 0.04
2.07
0.10
15

---

RMSE to Test Data
RMSE to Exact GP
Euclidean Error
RKHS Error
1.2
45
20
1
pol
0.8
30
0
0
600
2.5
0.4
0.0
0
400
1.0
3.6
0.9
3.5
0.6
elevators
2.4
0.6
0.0
0.0
1.2
0.3
160
5
0.0
0.0
30
0
1.2
1.8
100
1.5
bike
0.8
1.2
0
0.0
0.4
0.6
160
1.5
0.0
0.0
0.0
2.2
1.8
6
1
protein
1.6
1.2
0
0
1.0
0.6
25k
1.8
0.1
0.0
23k
0.6
10
10
10'
102 10
10
10
10'
10"
10a
10
10
10
10a
10a
10
10
10
10 10a
Time in sec
Time in sec
Time in sec
Time in sec
SGD (optimal noise)
SGD (low noise)
CG (optimal noise)
CG (low noise)
Figure 8: Convergence of GP posterior mean with SGD and CG as a function of time (on an A100
GPU) the POL (N 15k), ELEVATORS (N 16k), BIKE (N As 17k) and PROTEIN (N 46k) datasets while setting the noise scale to (i) maximize exact GP marginal likelihood and (ii) to 10
labeled low noise. We plot, in left-to-right order, test RMSE, RMSE the exact GP mean at the test inputs, representer weight error a| v" l1 and RKHS error || - PSGD || H2.
Table 3: UCI benchmark with subset of data on the smaller datasets, using a randomly-sampled subset with size that varies as percentage of the full data.
Dataset N
POL.
BLEVATORS
BIKE 17379
PROTEIN
KEGGDIR
15000
16599
45730
48827
5%
0.20
0.18 0.00
0.71 0.00 0.66 0.10
RMSE
10%
0.16 0.00
0.11
5%
-0.43 = 0.02
0.01
0.09
0.90 -0.7
NLL
10%
-0.66 = 0.02
0.51 0.01
0.07
0.91 0.01 0.08
convergence leaves performance unaffected. This suggests the dataset can be summarized well by small number of points. Indeed, SVGP obtains almost as strong performance as SGD in terms of RMSE with only 4096 inducing points. SVGP's NLL is weaker however, which is consistent with known issues of uncertainty overestimation when using a too small amount of inducing points. On the other hand, the large and potentially redundant nature of this dataset makes the corresponding optimization problem ill-conditioned, hurting CG's performance. Tables 3 and 4 provide additional baseline results for the subset data approximation, which behaves similarly SVGP.
Finally, Figure D shows all methods' Thompson sampling performance as a function of both compute time and acquisition steps for each individual true function length scale considered in our experiments. Differences among methods are more pronounced in the small compute budget setting. Here, SVGP performs well-on par with SGD-in in the large length scale setting, where many observations can likely be summarized with 1024 inducing points. CG suffers from slow convergence due to ill- conditioning here. On the other hand, CG performs length scale setting, while SVGP suffers. In the large compute setting, all methods perform similarly per acquisition step for all length scales except the small one, where SVGP suffers.
par with SGD in the
small
16

---

Table 4: UCI benchmark with subset of data the larger datasets, using randomly-sampled subsets of fixed size which are chosen to reflect typical GPU memory limitations. Numbers which are better than or competitive with respect to Table are printed in boldface.
Dataset N
3DROAD
SONG
BUZZ
HOUSEELEC
434874
515345
583250
2049280
25k
0.23 + 0.00
0.82 0.00
0.33
0.06 0.01
RMSE
50k
0.06
0.81 0.32 0.00
0.00
0.06 0.00
25k
-0.42
= 0.01
1.22
0.30 0.06
-0.53 + 0.33
NLL
50k
-0.78
0.01
1.20 + 0.00
0.26 + 0.06
-0.49 + 0.39
Small Compute Budget
Large Compute Budget
7.0
Max value
1= 0.1
4.5
4.0
Max value
6.0
1= 0.2
5.0
4.0
6.5
Max value
6.0
= 0.3
5.5
4.5
6.0
Max value
= 0.4
5.0
4.0
3.5
5.5
Max value
I= 0.5
5.0
4.5
3.5
50k
60k
70k
80kc
0h
1h
2h
3h
50k
60k
70k
80k
0h
5h
10h
15h
Acquisitions
Time elapsed
Acquisitions
Time elapsed
Random Search
SGD
SVGP
Figure 9: Maximum function values (mean and std. err.) obtained by Thompson sampling with our approximate inference methods function of acquisition steps and of compute time on an A100 GPU. The latter includes time spent drawing function samples and finding their maxima. All methods share starting dataset of 50k points and we take Thompson steps, acquiring 1000 points in each. Different length scales / exhibit different behaviors: CG performs better in settings with smaller length scales and better conditioning, SVGP tends to perform better in settings with larger setting and increased smoothness, SGD performs well in both settings.
17

---

C Inducing Points
Let z XMM be set of M inducing point locations. The variational inducing point framework
of Titsias 43 44] substitutes our observed targets with the inducing targets u E IR' Ad Our Gaussian
process, conditional on a set of inducing locations and targets, given by the mean and covariance functions
p))(() K (-)x K
k kg(((,-) k K2 x1)-
(17)
Titsias [43 provides a closed-form expression for the variational distribution N(u K(*))
over inducing targets which minimizes the associated Kullback-Leiberr divergence, namely
" K= (K KI E E- K2 E-
(18)
K (2) = (Kez + K Σ- "KZz) K
(19)
In turn, Matthews et al. (30] shows that these expressions lift respective infnite-ddimensional Kullback-Leible minimization problem between the variational and true posterior Gaussian pro- cesses. Using this, we obtain a variational posterior with mean and covariance functions
F (2) ()) K ())((=== + E 'K "Kz)) E
(20)
k (-) fly ((.)) ()= + E= - -K=) JK 20))
(21)
These expressions will be our starting point. We now write this Gaussian process posterior in pathwise manner and derive stochastic estimators for the corresponding representer weights.
C.1 Pathwise Representation of the Kullback-Leibler-pptinn Variational Distribution
The pathwise expression of Section 3.3 can be derived by restricting the domain of the respective optimization problem which describes the exact posterior samples. This is done by replacing the reproducing kernel Hilbert space with subset consisting of the span of set of canonical basis functions centered at the inducing points, as described by Wild et al. (47). Theorem for the posterior mean. Rather than deriving the respective result from scratch, we will verify its correctness by means of computing the relevant means and covariances. Define
10n-f(a1 ( - E) (f(=) v-) = f() + K K K K Σ) -
(22)
N(O,E) GP(O f(=)()) K
We now calculate the moments of this Gaussian process' marginal distributions and show them to match those of the KL-optimal variational Gaussian process given 20
Write
[((((z) (.) K() K K--Kaa z (K= K
z Σ) -
(23)
K(j K Σ-1 (Kez =z
Kezk
K
+ +I)
(24)
-()z (ax E- -
Kus + +x,z)
''K
(25)
(26)
and
Cov((f(*) v- - ( (1)
(27)
= [(()(*) v)()) - 1 (0) ((x) y)(,) (1))
(28)
'(())) (() Ko>=K_IF K- z K mz K K- Kze xx + Σ) K K Kz())
(29)
= -K(' + KorKK, (-I+I K E-1 (Kus K_ K E- +I) -K K-] Kz())
(30)
= -K(i + K()=K,] (-I+ (K E Krzk k= zz +I)-1) Kz(1)
(31)
=K('") + K()= (-K + (K E Kez + Kzz)-1) Kz(1)
(32)
= -A(e)(-.")
(33)
which recovers 21. as claimed,
18

---

C.2
Derivation of inducing point optimization objectives
We now derive the optimization objectives we use for inducing points.
The MAP objective (22) To find the optimization objective for the inducing point mean function's 's
representer weights, we apply 23 and begin from
f1y (-) (() (Kz E- KKn
Kzz) 22 Σ E (()
(34)
v² as the expression for the optimizer of ridge-regularized linear regression problem see Bishop (7), Chapter 3-with parameters
Now, we recognize 2 "=z)
features Kgz. Gaussian noise of covariance Σ, and regularizer Kzz. This means that its respective optimization problem
arg min vERM
K ||
(35)
The sampling objective 13 For the uncertainty reduction term's representer weights with induc- giving
ing points, also leverage
(())= Ku Kz (Ker K, + E) +X)))(((((())) +
(36)
= ()))((Ezn E Kzz -'Kz k=z)-KKa E E (((()((()))+
(37)
= =K
(38)
This is again therefore the solution to regge-regularized quadratic problem, but now the targets are given by the random variable (f(=) + t). We thus write the objective
arg min RERa
( (z) + K 1 + ||ae| K...
(39)
C.3 On the Error in the Nyström Approximation K2
Kaz
Figure I0 compares the KL-optimal inducing point posterior GP with that obtained when applying
the approximation presented Section 3.3 That is, taking the prior function samples which we fit with the representer weighed canonical basis functions to be f(x) with f GP k) instead of ((= = K= This amounts to approximating the Nyström-type matrix K K=z K with its exact-posterior counterpart K The difference between them exactly equal to posterior covariance: by posterior contraction, both of these matrices become very similar if there is an inducing point placed sufficiently close to every data point. In practice, this tends to occur when an inducing point is placed within roughly helf-leeggh--cale of every observation. This is effectively what needed for inducing point methods to provide good approximation of the exact GP. This is reflected in Figure TO where see that our approximate inducing point posterior differs from the exact inducing point posterior only in situations where the latter fails to good approximation to the exact GP in the first place. This manifests the approximate method providing larger error bars. When the number of inducing points increases, both methods become indistinguishable from each other and the exact GP. Fortunately, the linear cost of SGD in the number of inducing points allows us to use a very large number of these in practice.
D Low-variance Sampling Objectives
Here, we show that the modified sampling objective (10) matches the standard sample-then-optimize objective (9) up constant. We then compare both objectives' minibatch gradient variance.
D.1 Standard and Modified Chapple-then-Optimiie Objectives: Equivalent Gradients
Let LLT = Σ be the Cholesky factorization of the noise covariance, let f(x) ~ Kew), and let N I). Consider the objectives hand, namely
f + Ls Ke 3 1 ||< I
(40)
19

---

10 Inducing Points
15 Inducing Points
20 Inducing Points
Infill asymptotics
-3
1000 Inducing Points
2000 Inducing Points
3000 Inducing Points
M
Large-domain
asymptotics
-7
-3.5
0
3.5
7
-7
-3.5
3.5
7 -7 -3.5
3.5
7
Exact GP
Approx. Inducing GP
Inducing GP
Figure 10: Comparison of exact and approximate inducing point posteriors for a GP with squared exponential kernel and 10k data points generated using the true regression function sin(2x) cos under two different data-generation schemes: infill asymptotics, which considers T; ~ N(0, 1), and large-domain asymptotics, which considers I; on evenly spaced grid with fixed spacing. We see that the approximation needed to apply inducing points is only inaccurate in situations where the inducing point posterior itself has significant error, which generally manifests itself bars that are larger than those of the exact posterior.
and
[[f(z) - 1 -
(41)
K,n
We show these are equal up to a constant by showing they have the same derivatives. Taking derivative with respect to a, we have
0 [1(x)) + Ls K 0 ||ax|
(42)
da
=-2K E (f(z) + Ls Kx + 2K 2Kaza
(43)
= 2K [ ( - E -E'K + +L ox)
(44)
and
0
da 1f(x)) K
s||
(45)
--2Ka W (f(z) a + 2K (a L-T's
(46)
= =-2K E f((((x) E - L - aa)
(47)
respectively. These expressions match, giving the claim.
D.2 Minibatch Gradient Variance
Following Antorán et al. (3). we consider single sample minibatch gradient estimators for the data fit
term in ( that D For each of the two sampling objectives under consideration, the gradient estimators are
Hoad -N E,.U (I....) Kn... + - Km)
(48)
ØT 0m - N Eiott( N) K (f(xz) a)
(49)
20

---

for 99 and (TO), respectively. Since both objectives use the same Fourier feature estimator for the gradients of the regularizer, we omit these from our analysis. The variances for single sample gradient estimators of both expressions are given by
Cov(@Lola) = N Cov(K (f(x) + Ke,+a))
(50)
Cov(OLnew) N Cov(K (f(x) K a)),
(51)
respectively, Expanding both expressions and subtracting them arrive
Cov(OLuua Cov( N Cov(K 2N Ke K
(52)
Whether proposed sampling objective presents lower variance rests the positivity of the above matrix. To analyze this, we must give concrete values to a. Following Antorán et al. we consider
two settings, initialization and convergence.
(I) Initialization. Here, a and therefore
(o Cov((LLew) N Cov N (K EK xxx)
(53)
a positive definite matrix. Thus, our proposed objective has lower variance.
(II) Convergence. Here, a = (Kuu + Σ) - (f(x)) + E), so we have
Cov(OLoia) - Cov = Cov (ov(KæzE) - 22V Cov Kar K
(54)
= NK EK 2 2NK mra (Kar + Σ) "EEK
ΣK
(55)
Letting E o²I and UAUT = Ke be the eigendecomposition of the kernel matrix, have
NK EK ax 2NK2 (Kaz nm + -EK
(56)
NUA 2A(A + +
(57)
Thus, convergence, whether our proposed objective reduces variance rests whether not the matrix I- - 2A(A + +o²1) is positive definite.
Following Section 3.4 in practice, the representer weights only converge along the top few principal directions. For these, the corresponding eigenvalues will likely much larger than the noise and thus 2x 0. However, in the vast majority of directions, the opposite true: this yields lower variance the minibatch estimator of our proposed objective, which is shown in Figure 2
E
The Implicit Bias of SGD GPs
This section provides our main thcoretical results.
E.1
Convergence of Stochastic Gradient Descent
Before studying what happens in Gaussian processes, we first prove a fairly standard result on the convergence of SGD for a quadratic objective appearing in kernel ridge regression, which represents the posterior mean. For simplicity, we do not analyze Nesterov momentum, nor aspects such gradient clipping. We will take the Gaussian observation noise covariance be a constant multiple of the identity in this subsection. Recall that random variable 2 is called G-sub- Suub-Gaussian if E(exp(X(e (x()=)()))) S cxp(3 GX?) holds for all Then, random vector is called its dot product with any unit vector is E-sub-Gaussiae to having tails that are no heavier than those of a Gaussian random vector, formulated appropriately. Let UA U K be the eigenvalue decomposition of the kernel matrix. The eigenvectors and eigenvalues X; are given in descending order: X₁ 2
One can show this condition is equivalent
AN 0.
Lemma 2. Let ò > and Σ o²1 o² > 0. Let & be a sufficiently small learning rate of ?²
0 < イイ < A1(A Let a² E R be the solution of the respective linear system, and let an be the Polyak-aves averaged SGD iterate after t steps, starting from an initial condition of ax = Assume that the stochastic estimate of the gradient is unbiased and G-sub Sauu-Gasssian Then, with probability 1 8, we have for any I C (1,..,N N} and all E I that
|1 + G n
[
|u (a ax)l < A
log
(58)
21

---

Proof. Consider the objective
1 N [(a)= 20²
E (v= K = a a'KKza K
(59)
2
and respective gradient
0L a
1 a²
1 o² (K (Kza + +ooI)a -Kaay)
(K (k a- K + K
(60)
da
Let us first look at non-stochastic gradient optimization of I,, without Polyak averaging. The iterate O4 is given by
0L = da (at-1) (I - K (Kaz + o²1) a-
2
Kazy
(61)
0²
0²
Writing M = (I- t1 Kzz (Ee + oo))),, and recalling that an = 0, we thus have that
1-1
a== 0² M'Keoy= 7) (T (I-M) M K = a (Kau o'I)-IM'' (62)
"7
Σ 1=0
where, for our choice of learning rate, (i) by direct calculation using simultaneous diagonalizability of M and Ke the matrix M has eigenvalues whose absolute values are strictly less than 1, (ii) thus, the geometric series converges, (iii) and I M is positive definite. Examining the error in direction u, we use simultaneous diagonalizability see that
(1 nd (A3 + +
luf(a" (a" a:)l = |ui (Kev + 1-1)--M''
1u y
(63)
d 0²
nd (1: + +
||
(64)
A o2
which applies to ordinary gradient descent without stochastic gradients or Polyak averaging. Next, consider stochastic gradient optimization. We will first consider the case where the gradient is independently perturbed by C, ~ N(0, I) for each step t >> 0, and then relax this to sub-Gaussian noise in the sequel. Specifically, consider
a = -0=-1-1- -
(al-1) + (o) 0z,
AL do
Σ M
1-1
(65)
1=0
with 01j = 0kg = 0. For these iterates, consider the respective Polyak-averaged iterates denoted by
σ; -1 a's. We have
|u/(a' (a' - -a.)l = Σ uf (a (a- -aú)) < Σ u (a" -aj))
Σ u (a, - az)
(66)
B,, by expanding the Polyak averages and applying the triangle inequality, and where we have introduced the notation Ar, and B Using (63), the triangle inequality and another geometric series, we can bound the first sum as
ALe
1 (A + 02)
(1- n 1-
0² 1111
Ar.l
(67)
A); + +
For the second sum, note that we can re-index the order of summation to to count each C, only once, rather than once for each Polyak average. This gives
t "
Big = Σ
u
M
C j-4
1
Σ
u
M )s
(68)
t 5=1
f=
9=0
7
ΣΣ
1-
n71 g²
(d: (2: "uu
(69)
22

---

where u.C, ~ N(0, 1) are IID. The variance of B; is bounded using another geometric series
Var(B,,)) Var
02
(M++
ut
(70)
112
Var
Σ
1
(^; + 02)
u u.:
(71)
t2
o²
9=0
1
n
(1:+
Var un (5)
(72)
12
9=0
1
0²
(73)
X²(); + 02)
Thus, from standard tail inequalities for Gaussian random variables, for any fixed 8' we have
20 log 8 2 21-8'
IP B S V x + 02)21
(74)
We then take 8' 8/] and apply the union bound for all indices in I. Finally, by comparing moment generating functions, we can relax the Gaussian assumption to Suuu-aaasii random variables. To complete the claim, combine the bounds for the two sums by writing
a² 1
20
Z log
|ui (a" Tx.) S |A |1rrdl
GV
(75)
mi: + 02)2t
2E(), + 02)2
0² Ai(di + 02)
1|
+ +GV =//2 log
[
(76)
n(A: + 02
|1ull2 A no²t
2
+ G log
(77)
which gives the claim.
E.2 Convergence of SGD in RKHS Subspaces
The idea is to use the preceding result to show that SGD converges fast with respect to certain seminorm. Let & be the kernel and let H& be its associated reproducing kernel Hilbert space. We now define the spectral basis functions from Section 3.4
Definition Let K UAUT be the respective eigendecomposition. Define the SPECTRAL BASIS FUNCTIONS
N Σ
Un k( VA
2 (0()))
(78)
These are precisely the basis functions which arise in kernel principal component analysis. We also introduce two subspaces of the RKHS: the span of the representer weights, and the span of a subset of spectral basis functions. These are defined as follows.
Definition 4. Define the REPRESENTER WEIGHT SPACE
R span :) :i= 1,.,,NJ NJ.. C
(79)
equipped with the subspace inner product.
Definition 5. Let I S N} be an arbitrary set of indices. Define the INTERPOLATION
SUBSPACE
R spant u(() : I R
(80)
23

---

We can use these subspaces to define a seminorm on the RKHS, which we will use to measure convergence rates of SGD. Definition 6. Define the INTERPOLATION SEMINORM
If R IIproj fl
(81)
Here, projs denotes orthogonal projection of function f onto the subspace S of Hilbert space: to ease notation, in cases where we project onto the span of single vector, we write the vector in place of S. The first order of business functions and representer weight space.
understand the relationship between the spectral basis
Lemma 7. The functions u for I,,,, form an orthonormal basis of Ra,
Proof. Let i / j. Using the reproducing property, definition of an eigenvector, the assumption that u an eigenvector of Kaz, and orthonormality of eigenvectors, we have
Ui K K
A VAA
( Hr.
u u, -
(82)
Viyyg
The fact that (9 form basis follows from the fact that they are linearly independent, that there N them in total, and that by definition this equals the dimension of R Finally, we calculate the norm:
2c U;KKa u.
1 1
A =1 1.
(83)
Next, we compute a change-of-basis formula between the canonical basis functions and spectral basis functions.
Lemma 8. For some vector 0 E RN, let
N
0(.) = Σ 0ch(ai,)) Σ ()
(84)
Then we have
w = A 1/2
|101 2
= w w.
(85)
H2
Proof. By expanding ( have
N
N
N
Us
8(-) Σ
w.
Wj
k(x).)) Vxg
(86)
j=1
i=1
j=1
which since k((2,)) is a basis means 0,
Σ
N
Wy
Ug V
or in matrix-vector notation 0
UA w1/ The final claim about norms follows from orthonormality of ( by writing
N
N
||01 (0,0)
wwwy(ul
(1) H2 = Σ w = w w.
(87)
This allows us to get an explicit expression for the previously-introduced seminorm.
Lemma 9. For O E Rk.m letting Az.
:., Ax ANIINeT), we have
10 R = - U
(88)
24

---

Proof. Decompose f in the above orthonormal basis, giving
0(-) = Σ wy (.))
(89)
By orthonormality, the definition of R and properties of projections, we have
(proj 0)(-) Σ
wyu(i))
(90)
jEI
Therefore, again using orthonormality, we obtain
10
Σn UAAU
(91)
JEI
which gives the claim.
With this at hand, we claim that our variant SGD converges quickly with respect to this seminorm.
Proposition 10. Under the assumptions of Lemma 2 for any I we have
|3112 1o²t
2 log
4 5
|f MSGD R
G
(92)
A
Proof. Using Lemma 9 and linearity, we have
1s MSGD
Σ co 2z (())(())))))
= - -' v*)
(93)
uf (v w')//>0=
|1w/l12 no²t
GV.
/
log
[
A
(94)
V
A
8
NET
|1 no²t
2
I
+ GV.
log
(95)
t
A
where the final inequality comes from substituting in the result of Lemma 2, yielding the claim.
Our main claim follows directly from the established framework.
Proposition 1. Let 5 > O. Let Σ = I f o² > 0. Let μSGD be the predictive mean obtained by Polyak-averaged SGD after steps, starting from an initial set of representer weights equal to zero, and using sufficiently small learning rate of < < 20 Assume the stochastic estimate of the gradient is SassubG-uussin Then, with probability 1 - ô, we have for = 1, N that
1
l2
2 t
ppoojat Pfw - projuco PSGDD)
+
G
log
(16)
VA
2t
5
Proof. Choose to be a singleton, and apply Proposition where by replacing [Z with N in the final bound the claim can be made to hold with probability 1-0 for all i = I,.., N simultaneously. O
This concludes the first part of the argument for why SGD a good idea: it converges fast with respect to this seminorm. In particular, taking I be the full index set, position-dependent pointwise convergence bound follows.
Corollary 11. Under the assumptions of Lemma
we have
||æ/|12
I
N
(usCD(-) - ())(())
G
log
w ( (-).
(96)
1
25

---

Proof. Using Lemma 8 write
N -v)). where w A2/2 (v""". Applying Lemma 2 with I {u,.,, we conclude
N
PsGD()) ())v()) Σ
(1) 1z ())k(())))
Σ
w w )u (-)
(97)
N wi (+)
(MGGD()) 11/1
wy
v (v
(98)
N
|1u1|2
[
A
+ G
log
(99)
==1
-llella + GV
[
N
log
(0)
(100)
8
Vi:
The claim follows.
Examining the expression, the approximation error will be high at locations E X where (x) corresponding to tail cigenfunctions is large. We proceed to analyze when this occurs.
E.3
Courant-Fischer in the Reproducing Kernel Hilbert Space
Since the preceding seminorm induced by projection within the RKHS, the step is understand what kind of functions the corresponding subspace contains. To get qualitative view, the basic idea is to lift the Courant-Fischer characterization of eigenvalues and eigenvectors to the RKHS.
We will need the following variant the Min-Max Theorem, which is useful for describing multiple eigenvectors at once-as solutions to sequence of Raleigh-quotient-tye optimization problems, where the next optimization problem performed on a subspace which ensures its solution is orthogonal to all previous solutions. Mirroring the rest of this work, consider eigenvalues in decreasing order, namely λ1 2 2 AN= 2 O.
Result 12. Let A be a positive semi-defiitte matrix. Then
n Au
A
max (OO)
u² u
(101)
U E.jjsi
where the eigenvectors u, are the respective maximizers.
Proof. Horn and Johnson 226.. Theorem 4.2.2, modified uppropriately handle eigenvalues in decreasing order, and where we choose i,, ... d = I,.., - +
We will prove the main claim below in two stages: first, on representer weight space Rk..r then on the full RKHS H₃. We begin with the former.
Lemma 13. We have
a ()) arg max
u(xx)) ||ull 1, (u, 6) H2 0.Vj <i
(102)
Proof. Define wy -A21 1/, where we recall that u, are the respective cigenvectors of the kernel matrix. From the Result form of the CouranttFFischer Min-Max Theorem, the reproducing property, and the explicit form of the RKHS inner product on R in the basis defined by the canonical basis functions,
conclude that
u ''KKnu
W
K mm
w
A MER (()
(103)
n" u
WER (())
wt K
" ViVjeii
wn K =.Vis
26

---

Ein wak(wrs xt)
max WERNMO (0)
max WER N ((0)
(104)
wI k
Σ wyk(Is;) :)
w" n Vi
Knn =.,jii
Hx.
llu(x))
max = 100 (u.a Ha 0,95s
max uERe,
Σ u(2x)²
(105)
|luli H
|| My (a,u")),, =0,ijss
which gives the claim.
Next, we use a Representer- Siiieenter-TThromm-typeee orthogonality argument to extend the optimization problem to all of Hg. The argument will rely heavily the Projection Theorem for Hilbert spaces: for overview of the setting, see Lang (28). Chapter
Proposition 14. We have
u (+) = arg max
Σ u(x;)² || =1, (u, u0) =.Wij<i
(106)
Proof. Let (.()) be the orthogonal complement of the finiteddimensional subspace span (u0), , 1, 1..., i - 1} in Hg and note by feasibility that ( E H Consider the decomposition of (
into its projection onto an orthogonal complement with respect to the finite-dimensional subspace k n RR... namely
(.) (.)(+)+ (-)
(107)
where this notation is defined as ull proj (( Ru, alt) and u- ul() -ll. Observe that
w- (x;) = (w+ k(x;, )) hr = (proj w-, ,,((,, = (wt,, projproj H k(xx,,) H) = (108)
where the first equality follows from the reproducing property, the second equality follows because w- ( the third equality follows because the projection is orthogonal and therefore self-adjoint, and last equality follows since u- by definition lives in the orthogonal complement of H n Re within H!)) with the subspace inner product. Moreover, by the Projection Theorem and the fact that
H inherits its norm from Hę, note that || | + ||u-112
(109)
We now argue by contradiction: suppose that u' # 0. This means that ||u || above, we have
0. From the
N
(zx) )2
u (x,)²
| (1(
(110)
1=1
1=1
Define the function
[lu(() |1e (.)
ul (-)
(111)
|uul Hx
and note that ||,'||| 1, and (u', (()) = 0 for all j < since ulll E ( which makes w' feasible solution to the optimization problem considered. At the same time, satisfies
[ 12 Σ ulll (xx
N
Σ 2
Σ 2x.2
(
(112)
which shows that (1) contrary to its definition, is not actually optimal-contraaictinnn contradiction. We conclude that u-t which means that (i) E H!!0) Rk, and the claim follows from Lemma
This means that the top eigenvalues can be understood as the largest possible squared sums of canonical basis functions evaluated at the data, subject to RKHS-norm constraints. situations where k(r,,) :) is continuous, non-negative, and decays as one moves away from the data, is intuitively clear that such sums will maximized by placing weight on canonical basis functions which cover much of the data possible. This mirrors the RKHS view of kernel principal component analysis.
27


Reviews:
[{"review_id": "xXJ4cSkmBO", "reviewer": "Reviewer_Hdii", "summary": "This paper introduces a method for fast approximating a Gaussian process posterior when the data size is large. Exact computation complexity would be cubic in the data size, while this method is linear. It originates from the idea of pairwise conditioning of Gaussian process, where the law of Gaussian process posterior is expressed in terms of the law of Gaussian process prior. Decomposition into eigenfunctions gives a way of approximating the Gaussian process posterior, so the posterior inference transforms into appropriately choosing the coefficients of decomposition. Objectives are formed quadratic in the coefficients, which can be optimized with SGD for the sake of lower computational complexity compared to conjugate gradient methods. Ideas of inducing points are also discussed for reducing the data size needed. Intuitive discussion of errors in different regions are companied by figure illustrations as well as some theoretical results. Adequate numerical experiments are presented to support the method.", "questions": "To control the approximation error below a constant threshold, how does the number of basis L (number of components) needed scales with the data size and dimension? Either theoretical or numerical result could be interesting.", "limitations": "See Weaknesses and Questions.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1. The paper has clear descriptions and is well written.\n2. The ideas are mostly original and combine the advantages of multiple methods.\n3. Many figure illustrations are present, making the ideas easily understood.\n4. Supportive numerical experiments are conducted.\n5. The problem of reducing the computational complexity of Gaussian process posteriors is itself very important and meaningful.", "weaknesses": "It would be nice to discuss why the Fourier basis are used for eigen decomposition and how it performs compared to other basis such as wavelets.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nThank you very much for your review! We are delighted that you found our descriptions “clear” and our paper “well-written” - thank you for these comments!\n\n> “It would be nice to discuss why the Fourier basis are used for eigen decomposition and how it performs compared to other basis such as wavelets.”* and asks the question *“To control the approximation error below a constant threshold, how does the number of basis L (number of components) needed scales with the data size and dimension? Either theoretical or numerical result could be interesting.”\n\nThank you very much for these two closely-related questions! Due to the use of the terms “Fourier basis” and “spectral basis” in a non-synonymous way in our work, these questions can be interpreted in two ways: either (1) referring to the number of Fourier features, or to (2) the number of spectral basis functions along which we examine the convergence of SGD. Both questions are interesting, including potentially to other referees, so we will answer each of them.\n\n1. **Fourier features.**  \nIn our work, Fourier features are used to (a) approximate the prior samples needed for pathwise conditioning, and (b) to approximate the regularization term $ || \\alpha ||_K$  which appears in the quadratic objective used by SGD. Neither requires Fourier features explicitly: any finite basis function approximation of the prior kernel will work. None of our techniques are limited to Fourier features - we use them because they are convenient and work well for stationary kernels. For other kernels, including potentially non-stationary kernels, other bases such as wavelets or random hashes (Tripp et al. 2023, \"Tanimoto Random Features for Scalable Molecular Machine Learning\") could be used instead of Fourier features. This can be especially interesting for non-stationary kernels and kernels on boundary-constrained domains.\n\n* **On random (Fourier) feature approximation error**: We only use random features to approximate quantities that do not dependent on the targets: (a) prior function samples and (b) norms in the metric induced by the kernel matrix $\\|\\alpha\\|_K$. In (a), the approximation error is also independent of the number of observations. We found this also to be the case empirically in (b).\n\n* Indeed, we use the same number of random features for across all experiments; 2000 for prior sampling and 100 for kernel matrix norms.   \n\n* Crucially, we do *not* approximate any conditional distributions (i.e. matrix inverses) using random features; we use SGD for this. As the reviewer suggests, approximating conditionals with random features requires a number of features increasing in the number of observations and is prone to variance starvation. See Sutherland and Schneider (2015), \"On the Error of Random Fourier Features\" and Wilson et al. (2021), \"Pathwise Conditioning of Gaussian Processes\" for details including explicit error analysis.  \n\n\n\n\n2. **Spectral basis functions.**  \nOne can ask how many spectral basis functions one needs to look at in the convergence bound to ensure good performance. Recall that these are defined as $u^{(i)}(\\cdot) = \\sum_{j=1}^N U_{ji} / \\sqrt{\\lambda_i} k(x_j, \\cdot)$, where $U$ is the U-matrix in the eigendecomposition of $K_{xx}$ and $\\lambda_i$ are the respective eigenvalues. We chose the term “spectral” for these basis functions since they arise from the kernel matrix’s eigendecomposition. They define a “kernel-specific” basis that in some intuitive sense resembles the Fourier basis, but is different because it is data-dependent. As revealed by Proposition 1, the question of how many such basis functions are needed to ensure good performance is central to the performance of SGD in GPs, because it determines how many iterations are needed to ensure convergence. We do not address this question from a theoretical standpoint, because it is very technically difficult - please see our response to Referee Hc1f part (3) for more on this. "}, {"review_id": "N5XUPdw7hG", "reviewer": "Reviewer_Hc1f", "summary": "This work proposes SGD GP, a method based on stochastic gradient descent to efficiently compute the GP posterior samples given fixed hyperparameters. The method relies on the pathwise conditioning GP posterior formulation and random Fourier features (RFF) approximation. The key idea is to express the GP posterior quantities as solutions to quadratic optimization problems whose objective is a sum over data points and hence SGD can be applied. \n\nThe paper shows that SGD GP produce accurate predictions. However, SGD GP can converge slowly or converge to sub-optimum. However, non-convergence behaviors of SGD only occur in region closed to the data boundary. SGD GP performs comparably to SVGP and conjugate gradients (CG) in most settings, and can outperform them in large-scale systems or ill-conditioned problems.", "questions": "- Related to the first point above, can the method be adapted for broader settings where the hyperparameters need to be learned, or the likelihood is not Gaussian? If not, what is the limiting factor there?\n\n- Regarding the experiment evaluation, sec 4.1, is the predictive RMSE and NLL good metrics for evaluation? If I understand it correctly, SGD and CG are based on the same set of learned hyperpararmeters (and SVGP is based on a separately learned variational model). The core goal of the comparison is whether SGD recovers more faithful inference approximation compared to CG as a baseline. So I thought the valuable metric should be a distance against the \"exact\" inference result (e.g. CG with max. number of iterations and with high numeric precision). So I am not sure how to interpret Table 1. \n\n- In Fig 3 and Fig 5, do you have a sense why CG errors first go up and then go down? In Fig 5 last panel (houseelectric) would you expect CG error to match SGD performance eventually if running for more time? In both figures, I think it would be helpful to provide an exact baseline (e.g. CG ran to reach tolerance like1e-3). \n\nMinor comments. \n\n- A few notations and figures are not very clear. E.g. Fig 4. top left panel, indicate the size of error bands (blue shaded area) and dotted black line; top right corner, indicate the black dots (observations). In Proposition 1, define G-sub-gaussian.", "limitations": "The limitation discussion is missing. I don't think the paper would have negative societal impact. My main concern on technical limitaiton is weakness point 1.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- The paper is overall well-written and easy to follow. \n\n- The proposed method is novel and sound. It is shown to provide better predictive performance given the same inference time compared to SVGP and CG. What I found is most compelling is that SGD-GP seems to be a stronger alternative in large-scale or ill-conditioned systems.\n\n- I also find the spectral analysis of SGD convergence of three different regions to be very insightful.", "weaknesses": "- The major weakness is that the setting the paper considers is quite limited, which is the posterior inference given fixed learned hyperparameter. I would appreciate if the authors can elaborate on the significance of the setting and why the proposed methodology is in particular important. For example, how often would the ill-conditioned systems arise in practice and how common is the case that hyperparameters are known in advance. It also seems like the method is applied to a isotropic Gaussian likelihood (see the question section below). \n\n- The fact that SGD converges slowly in extrapolation region is a bit concerning. Especially in applications like Bayesian optimization, this is the region of high interests for exploration. Or in settings where there are distribution shifts, the under-calibrated uncertainty in this region can be an issue. In general, I found the \"benign non-convergence\" argument in the extrapolation region not convincing. Would appreciate if the authors can elaborate on this issue, and if if there is any potential way to alleviate the non-convergence property? \n\n- The SGD convergence analysis is insightful. But I would appreciate a more formal mathematical characterization of three regions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We wanted to start by thanking you for your time in reading our work and providing very helpful comments! We are thrilled you found our work “well-written and easy to follow” and agreed that our results are interesting because they provide a “stronger alternative in large-scale or ill-conditioned systems”, which as we will argue below include most systems of sufficient scale. \n\n----\n\nBelow we address the weaknesses and questions:\n\n1.  **Significance of setting**\n\n\n* *Significance*: You are correct that our proposed approach can not be used to tune hyperparameters and requires a Gaussian (*but not necessarily isotropic*) likelihood. In our view, presenting a new, scalable, way to perform posterior inference in GPs is a valuable contribution on its own which can pave the path for future work on hyperparameter selection and non-conjugate inference.\n\n\n* *The fixed hyperparameter* setting occurs in large-scale Bayesian optimization, particularly in industrial settings. Here, hyperparameters are often selected using historical (offline) data, and no updates are made online thereafter. Combining online hyperparameter updates and closed-loop systems is very challenging in practice and rigorous theory of these updates is scarce.\n\n\n* *Ill conditioned systems* appear in almost all all moderate-to-large-scale GP problems. See Foster et al. (2009, JMLR), \"Stable and Efficient Gaussian Process Calculations\", or Terenin et al. (2023, arXiv:2210.07893), \"Numerically Stable Sparse Gaussian Processes via Minimum Separation using Cover Trees\". Regularization through for instance inducing point choice helps, but at the cost of bias and performance. Our work suggests an orthogonal way to handle instability is to design algorithms that tolerate ill-conditioning well.\n\n\n\n2. **Benign non-convergence**  \n\n* > “That SGD converges slowly in extrapolation region is a bit concerning. Especially in applications like Bayesian optimization...high interest for exploration.\"  \n\n* This is a very good comment; it illustrates why we find this work so exciting: a priori, we expected the same thing.  \nOur empirical results instead showed that SGD can achieve strong performance in spite of non-convergence. In particular, SGD produces error bars in the extrapolation region which are closer to the prior than the true GP, and thus SGD overestimates uncertainty here. This may cause over-exploration, making convergence somewhat slower. We consider this \"benign\" compared to underestimating uncertainty which may cause catastrophic failure in Bayesian optimization (convergence to a local optimum).  \n\n* The distribution-shift setting is hard to make reasoned arguments about, it is a rather ill-posed problem and there is no guarantee that the exact Bayesian model will perform best in such a setting.\n\n\n\n\n3.  **Formal characterization of region-specific error.**  \nA full mathematical characterization of the part of state-space where non-convergence occurs is what we initially aimed for; however, it proved too difficult. Such a characterization would require one to understand where in space do eigenfunctions corresponding to intermediate eigenvalues occur, which is non-trivial because it is a non-asymptotic question. We believe the amount of work required for this would warrant a separate submission. Alike the reviewer, we think the current analysis is insightful.\n\n4. **The suggested extensions** are good ideas! *Non-conjugate inference* can be immediately achieved via the Laplace approximation as in Antorán et al. (2023), \"Sampling-based inference for large linear models, with application to linearized Laplace\". *Hyperparameter optimization* would require bi-level optimization: an outer loop for the hyperparameters and/or variational parameters, along with an inner loop for the linear systems. Since we cannot expect the inner loop to converge, one would need to study how to ensure that the outer loop behaves well even if the inner loop is not at the optimum. \n\n5. **Questions on experiments**\n\n* You are correct:  **SVGP shares model hyperparameters with CG and SGD**, and also has some variational parameters where applicable. \n\n* We would expect **CG run for sufficiently long** to outperform all alternatives on all data sets. However, for houseelectric (2M+ points), this might take weeks. We are unable to commit that much compute. \n\n\n* **We do provide \"RMSE to exact GP\"** on our four small datasets in figures 3 and 8. In Table 1 we use test RMSE and NLL since for datasets with more than 50k observations (the focus of our work), we cannot do exact GP inference. In Table 1, for the four smallest datasets, CG converges to $10^{-2}$ tolerance and thus can be thought of as an **exact GP baseline**. \n\n\n\n6. **Question: CG non-monotonicity**.  \nThis is a good question which we also asked ourselves for some time. Our best explanation: CG converges monotonically in the RKHS norm induced by the chosen kernel. For the Matérn kernel, this RKHS norm is (effectively, i.e. with certain parameter choices and up to norm equivalence) a weighted sum of the $L^2$ norms of the first $m$ derivatives plus the $L^2$ norm of the function itself. CG is thus trading off minimizing the norm of the $m$ derivatives of the function at the expense of the 0th order term: the $L^2$ error in the fit. The derivative norms being minimized first yields the divergence when looking at just $L^2$ error. Of course, since the RKHS norm is a sum of said $m+1$ non-negative terms, minimizing it will eventually force the 0th order term (error in $L^2$ norm) to go to zero too.\n\n\n7. **Other**. Thanks for the suggestions! We have reworked Fig 4 and Proposition 1 to address your comments and add all required definitions.\n\n----\n\nWith these responses in mind, we gently and politely request that you please consider increasing your score towards firm acceptance."}, {"review_id": "UpQI6WS5Md", "reviewer": "Reviewer_w7c5", "summary": "The paper presents a novel approach for sampling from GP posteriors based on SGD, which bypasses the need to solve the typical linear system of equations that is prevalent in both the exact variant (cubic in the number of query points) or the pathwise approximation (cubic in the number of data points). The SGD approximation is provided for both exact GPs and for inducing point approximations. Moreover, the approximation quality is investigated for three regions of varying data density, and explanations for the (occasionally poor) approximation quality are given.", "questions": "- __Comparison to Pathwise Sampling:__ Does the proposed method hold any advantages to pathwise sampling in a low-to-moderate data regime in terms of accuracy or complexity? At which point (in #data points) does SGD start becoming beneficial?\n- __Computing the predicitve uncertainty:__ Perhaps a trivial question, but how is the predicitve uncertainty computed when one can only access the posterior mean and samples from the posterior (and not the posterior variance) in Section 4.1?", "limitations": "Limitations are adequately addressed.", "rating": 9, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- __Novel, interesting idea:__ The use of SGD for atypical objectives is interesting, and the analysis of the approach is detailed from both a theoretical and empirical perspective.\n\n- __Clearly addressed limitations:__ The SGD approach is _not_ a silver bullet, and the authors make this clear by highlighting the approximation quality in data-dense regions (good), faraway regions (good), and interpolation regions (not as good). \n\n- __Informative, well-designed figures:__ The various figures are not only visually appealing, but provide . Figure 1 and 4 in particular highlight the strengths and shortcomings of the method nicely, and seamlessly add intuition as to why that is. \n\n- __Diverse Experiments:__ Experiments from both large-scale GP regression and BO are included, which demonstrates that the method is applicable and potent in both domains.\n\n- __Clear writing:__ The paper is consistently well-formulated, pedagogical and as far as I could tell, correct. Moreover, I beieve that there has been substantial effort to provide the reader with additional intuition for why the approach is effective.", "weaknesses": "I struggle to find weaknesses with the paper, but remain unconvinced on its potential impact due to the relatively small niche (sampling for GPs in large data regime) that is addressed. However, I am not overly confident in this assessment, and invite the authors to challenge my opinion on this topic. For example, do the authors see opportunities for impactful follow-up work which spans other areas of ML?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank you their time in reading our work! We are thrilled that you found our writing and plots “Clear\" and \"Informative” and agreed that our paper \"Clearly addressed limitations\". We now discuss weaknesses and respond to your questions:\n\n\n**Significance:**\n> \"I ... remain unconvinced on its potential impact due to the relatively small niche (sampling for GPs in large data regime) that is addressed.\"\n\n* **Our main motivating setting is large-scale Bayesian optimization**. Sampling from large-data GPs with fixed hyperparameters is a key component in large-scale Bayesian optimization, particularly in industrial settings. Bayesian optimization (whether under this name, or that of GP bandit algorithms) is a strong approach for the optimization of black-box systems. In particular, sampling from GPs is a core element of the Thompson sampling algorithm. Historically, due to the cost of fitting GPs, Bayesian optimization was limited to small—perhaps even toy—systems; however, work undertaken over the last 5-10 years on scaling inference in GPs, which our method is a part of, now allows for its use on an industrial scale (e.g., optimizing stock levels and recommendations at Amazon). The Thompson-sampling-based approach to Bayesian optimization is particularly well suited to parallelization and asynchronous processing, and thus combines well with such large-scale systems. We believe our approach to sampling to be particularly well-suited for Thompson sampling and easy to use due to its robustness to ill-conditioning. Thus, it has the potential to be adopted by industry users at the usual big-name tech companies (which all provide online recommendations to users and tackle other similar bandit problems, well-suited for Bayesian optimization). \n\n* More generally, there is growing interest in applying GPs to **spatiotemporal modeling** (Howes et al., PLOS Global Public Health 2023, \"Spatio-temporal estimates of HIV risk group proportions for adolescent girls and young women across 13 priority countries in sub-Saharan Africa\"), applications in the **physical and natural sciences** (Goḿez-Bombarelli et al., ACS Central Science 2018, \"Automatic Chemical Design Using a Data-Driven Continuous Representation of Molecules\") and **climate modeling** (Thompson et al., Environmental Data Science 2022, \"A dependent multimodel approach to climate prediction with Gaussian processes)\". Here datasets tend to be large, and as shown by Foster et al. (2009, JMLR, \"Stable and Efficient Gaussian Process Calculations\"), and Terenin et al. (2023, arXiv:2210.07893, \"Numerically Stable Sparse Gaussian Processes via Minimum Separation using Cover Trees\"), **ill-conditioned systems appear in almost all moderate-to-large-scale GP problems**. Regularization through for instance inducing point choice helps, but at the cost of bias and performance. Our work suggests an orthogonal way to handle instability is to design algorithms that tolerate ill-conditioning well.\n\n\n**Questions:**\n\n1. **Comparison to Pathwise Sampling:**  \n> \"Does the proposed method hold any advantages to pathwise sampling in a low-to-moderate data regime in terms of accuracy or complexity? At which point (in #data points) does SGD start becoming beneficial?\"\n\n* This is a great question! Since our approach is an approximation to efficient sampling (namely pathwise conditioning with no approximations except for the prior term) we expect it to perform worse whenever solving the involved linear systems exactly is tractable, for instance in the low-data, well-conditioned regime. From our experiments (Table 1), we estimate that the transition where SGD may start to become better than CG occurs in the 50k-100k datapoint range. Where in that range depends on kernel-matrix conditioning. For very poorly conditioned kernel matrices, SGD can perform better with only a couple of tens of thousands of points.  \n\n2. **Computing the predictive uncertainty:**  \n> \"How is the predictive uncertainty computed when one can only access the posterior mean and samples from the posterior (and not the posterior variance) in Section 4.1?\"\n\n* For each test-point, we estimate the scalar predictive variance from 64 0-mean posterior samples $f_i(x)$ as $\\frac{1}{64}\\sum_{i=1}^{64} f_i(x)^2$. We do this for all methods under consideration. This is tucked away in the second paragraph of Section 4.1 - thank you for pointing this out, we will look into making this point easier to find."}, {"review_id": "oa8k5B8cf6", "reviewer": "Reviewer_xja3", "summary": "I thank the authors for supplying the additional experiments and making the comparisons. It appears the technique is a novel and competitive method that can perform very well even for difficult datasets in the large-scale regime. \n\nBelow is the initial review, unchanged.\n-------------------------------------------------------------------\n\nThe authors propose to sample from the posterior of a GP via optimisation using SGD. The authors rephrase sampling as a quadratic optimisation problem that allows an efficient approximation of the gradient using random kernel features. Then, SGD is used to find the solution. Theoretical derivations show tht the proposed algorithm gives good variance estimates in densely sampled areas as well as areas far outside the sampled regions. Further, the algorithm performs better than CG and SVGP in certain settings.", "questions": "Question: Can you provide plots including a parameter study of SVGP including a stronger optimizer?\n\nSuggestion: I would suggest to introduce SoP method in the related work, but also compare them to the proposed method in the experimental section. I would suggest to compare them empirically to the proposed method by given them a dataset size with computation budget the same as the proposed method to allow a fair comparison on wallclock time. ", "limitations": "-", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The article provides a reltive complete package with an understandable derivation, good theoretical results and experimental evaluations.\n\nThe observations of the good performance of SGD in densely sampled regions is very interesting.", "weaknesses": "The article omits subset of data (SoD) methods completely, both in the related work and the experimental section. [1] introduces them as as category, [2] establish them empirically as competitive method in sampling and [3] investigates the size of the subsampling dataset.\n\n[1] A Unifying View of Sparse Approximate Gaussian Process Regression, Qui~nonero-Candela and Rasmussen (2005) \n\n[2] A Framework for Evaluating Approximation Methods for Gaussian Process Regression, Chalupka and Murray and Williams (2013)\n\n[3] Adaptive Cholesky Gaussian Processes, Bartels et al. (2023)\n\nSubsampling methods perform similarly to the proposed method as they work well in densely sampled regions and also well in regions far away from the data distributions with only region of large error in the low sampled tails. This makes comparison mandatory.\n\nThe compairosn with SVGP is unfair due to the small number of inducing points (1024) and a bad optimisation algorithm for them (ADAM). It seems like the authors did not tune the competing methods well, making the resulting baselines weak. \n\nTheoretically, the authors omit in their complexity class evaluations the number of SGD steps, which might become very large. \nFurther, the authors omit the standard work on GP, while refering to standard notation introduced by it:\n\n[4] Gaussian Processes for Machine Learning, Carl Edward Rasmussen and Christopher K. I. Williams (2006)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review! Below we address the key points:\n\n**Comparison with Subset of Data Methods**\n\nThank you for bringing SoD methods to our attention. *We ran SoD on our regression datasets.* We randomly select subsets from the training data and build exact GP models with these points. We use the complete datasets for data normalization. We provide mean results and std. err. across dataset splits and subset seeds. \n\nWe consider 5% and 10% subsets for the small and medium datasets. For the larger datasets, we use 25k and 50k (largest possible on an 80GB A100 GPU) subsets. These represent roughly 5% and 10% of 3droad song and buzz, and 1.25 and 2.5% of houselectric.\n\nDue to the character limit, we do not reproduce the numbers from our paper here but we bold best results across this table and the paper. \n\n**SoD Small & Medium Datasets**\n|      |  SoD |        pol       |    elevators    |       bike      |      protein     | keggdir   |\n|:----:|:----:|:----------------:|:---------------:|:---------------:|:----------------:|-----------|\n|      |      |     N = 15000    |    N = 16599    |    N = 17379    |     N = 45730    | N = 48827 |\n| RMSE |   5%  | 0.20 ± 0.01  | 0.44 ± 0.00 | 0.18 ± 0.01  | 0.71 ± 0.01 | 0.11 ± 0.00  |\n|      |  10%  | 0.16 ± 0.00  | 0.41 ± 0.01 | 0.13 ± 0.00  | 0.66 ± 0.00 | 0.10 ± 0.00  |\n|  NLL |   5%  | -0.43 ± 0.02 | 0.57 ± 0.01 | -0.61 ± 0.09 | 0.99 ± 0.01 | -0.76 ± 0.10 |\n|      |  10%  | -0.66 ± 0.02 | 0.51 ± 0.01 | -1.28 ± 0.07 | 0.91 ± 0.01 | -0.82 ± 0.09 |\n\n**SoD Large Datasets**\n|      |  SoD |      3droad      |       song      |       buzz      |     houseelec    |\n|:----:|:----:|:----------------:|:---------------:|:---------------:|:----------------:|\n|      |      |    N = 434874    |    N = 515345   |    N = 583250   |    N = 2049280   |\n| RMSE |  25k | 0.23 ± 0.00  | 0.82 ± 0.00 | 0.33 ± 0.00 | 0.06 ± 0.01  |\n|      |  50k | 0.16 ± 0.00  | 0.81 ± 0.00 | **0.32 ± 0.00** | **0.06 ± 0.00**  |\n|  NLL |  25k | -0.42 ± 0.01 | 1.22 ± 0.00 | 0.30 ± 0.06 | -0.53 ± 0.33 |\n|      |  50k | **-0.78 ± 0.01** | **1.20 ± 0.00** | **0.26 ± 0.06** | -0.49 ± 0.39 |\n\nIn terms of mean prediction, SoD does not perform best on any small or medium datasets. We also tried higher percentages: the results are similar. On the large datasets, the 50k subset performs best on buzz, where SVGP previously was the best method, and on NLL for 3droad and Buzz. Note that an exact GP on 50k points requires 40 GB of GPU memory; it is not accessible to most practitioners.\n\nIn terms of error bar geometry, SoD performance is heavily dataset dependent: it performs well when the observations are largely redundant. **We graphically illustrate this with Figure 1(r) in the rebuttal PDF attached to the summary post**.\n\nTo conclude, we would like to emphasize that our contribution is presenting a **novel approach for scaling up Gaussian processes**, quite dissimilar from existing ones. This is valuable due to its opportunity to unlock new avenues for GP research, not because our technique outperforms all baselines across all tasks, which it does not - and, *this is not one of our claims*, as Table 1 (in paper) clearly shows instances where SGD is outperformed by CG and SVGP.\n\n\n**SVGP Hyperparameters and optimization algorithm (ADAM)**\n\n* **Hyperparameters.** To facilitate comparisons with prior work, **our paper used the same SVGP and CG  hyperparameters as Wang et. al. 2019** (a very well-cited NeurIPS2019 paper), which studies conjugate-gradient-based Gaussian processes, and whose techniques are now used in GPyTorch. \n\nWe have also trained SVGP with $M=4096$ inducing points on our 4 largest datasets (we will include the rest in the camera ready), which is rather expensive because 10k optimization steps are needed for convergence and each step's cost is cubic in the number of inducing points $\\mathcal{O}(M^3)$. \n\n|   SVGP 4096   |      3droad     |       song      |       buzz      |     houseelec    |\n|:----:|:---------------:|:---------------:|:---------------:|:----------------:|\n| RMSE | 0.49 ± 0.01 | 0.81 ± 0.00 | **0.33 ± 0.00** |  0.11 ± 0.01 |\n|  NLL | 0.67 ± 0.02 | 1.22 ± 0.00 | **0.25 ± 0.04** | -0.90 ± 0.10 |\n\nAlthough increasing the inducing points improves SVGP's performance on all large datasets, SVGP 4096 only performs best on buzz, where the 1024 point version was already best.\n\n* **Adam.** Our understanding is that ADAM is currently the best optimizer for SVGP. It is the default in GPyTorch, GPFlow, and GPJax. We suspect the reviewer may be thinking of *full-batch* methods such as LBFGS, which **cannot be used with stochastic approximation (i.e. minibatching)**. These optimizers are SoTA for full-batch inducing point approaches such as Titsias 2009, but are incompatible with SVGP (namely, *stochastic* variational GP) which is mini-batch based.\n\n**SGD steps in the complexity analysis**\n\nThis is a good question! Our work's key finding is that *SGD does not need to be run to convergence to produce good empirical performance*. On this basis, it makes sense to view the number of SGD steps as a hyperparameter. This is reflected in our paper's language; **We make no claims about the complexity of SGD Inference, only about the cost of a single optimization step** - e.g. in line 106: “(7)[our unbiased estimator] presents O(N) complexity, in contrast with the O(N^2) complexity of one CG step.”\n\n**Empirically, the number of steps needed to obtain a given performance level is, roughly independent of the dataset size**. We use 100k SGD steps across all experiments, except the Bayesian optimization experiment where we vary this parameter as part of the experiment. We find not only that this is sufficient to obtain reasonable results in all cases but also that the convergence plots (Figures 3, 5, 8) present a similar shape for numbers of observations running from 10k to 2M. \n\n**We have added citation of Rasmussen and Williams (2006)**"}]

Author rebuttal:
We thank all reviewers for the time taken to read our paper and for their insightful and helpful comments. We are pleased the reviewers unanimously found our paper to be well-written, novel, and interesting.

----

The two most pressing concerns come from:

* Reviewers Hc1f and w7c5 ask about the **significance of the setting covered**.



The technique presented is immediately useful for **large-scale Bayesian optimization**. This problem appears most often in industrial settings and is the focus of our paper's second experiment

Our work also presents a novel approach to deal with **ill-conditioning in GP kernel matrices**. Ill-conditioned systems appear in almost all moderate-to-large-scale GP problems. 

Further details on both of these points are given in the rebuttals for Reviewers Hc1f and w7c5.

----

* Reviewer xja3 is concerned about our **lack of comparison against Subset of Data (SoD) inference and our choice of SVGP hyperparameters**.

We have run **additional regression experiments using SoD methods**. 

We further illustrate the strengths and weaknesses of SoD in **Figure 1(r) in the attached PDF**: SoD performs strongly only when the data is very redundant.
We have also **run the SVGP baseline with the number of inducing points increased to 4096.**

Quantitative results are provided in the individual response to xja3.

Meta-review: This paper considers the use of stochastic gradient descent to approximately solve linear systems in Gaussian process (GP) modeling tasks, and this method is explored in a number of experiments for prediction and decision making. Overall, reviewers found the contribution to be novel and leading to interesting insights, and reviewers also appreciated presentation of the ideas and clean figures. Understanding a new approach of training GP models in large data settings is a topic of interest to a broad audience from researchers interested in methodology/theory to ones interested in applying scalable GP methods in practice.  

We encourage the authors to revise the paper to incorporate suggestions from the reviewers and points made during the discussion phase, e.g., emphasizing the connection to broader areas where this work would be impactful, including the additional baselines presented during the rebuttal, and adding more in depth discussion of limitations and potential follow-up directions.
Decision: Accept (oral)