Computation-Aware Gaussian Processes: Model Selection And Linear-Time Inference
Authors: Jonathan Wenger
Kaiwen Wu
Philipp Hennig
Jacob R. Gardner
Geoff Pleiss
John Patrick Cunningham
Venue: NeurIPS 2024 poster (2024)
Keywords: Gaussian Processes
Model Selection
Approximate Inference
Variational Inference
Probabilistic Numerics
Primary area: probabilistic_methods
TL;DR: We demonstrate how to perform model selection for computation-aware Gaussian processes enabling training on over a million data points on a single GPU.

Abstract: Model selection in Gaussian processes scales prohibitively with the size of the training dataset, both in time and memory.
While many approximations exist, all incur inevitable approximation error.
Recent work accounts for this error in the form of computational uncertainty, which enables---at the cost of quadratic complexity---an explicit tradeoff between computational efficiency and precision.
Here we extend this development to model selection, which requires significant enhancements to the existing approach, including linear-time scaling in the size of the dataset.
We propose a novel training loss for hyperparameter optimization and demonstrate empirically that the resulting method can outperform SGPR, CGGP and SVGP, state-of-the-art methods for GP model selection, on medium to large-scale datasets.
Our experiments show that model selection for computation-aware GPs trained on 1.8 million data points can be done within a few hours on a single GPU.
As a result of this work, Gaussian processes can be trained on large-scale datasets without significantly compromising their ability to quantify uncertainty---a fundamental prerequisite for optimal decision-making.

Computation-Aware Gaussian Processes: Model Selection And Sinear-Time Inference
Jonathan Wenger'
Kaiwen Wu²
Philipp Hennig
Jacob Gardner"
Geoff Pleiss
John P. Cunningham
Columbia University
University of Pennsylvania
University of Tübingen, Tübingen AI Center University of British Columbia, Vector Institute
Abstract
Model selection in Gaussian processes scales prohibitively with the size of the training dataset, both in time and memory. While many approximations exist, all incur inevitable approximation error. Recent work accounts for this error the form of computational uncertainty, which enables-at the cost of quadratic complexity-an explicit tradeoff between computational efficiency and precision. Here we extend this development to model selection, which requires significant enhancements the existing approach, including linear-time scaling in the size of the dataset. We propose a novel training loss for hyperparameter optimization and demonstrate empirically that the resulting method can outperform SGPR, CGGP and SVGP, state-of-the th-tthartt art methods for GP model selection, on medium to large- scale datasets. Our experiments show that model selection for computation-awe aware GPs trained on 1.8 million data points can be done within few hours single GPU. As result of this work, Gaussian processes can be trained on large-scale datasets without significantly compromising their ability to quantify uncertainty- a fundamental prerequisite for optimal decision-making.
Introduction
Gaussian Processes (GPs) remain a popular probabilistic model class, despite the challenges in scal- ing them to large datasets. Since both computational and memory resources limited in practice, approximations are necessary for both inference and model selection. Among the many approx- imation methods, perhaps the representation. The resulting posterior approximations typically have a functional form similar to the exact GP posterior, except where posterior mean and covariance feature low-rank updates. This strategy can be explicit-by either defining feature functions (e.g. Nyström [1]. RFF [2])-or lower-dimensional latent inducing point space (e.g. SoR, DTC. FTTC [3], SGPR [4], SVGP [5]), or implicit- using an iterative numerical method (e.g. CGGP [6-10]). All of these methods then compute coefficients for this lower-dimensional representation from the full set of observations by direct projection (e.g. CGGP) or via an optimization objective (e.g. SGPR, SVGP).
common approach is to map the data
lower-dimensional
While effective and widely used in practice, the inevitable approximation enror adversely impacts predictions, uncertainty quantification, and ultimately downstream decision-makinn Many pro- posed methods come with theoretical error bounds [c.g. 2. 11-14], offering insights into the scaling and asymptotic properties of each method. However, theoretical bounds often require too many as- sumptions about the data-generatin genareeinnennggggg process offer "real-world" guarantees [15], and in practice, the fidelity of the approximation is ultimately determined by the available computational resources.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

CholeskyGP
SVGP
CaGP-CG
CaGP-Opt
2-
Vuriance Predictive
0-
0.2
0.0
0.0
0.5
1.0
0.0
0.5
1.0
0.D
0.5
1.0
0.0
0.5
1.0
meon
uncertainty (latent noise)
Training data
for Data Generating Hyperparametees
Figure 1: Comparison of an exact GP posterior (CholeskyGP) and three scalable approximations: SVGP, CaGP-CG and CaGP-Opt (ours), Hyperparameters for each model were optimized using model selection strategies specific to each approximation. The posterior predictive given the data- generating hyperparameters is denoted by gray lines and for each method the posterior (dark-shaded) and the posterior predictive are shown (light-shaded). While all methods, including the exact GP, do not recover the geterrenngg process, CaGP-CG and CaGP-Opt are much closer than SVGP. SVGP expresses almost no posterior variance near the inducing point in the data-sparse region and thus almost all deviation from the posterior mean is considered to be observational noise. In contrast, CaGP-CG and CaGP-Opt express significant posterior variance in regions with no data.
One central pathology is overconfidence, which has been shown to be detrimental in key applica- tions of GPs such as Bayesian optimization [e.g. variance starvation of RFF, 16]. and manifests itself even in maae-oof-the-at variational methods like SVGP. SVGP, because it treats inducing variables as "virtual observations", can overconfident the locations of the inducing points they are not in close proximity to training data, which becomes increasingly likely in higher dimensions. This phenomenon can seen in a toy example in Figure where SVGP has near zero posterior variance at the inducing point away from the data. See also Section S5.1 for a more detailed analysis.
These approximation errors a central issue in inference, but they are exacerbated in model selec- tion, where eтrоrs compound and result in biased selections of hyperparameters [12, 17, 18). Con- tinuing the example, SVGP has been observed overestimate the observation noise [18]. which can lead to oversmoothing. This issue can also be seen in Figure 1, where the SVGP model produces smoother posterior mean than the exact (Cholesky)GP and attributes most variation from the pos- terior mean to observational noise (see also Figure S3(b)). There have been efforts to understand these biases [18] and to mitigate the impact of approximation error on model selection for certain approximations [e.g. CGGP, 12], but overcoming these issues for SVGP remains a challenge.
Recently, Wenger al. [19] introduced computation-avee GP approximation methods which-for a fixed set of pyperparametees-provaaee does not suffer from overconfidence. Like SVGP and the other approximations mentioned above, CaGP also re- lies on low-rank posterior updates. Unlike these other methods, however, CaGP's posterior updates are constructed to guarantee that its posterior variance is always larger than the exact GP variance. This conservative estimate can be interpreted as additional uncertainty quantifying the approxima- tion error due limited computation; i.e. computational uncertainty. However, so far CaGP has fallen short in demonstrating wallclock time improvements for posterior inference over variational methods and model selection has so far remained an open problem.
Gaussian processes (CaGP), class of
Contributions In this work, we extend computation-aware Gaussian processes by demonstrating how to perform inference in linear time in the number of training data, while maintaining its the- oretical guarantees. Second, we propose a novel objective that allows model selection without a significant bias that would arise from naively conducting model selection on the projected GP. In detail, we enforce a sparsity constraint on the "actions" of the method, which unlocks lincar-time inference, in a way that is amenable hardware acceleration. We optimize these actions end-to-end alongside the hyperparameters with respect to custom training loss, to optimally retain as much information from the data as possible given a limited computational budget. The resulting hyperpa- rameters are less prone to oversmoothing and attributing variation to observational noise, as can be seen in Figure 1, when compared to SVGP. We demonstrate that our approach strongly competi- tive on large-scale data with sate-of-thhe-rt variational methods, such as SVGP, without inheriting their pathologies. As a consequence of our work, one can train GPs on up to 1.8 million data points in few hours on single GPU without adversely impacting uncertainty quantification.

---

2 Background
We aim learn a latent function mapping from X C Rd to Y C IK given training dataset X (I) E Raxd of n inputs 2j E IR and corresponding targets y = (yI,...; W.) E R"
Gaussian Processes A Gaussian process f (),,, Ko) is stochastic process with mean func- tion j and kernel Kg such that f = f(x)= = with mean μ4 = (u(z;) and covariance K; (z x,). The kemel Kg depends on hyperparam- eters 0 I2, which we omit in our notation. Assuming y f(X) ~ (f
~ ([,, Ko) is jointly Gaussian
o² I), the posterior is Gaussian process GP(,,, K,) where the mean and covariance functions evaluated test input To E I& are given by
Ne (zxo)) - m(xz)) + K X X)u,,
(1)
K. (f(o))) (x., xo) K (Io. X) X)K
K (X,
where K = K + and the representer weights are defined
K
(y m))
In model selection, the computational bottleneck when optimizing kernel hyperparameters the repeated evaluation of the negative log-marginal likelihood
NLLL (0) - -oggg(y 0) ( - 11) k p1) + log det (K) +n log (2g(2)))
quadratic loss: promotes data fit model complexity
and its gradient. Computing (2) and its gradient via a Cholesky decomposition has time complexity Onn*) and requires (n²)) memory, which prohibitive for large n.
Sparse Gaussian Process Regression (SGPR) [4] Given a set of m << n inducing points Z = (Z1....., I and defining u f(Z) (f( f(mm))) SGPR defines varia-
tional approximation to the posterior through the factorization (f(()) 1 ~ a((()) = p((())
u) qu))du, where ( is m-dimensional multivariate Gaussian. The mean and covariance of q(u) (denoted as m == E,(2), : Σ == Cov, Co))(u)) are jointly optimized alongside the kernel hyperpa-
rameters 0 using the evidence lower bound (ELBO) as an objective:
m, Σ, 0, Z arg minm, E,0,0 2. (ELBO.
(3)
ELBO m. 2)
= NT. (0) + (g( Il (f v))))
(4)
-- (log ply f)) + KL(((u) p(u)) 2 logg(y |)..
The inducing point locations Z can be either optimized as additional parameters during training or chosen a-priori, typically in data-dependent way [see e.g. Sec. 7.2 of 20). Following Titsias [4], ELBO optimization and posterior inference both require (nn²² computation and ((m)) memory. significant improvement the costs of exact GPs.
Stochastic Variational Gaussian Processes (SVGP) [S] SVGP extends SGPR to reduce com- plexity further to O(m²) computation and O(m²) memory. It accomplishes this reduction by placing the first term in Equation (4) with an unbiased approximation
En (log p f)) = E,,r, (Ei P(u f(xx)) = n (ey((()))) (log f(x)))).
Following Hensman et al. [5]. optimize m. Σ alongside 0 Z through joint gradient updates. Because the asymptotic complexities no longer depend on n, SVGP can seale to extremely large datasets that would not be able to fit into computer/GPU memory.
Computation-awar Gaussian Process Inference (CaGP) [19] CaGP maps the data y into a lower dimensional subspace defined by its actions S; E Raxt on the data, which defines an approx- imate posterior GP(u, , K:) with
(xx)) = m(xx) + K( X)²; Ki(#,, v.) K (x., x.) K (x., X)C,K(X,
(5)
'Wenger et al. [19] named their computation-awan aware inference algorithm "IterGP", to emphasize its itera- tive nature (see Algorithm S1). We adopt the naming "CaGP" instead, since if the actions S are not chosen sequentially, is more efficient to compute the computation aware posterior non-iteratively (see Algorithm S2).

---

where C, ( KS.) = = K-' is a low-rank approximation of the precision matrix and v; = Ci(y - μ1) = v, approximates the representer weights. Both converge to the corresponding exact quantity as the number of iterations, equivalently the downdate rank, + n. Note that the CaGP posterior only depends on the space spanned by the columns of S., not the actual matrix (Lemma S4). Finally, the CaGP posterior can be computed in (n²i) time and (nii) memory.
CaGP captures the approximation error incurred by limited computational resources as additional uncertainty about the latent function. More precisely, for any data-generating function E I in the RKHS HK- defined by the kernel, the worst-case squared error of the corresponding approxi- mate posterior mean n.l is equal to the approximate predictive variance (see Theorem S2):
sup lell
(u(x²) 14 (xx))) = (2 Io) + +
ME
This guarantee is identical to one for the exact GP posterior mean and variance (see Theorem S1), except with the approximate quantities instead. Additionally, it holds that CaGP's marginal (predic- tive) variance is always larger or equal to the (predictive) variance of the exact GP and monotonically decreasing, i.e. (2 2o) Ký(2z (xto wo) 2 (2,,,, sition S1). Therefore, given fixed hyperparameters, CaGP is guaranteed to never be overconfident and as the computational budget increases. the precision of its estimate increases. Here we call such a posterior iittttattio--aare and we will extend the use of this object to model selection.
K. (30, Xo) for i S j S n (see Propo-
3 Model Selection in Computation-Aware Gaussian Processes
Model selection for GPs most commonly entails maximizing the evidence log 0) a function of the kernel hyperparameters E IRP. As with posterior inference, evaluating this objective and its gradient is computationally prohibitive in the large-scale setting. Therefore, our central goal will be to perform model selection for computation-awe -aware Gaussian processes in order to scale to a large number of training data while avoiding the introduction of pathologies via approximation.
We begin by viewing the computation-aware posterior as exact inference assuming we can only observe i linear projections of the data defined by (linearly independent) actions S, E Raxi i.e.
ÿ = S y E R', where 5 S, chol(S," (.))- Rxxi
is the action matrix with orthonormalized columns. The corresponding likelihood is given by p( f(X)) N (ū: S f(X),o²I)
(8)
As we show in Lemma S1. the resulting Bayesian posterior is then given by Equation (5). Recall that the CaGP posterior only depends on the column space of the actions S, (see Lemma S4), which is why Equation (5) can be written in terms of S; directly rather than using SY.3
Projected-Data Log Marginal Likelihood The reinterpretation of the computation-awe aware poste- rior as exact 1 Bayesian inference immediately suggests using evidence maximization for the projected data ÿ = 5 y E R' as the model selection criterion, leading to the following loss (see Lemma S2):
NNL (0) = log p 0) = - log / ÿ S. ,4 M. S KS)) proj
penalizes ncar-colincar actions
2 (y- 14) (S! KS.) s7 -s. (y 11) + log dSt K - log det S.) +i logg(222) (2T))
quadratic loss: promotes fitting projected data 9
projected model complexity
Equation (9) involves Gaussian random variable of dimension i < n, and so all previously in- tractable quantities in Equation (2) (i.e. the inverse and determinant) are now cheap to compute. Analogous to the CaGP posterior, we can express the projected-data log marginal likelihood fully in terms of the actions S, without having to orthonormalize, which results in an additional term penal- izing colinearity. Unfortunately, this training loss does not lead to good generalization performance, as there is only training signal in the i-dimensional space spanned by the actions S, the data are projected onto. Specifically, the quadratic loss term only promotes fitting the projected data ÿ. not all observations y. See Figures SI and S2 for experimental validation of this claim.
? At first glance, SVGP satisfies similar result (Theorem S3). However, this statement does not express the variance in terms the worst-case squared error to the "true" latent function. See Section S1.2 for details. è GGiven this observation, we sometimes abuse terminology and refer to the actions as "projecting" the data a dimen-dima space, although S: does not need to have orthonormal columns.

---

ELBO Training Loss Motivated by this observation, we desire tractable training loss that en- courages maximizing the evidence for the entire set of targets y. Importantly though, we need to accomplish this evidence maximization without incurring prohibitive ((n)) computational cost. We define a variational objective, using the computation-awan
posterior q y. O) in Equation (5) define a variational family 2 == { = (f: (X), (X, X)) S E [2*x) j parametrized
by the action matrix S. We can then specify a (negative) evidence lower bound (ELBO) as follows:
PELBO CaGP (
(N)
KL P.). 2 log p(y 0).
(10)
halances and model complexity regularizes s.t. 9: R2 p-
This loss promotes learning the same hyperparameters as if we were to maximize the computa- tionally intractable evidence logr(y 0) while minimizing the error due posterior approximation (f) = (f 0). In the computation-awalt aware setting, this translates to minimizing computational uncertainty, which captures this inevitable error.
Although both the evidence and KL terms of the ELBO involve computationally intractable terms, these problematic terms cancel out when combined. This results in an objective that costs the same as evaluating CaGP's predictive distribution, i.e.
ELBO CaaP
(0)
= 1 2 02 (1y
M4 (X)1
+
Σ k,(23, (n (n--)) log((?) nlgg(22))
==1
(11)
+ -o.5KSS;
KSü; - tr((5) - KS) + log det KS) - log det (ST $)
1)) For a derivation of this expression of the loss see Lemma S3. If we compare the training loss PELBO CCGP in Equation (10) with the projected-data NI.L in Equation (9), there is an explicit squared loss penalty on the entire data y, rather than just for the projected data ý. resulting in better generalization as Figures S1 and $2 show on synthetic data. In our experiments, this objective was critical to achieving tfat-oof-the-ar performance.
where é; o; = (ST
4 Choice of Actions
So far we have not yet specified the actions S E Rxx mapping the data to a lower-dimensional space. Ideally, we would want to optimally compress the data both for inference and model selection.
Posterior Entropy Minimization We can interpret choosing actions S as a form of active learn- ing, where instead of just observing individual datapoints, allow ourselves to observe linear combinations of the data y. Taking an information-theercee to choose actions such that uncertainty about the latent function is minimized. In fact, show in Lemma S5 that given a prior ff ~ GP( K) for the latent function and a budget of i actions S, the actions that minimize the entropy of the posterior at the training data
perspective [21], we would then aim
($1.... 8,) = arg min H () SER=xx
(12)
are the top-i cigenvectors 81...., 8 of K in descending order of the eigenvalue magnitude (see also Zhu et al. [22]). Unfortunately, computing these actions is just as prohibitive computationally as computing the intractable GP posterior.
(Conjugate) Gradient / Residual Actions Due to the intractability of choosing actions to min- imize posterior entropy, we could try to do so approximately. The Lanczos algorithm [23] is an iterative method to approximate the cigenvalues and eigenvectors of symmetric positive-definite matrices. Given an appropriate seed vector, the space spanned by its approximate eigenvectors is equivalent to the span of the gradients residuals r; = (y p1) Kv,, of the method of Conjugate Gradients (CG) [24] when used to iteratively compute an approximation i. e, K K- 11) to the representer weights V.. We show in Lemma S4 that the CaGP posterior only depends on the span of its actions. Therefore choosing approximate eigenvectors computed via the Lanezos process as actions is equivalent to using CG residuals. This allows us to reinterpret CaGP-CG, as introduced by Wenger et al. [19]. as approximately minimizing posterior entropy.* See Section S3.1 for details.
'Wenger et al. [Sec. 2.1 of 19] showed that CaGP using CG residuals as actions recovers Conjugate- Graadeen-bbased GPs (CGGP) [7, 9, 10] in its posterior mean, extending this observation to CGGP.
5

---

Action S
Action Sa
Action S3
Action S₁
Action Ss
CaGp-0G Ligenvectors
Grassman distance di.S Vownn
0.0
1.0
Method
Random
CaGP-CG
CaGP-Opt
4-
1.0
CaCP-Opt
100
Epoch
Figure 2: Visualization of action vectors defining the data projection. We perform model selection using two CaGP variants, with CG and learned sparse actions denoted as CaGP-CG, and CaGP- Opt-on toy 2-dimensional dataset. Left: For each Ø2 E .t.......nn. we plot the magnitude of the entries of the top-5 cigenvectors of K and of the first five action vectors. Yellow denotes larger magnitudes; blue denotes smaller magnitudes. Right: We compare the span of the actions S against the top-i eigenspace throughout training by measuring the Grassman distance between the two subspaces (see also Section S5.2). CaGP-CG actions are closer to the kernel eigenvectors than the CaGP-Opt actions, both of which are more closely aligned than randomly chosen actions.
As Figure 2 illustrates, CG actions are similar to the top-i eigenspace all throughout hyperparameter optimization. However, this choice of actions focuses exclusively on posterior inference and incurs quadratic time complexity Onni).
So far in action choices have entirely ignored model selection and tried to choose optimal actions assuming, fixed kernel hyperparameters. The second contribution of this work, aside from demonstrating how perform model selection, is recognizing that the tions should be informed by the outer optimization loop for the hyperparameters. We thus optimize the actions alongside the hyperparameters end-to-end, meaning the training loss for model selection defines what data projections are informative. This way the actions are adaptive to the hyperpa- rameters without spending unnecessary budget on computing approximately optimal actions for the current choice of hyperparameters. Specifically, the actions are chosen by optimizing BELBO function of the hyperparameters and the actions S, s.t.
Learned Sparse Actions
CCaGP
(0., S,) arg min (..) e CaGP
EELBO (
(13)
Naively this approach introduces an n dimensional optimiza- tion problem, which in general is computationally prohibitive. To keep the computational cost low and to optimally leverage hardware acceleration via GPUs, we impose sparse block struc- ture on the actions (see Eq. 14) where each block column vector s R with h n/i entries such that the to- tal number of trainable action parameters, i.e. non-zero entries nnz(S) k.
0
0
8
S
(14)
0
equals the number of training data. Due to the sparsity, these actions cannot perfectly match the maxi- mum eigenvector actions. Nevertheless, conjunction with hyperparameter optimization produces a nontrivial alignment with optimal action choice minimizing posterior entropy. Importantly, the sparsity constraint not only reduces the di- mensionality of the optimization problem, but crucially also reduces the time complexity of posterior inference and model selection to linear in the number of training data points.
in Figure 2 that optimizing these sparse actions in
4.1 Algorithms and Computational Complesity
We give algorithms both for iteratively constructed dense actions (Algorithm S1), as used CaGP- CG, and for sparse batch actions (Algorithm S2), as used for CaGP-Opt, in the supplementary

---

material." 5 The time complexity is (nn²t) for dense actions and Oni max(i, (i, k)) for sparse actions, where k is the maximum number of non-zero entries per column of S, Both have the same linear memory requirement: (ni). Since the training loss ELBO CaCP only involves terms that are also present in the posterior predictive, both model selection and predictions incur the same complexity.
4.2 Related Work
Computational Uncertainty and Probabilistic Numerics All CaGP variants discussed in this paper fall into the category of probabilistic numerical methods [25-27], which aim quantify ap- proximation eттог arising from limited computational resources via additional uncertainty about the quantity of interest [e.g. 28-31]. Specifically, the iterative formulation of CaGP (i.e. Algorithm S1) originally proposed by Wenger et al. [19] employs a probabilisticc linear solver [32-35].
Scalable GPs with Lower-Bounded Log Marginal Likelihoods Numerous scalable GP approxi- mations beyond those in Sections and 2 exist; see Liu et al. [36] for a comprehensive review. Many GP models [e.g., 4, 5, 37-39] learn hyperparameters through maximizing variational lower bounds in the same spirit as SGPR, SVGP and our method. Similar to our work, interdomain inducing point methods [40-42] learn a variational posterior approximation on small set of linear functionals ap- plied to the latent GP. However, unlike our method, their resulting approximate posterior is usually prone to underestimating uncertainty the same manner as SGPR and SVGP. Finally, similar to our proposed training loss for CaGP-CG, Artemev et al. [43] demonstrate how one can use the method of conjugate gradients to obtain a tighter lower bound on the log marginal likelihood.
GP Approximation Biases and Computational Uncertainty Scalable GP methods inevitably introduce approximation error and thus yield biased hyperparameters and predictive distributions, with an exception of Potapczynski et al. [12] which trade bias for increased variance. Numerous works have studied pathologies associated with optimizing variational lower bounds, especially in the context of SVGP [12, 16-18], and various remedies have been proposed. In order to mitigate biases from approximation, several works alternatively propose replacing variational lower bounds with alternative model selection objectives, including leave-one-out cross-validation [44] and losses that directly target predictive performance [45, 46].
5 Experiments
We benchmark the generalization of computation-awant -wware GPs with two different action choices, CaGP-Opt (ours) and CaGP-CG [19], using our proposed training objective in Equation (10) on range of UCI datasets for regression [53]. We compare against SVGP [5], often considered to be state-of-the-are for large-scale GP regression. Per recommendations by Ober et al. [15], we also in- clude SGPR [4] as a strong baseline for all datasets where this is computationally feasible. We also train Conjugate Gradient-based GPs (CGGP) [e.g. 7, 9, 10] using the training procedure proposed by Wenger et al. [10]. Note that CaGP-CG recovers CGGP in its posterior mean and produces nearly identical predictive error at half the computational cost for inference [Sec. 2.1 & 4 of 19]. which is why the main difference between CaGP- CGG and CGGP in our experiments is the training objective. Finally, we also train an exact CholeskyGP on the smallest datasets, where this is still feasible.
Experimental Details All datasets were randomly partitioned into train and test sets using (0.9, 0.1) split for five random seeds. We used a zero-mean GP prior and a Matérn(3/2) 2) kernel with outputscale o² and one lengthscale per input dimension 13, as well as scalar observation noise for the likelihood o², s.t. 0 = (o. la...., I.,, σ) E Rd+2 We used the existing implementa- tions of SGPR, SVGP and CGGP in GPyTorch [7] and also implemented CaGP in this framework (see Section S4.2 for our open-source implementation). For SGPR and SVGP we used m = 1024 inducing points and for CGGP, CaGP-CG and CaGP-Opt we chose i 512. We optimized the hy- perparameters either with Adam [54] for a maximum of 1000 epochs in float32 or with LBFGS [55] for 100 epochs in float64, depending the method and problem scale. On the largest dataset "Power", we used 400 epochs for SVGP and 200 for CaGP-Opt due to resource constraints. For SVGP we used batch size of 1024 throughout. We scheduled the learning rate via PyTorch's [56] LinearLR factoror 1) scheduler for all methods and performed a hyperparameter sweep
'or a detailed analysis see Algorithms SI and S2 in the supplementary material, which contain line-by-line time complexity and memory analyses.
7

---

Table I: Generalization error (NLL, RMSE, and wall-clock time) on UCI benchmark datasets. The table shows the best results for all methods across learning rate sweeps, averaged across five random seeds. We report the epoch where each method obtained the lowest average test NLL, and all per- formance metrics (NLL,, RMSE, and wall-clock runtime) are reported for this epoch. Highlighted in bold and color are the best approximate methods per metric (difference > 1 standard deviation).
Dataset
n
Method
Optim.
L.R
Epoch
NLL
RMSE 4
Avg. Runtime
mean
std
Parkinsons [47]
5 $828
21
CholeskyGP
LBFGS
0.100
100
-3,645
0.002
0.001
0.000
Immn 3s
SGPR
Adam
0.100
268
-2.837
0.087
0.031
0.022
27s
LBFGS
1.000
100
-3.245
0,067
0.007
0.003
2min 14s
SVGP
Adam
0.100
1000
-2.858
0.016
0.006
0.002
2min 25s
CGGP
LBFGS
0,100
81
-2.663
0.141
0.019
0.013
Imin 12s
CaGP.CG
Adam
1.000
250
-2.936
0.007
0.009
0.006
Imin 44s
CaGP-Opt
Adam
1.000
956
-3.384
0,005
0.004
0.002
Imin 27s
LBFGS
0.010
37
-3.449
0.009
0.002
0.000
Imin 53s
Bike [48]
15:422
16
CholeskyGP
LBFGS
100
-3.472
7anin 15s
SGPR
Adam
0.100
948
-2.121
0.110
0.026
0.004
4min 3x
LBFGS
1.000
100
-3.017
0.022
0.009
0.002
4min
SVGP
Adam
0.010
1000
-2.256
0.020
0.020
0.002
6mnn 41s
CGGP
LBFGS
1.000
15
-1.952
0.078
0.024
0.004
2min
CaGP-CG
Adam
1.000
250
-2,042
0.024
0.024
0.002
5min 17s
CaGP-Opt
Adam
1.000
1000
-2.401
0.037
0.018
0.002
8min 10s
LBFGS
1.000
100
-2.438
0.038
0.017
0.001
14min 48s
Protcin [49]
457
9
SGPR
Adam
0,100
993
0.844
0.006
0.561
0.005
10min 25s
LBFGS
0.846
0.006
0.562
6min 56s
SVGP
Adam
0.010
996
0.851
0.006
0.564
0.005
17min 19s
CGGP
LBFGS
0.853
20min 15s
CaGP CG
Adam
1.000
27
0.820
0.006
0.542
0.004
Imiin 26s
CaGP-Opt
Adam
0.100
941
0.829
0.005
0.545
1min 48s
LBFGS
1.000
84
0.830
0.005
0.545
0.004
14min 29.
KEGGu [50]
572488 26
SGPR
Adam
0.100
143
-0.681
0.025
0.123
0.002
2min 4s
LBFGS
1.000
100
-0.712
0.028
0.118
0.003
Smin 58s
SVGP
Adam
0.010
988
-0.710
0.026
0.118
0.003
24min 21s
CGGP
LBFGS
0,100
-0.512
0.034
0.120
0.003
33min 55s
CaGP-CG
Adam
1.000
229
-0.699
0.026
0.120
0.003
39min 5s
CaGP-Opt
Adam
1.000
990
-0.693
0.026
0.120
0.003
22min
LBFGS
0.010
40
-0.694
0.026
0.120
0.003
22min Us
Road [51]
391 387
2
SVGP CaGP-Opt
Adam
0.001
998
0.149
0.007
0.277
0,002
7min 37s
Adam
0.100
1000
-0.291
0.011
0.159
0.003
2h 11min 31s
Power [52]
1 1544352
7
SVGP CaGP-Opt
Adam
0.010
399
-2,104
0.007
0.029
0.000
Sh 7anin 57s
Adam
0.100
200
-2.103
0.006
0.030
0.000
32min 48s
for the (initial) learning rate. All experiments were run on an NVIDIA Tesla V100-PCIE-33GBB GPU, except for "Power", where we used NVIDIA A100 80GB PCIe GPU to have sufficient memory for CaGP-Opt with 512. Our exact experiment configuration can be found in Table SI.
Evaluation Metries We evaluate the generalization performance once per epoch on the test set by computing the (average) negative log-likelihood (NLL) and the root mean squared error (RMSE), as well as recording the wallclock runtime. Runtime is measured at the epoch with the best average performance across random seeds.
CaGP-Opt Matches or Outperforms SVGP Table 1 and Figure show generalization perfor- mance of all methods for the best choice of learning rate. In terms of both NLL and RMSE, CaGP- Opt outperforms or matches the variational baselines SGPR and SVGP at comparable runtime (ex- cept on "Bike"). SGPR remains competitive for smaller datasets; however, it does not scale to the largest datasets. There are some datasets and metrics in which specific methods dominate, for ex- ample on "Bike" SGPR outperforms all other approximations, while on "Protein" methods based on CG, i.e. CGGP and CaGP-CG, perform the best. However, CaGP-Opr consistently performs either best or second best and scales to over a million datapoints. These results are quite remarkable for nu- merous reasons. First, CaGP is comparable in runtime to SVGP on large datasets despite the fact that incurs a linear-time computational complexity while SVGP is constant time. Second, while all of the methods we compare approximate the GP posterior with low-rank updates, CaGP-Opt (with
'While SVGP is arguably linear time since it will eventually loop through all training data, each computation of the ELBO uses constant time stochastic approximation.
8

---

Parkinsons
Bike
Protein
KEGGG
Road
Power
2.5
2.0
2
1.5
Train
2.5
-2
Text NLL
0.5
0.5
0.0
-2
109
10a
1.0
10a
TTest RMSS
0.7
0.4
10
10-2
0.8
10
10
0.2
Zmin Time
Smin
8min 16min
20hin 40min
1h
2h
1.5h
4.5h
Time
Time
Time
CholeskyGP
SGPR
SVGP
OGGP
CaGP-CG
CaGP-OPt
Figure 3: Learning curves of GP approximation methods UCI benchmark datasets. Rows show train and test loss as a function of wall-clock time for the best choice of learning rate per method. CaGP-Opt generally displays "ramp-up" phase early in training where performance is worse than that SVGP. As training progresses, CaGP-Opt matches or surpasses SVGP performance.
i = 512) here uses half the rank of SGPR/SVGP m = 1024. Nevertheless, CaGP-Opt is able to sub- stantially outperform SVGP
spatial datasets like 3DRoad where low-rank posterior updates are often poor [57]. These results suggest that CaGP-Opt can be more efficient approximation than inducing point methods, and that low-rank GP approximations may be more applicable than previously thought [58, 59]. Figure 3 shows the NLL and RMSE learning curves for the best choice of learning rate per method. CaGP-Opt often shows "ramp-up" phase, compared then improves or matches its generalization performance. This gap is particularly large on "Road", where CaGP-Opt initially worse than SVGP but dominates in the second half of training.
SVGP, but
SVGP Overestimates Observation Noise and (Often) Lengthscale In Figure S5 we show that SVGP typically learns larger observation noise than other methods as suggested by previous work syn-
10-1
Coverage Liror
[18, 45] and hinted at by observations thetic data in Figure and Figure S3(b). Addi- tionally larger datasets SVGP also often learns large lengthscales, which in combination with large observation noise can lead ing effect. In contrast, CaGP-Opt generally learns lower observational noise than SVGP. Of course, learning small observation noise, in particular, is important for achieving low RMSE and thus also NLL, and points to why we should expect CaGP- Opt to outperform SVGP. These hyperparameter results suggest that CaGP-Opt interprets more of the data as signal, while SVGP interprets more of the data noise.
10-2
oversmooth-
1.5h
3h
4.5h
Time
SVGP
CaGP-Opt
Figure 4: Uncertainty quantification for CaGP-Opt and SVGP. Difference between the desired coverage (95%) and the empirical cov- erage of the GP 95% credible interval on the "Power" dataset. After training, CaGP-Opl has better empirical coverage than SVGP.
CaGP Improves Uncertainty Quantification Over SVGP key advantage of CaGP-Opt and CaGP-CG is that their posterior uncertainty estimates capture both the uncertainty due to limited data and due to limited computation. To that end, we assess the frequentist coverage of CaGP-Opt's uncertainty estimates. We report the absolute difference between a desired coverage percentage

---

and the fraction of data that fall into the a credible interval of the CaGP-Opt posterior; i.e. coverage Ja Newe D (y T (. Figure 4 compares the 95% coverage emror for both CaGP-Opt and SVGP on the largest dataset ("Power"). From this plot, we see that the CaGP credible intervals are more aligned with the desired coverage. We hypothesize that these results reflect the different uncertainty properties of the methods: CaGP-Opt overestimates posterior uncertainty while SVGP is prone towards overconfidence.
6 Conclusion
In this work, we introduced strategics for model selection and posterior inference for computation aware Gaussian processes, which scale linearly with the number of training data rather than quadrat- ically. The key technical innovations being sparse projection of the data, which balances mini- mizing posterior entropy and computational cost, and scalable way optimize kernel hyperpa- rameters, both of which are amenable to GPU acceleration. All together, these advances enable competitive improved performance over previous approximate inference methods on large-scale datasets, in terms of generalization and uncertainty quantification. Remarkably, our method outper- forms SVGP-often considered the de-facto GP approximation standard- when compressing the data into space of half the dimension of the variational parameters. Finally, we also demon- strate that computation-aware GPs avoid many of the pathologies often observed in inducing point methods, such as overconfidence and oversmoothing.
Limitations While CaGP-Opt obtains the same linear time and memory costs SGPR, it is not amenable to stochastic minibatching and thus cannot achieve the constant time/memory capabilities of SVGP. In practice, this asymptotic difference does not result in substantially different wall clock times, as SVGP requires many more optimizer steps than CaGP-Opt due to batching. (Indeed, on many datasets we find that CaGP-Opt faster.) CaGP-Opt nevertheless requires larger GPUs than SVGP on datasets with more than million data points. Moreover, tuning CaGP-Opt requires choosing the appropriate number of actions (i.e. the rank of the approximate posterior update), though we note that most scalable GP approximations have a similar tunable parameter (e.g. number of inducing points). Perhaps the most obvious limitation is that CaGP, unlike SVGP, limited to GP regression with conjugate observational noise model. We leave extensions to classification and other non-conjugate likelihoods future work.
Outlook and Future Work An immediate consequence of this work the ability to apply approach solves CaGP's
computation-aween aware Gaussian processes to "real-world" problems, open problems of model selection and scalability. Looking forward, an exciting future vision is a general framework for problems involving Gaussian process model with downstream task where the actions are chosen optimally, given resource constraints, to solve said task. Future work will pursue this direction beyond Gaussian likelihoods to non-conjugate models and downstream tasks such as Bayesian optimization.
10

---

Acknowledgments and Disclosure of Funding
JW and JPC are supported by the Gatsby Charitable Foundation (GAT3708), the Simons Foundation (542963), the NSF AI Institute for Artificial and Natural Intelligence (ARNI: NSF DBI 2229929) and the Kavli Foundation. JG and KW are supported by the NSF 145444, DBI-2400135), PH gratefully acknowledges co-funding co-funding by the European Union (ERC. ANUBIS, 101 123955). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. PH a member of the Machine Learning Cluster of Excellence, funded by the Deutsche Forschungsgemminschall (DFG, German Research Foundation) under Germany's ' Excellence Strategy - EXC number 2064/1 - Project number 390727645; he also gratefully acknowledges the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: OIIS 8)399);; and funds from the Ministry of Science, Research and Arts of the State of Baden-W Württemberg. GP acknowledges support from NSERC and the Canada CIFAR Al Chair program.
References
[1] C. Williams and M. Seeger. "Using the Nyström method to speed up kernel machines". In: Advances in Neural Information Processing Systems (NeurIPS). 2001 (cit. on p. 1).
[2] A. Rahimi and B. Recht, "Random Features for Large-Scale Kernel Machines". In: Advances in Neural Information Processing Systems (NeurIPS). 2007 (cit. on p. 1).
[3] J. Quiñonero-Canded and C. E. Rasmussen. "A unifying view of sparse approximate Gaus- sian process regression". .. In: Journal of Machine Learning Research 6 (2005) (cit. on p. 1).
[4] M. Titsias. "Variational learning of inducing variables in sparse Gaussian processes". In: International Conference on Artificial Intelligence and Statistics (AISTATS). PMLR. 2009, pp. 567-574 (cit. pp- 1,3, 7).
[5] J. Hensman, N. Fusi, and N. Lawrence. "Gaussian processes for big data". In: Conference on Uncertainty in Artificial Intelligence (UAI). 2013, pp. 282-290 (cit. on pp. 1, 7). [6] M. Gibbs. "Bayesian Gaussian processes for classification and regression". PhD thesis. 1997
(cit. on p. 1).
[7] J. R. Gardner, G. Pleiss, D. Bindel, K. Q. Weinberger, and A. G. Wilson. "GPyTorch Black- box matrix-matrix Gaussian process inference with GPU acceleration". In: Advances in Neu- ral Information Processing Systems (NeurIPS). 2018, pp. 7576-7586 (cit. on pp. 1,5,7)
[8] G. Pleiss, J. Gardner, K. Weinberger, and A. G. Wilson. "Constant- predictive distribu- tions for Gaussian processes". In: International Conference on Machine Learning (ICML). 2018, pp. 4114-4122 23 (cit. on p. 1).
[9] 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 Neural Information Processing
Systems (NeurIPS) 32 (2019) (cit. on pp. 5, 7).
[10] J. Wenger. G. Pleiss. P. Hennig. P. Cunningham, and J. R. Gardner. "Preconditioning for Scalable Gaussian Process Hyperparameter Optimization". In: International Conference on Machine Learning (ICML). 2022 (cit. on pp. 1, 5, 7).
[11] F. Bach. "On the equivalence between kernel quadrature rules and random feature expan- sions". In: Journal of Machine Learning Research 18.21 (2017), pp. 1-38 (cit. on p. 1).
[12] A. Potapczynski, L. Wu, D. Biderman, G. Pleiss, and Cunningham. "Bias-Free Scalable Gaussian Processes via Randomized Truncations". In: International Conference on Machine Learning (ICML). 2021 (cit. on pp. 1, 7).
[13] Burt, C. E. Rasmussen, and M. van der Wilk. "Rates of Convergence for Sparse Vari- ational Gaussian Process Regression". In: International Conference on Machine Learning (ICML). 2019. URL: http:://araai www..orgg/bb........ org/abs/1903, 03571 (cit. on p. I).
[14] M. Kang, F. Schäfer, J. Guinness, and M. Katzfuss. Asymptotic properties of Vecchia ap- proximation for Gaussian processes. 2024. DOT: 10. 48550/ arxiv. 2401 15813 (cit. on p. 1).
11

---

[15] S. W. Ober, D. R. Burt, A. Artemev, and M. van der Wilk. "Recommendations for Baselines and Benchmarking Approximate Gaussian Processes". In: NeurIPS Workshop on Gaussian Processes, Spatiotemporal Modeling, and mecision-making Systems. 2022. URT.: https:// /// (cit. on pp. 1,7).
gP peseminar-seriee serriiie.gttuub... io/ assetsss///////s////////////e tersstssccmmraaa sateerrrrrddyyee6
[16] Z. Wang. C. Gehring, P. Kohli, and S. Jegelka. "Batched Large-scale Bayesian Optimiza- tion in High- -igh-dimensional Spaces". In: International Conference on Artificial Intelligence and URL: http: / arxiv.
Statistics (AISTATS). 2018. DOI: 10. 48550/ aar/iv..v.... 1706. 01445.
org/ hr//abs/1706. 01 01445 (cit. on pp. 2, 7).
[17] E. Turner and M. Sahani. "Two problems with variational expectation maximisation for time series models". In: Bayesian Time Series Models. Cambridge University Press, 2011,
pp. 104-124. DOI: 10. 1017/C009780011998000
006 (cit. on pp. 2. 7).
[18] M. Bauer, M. van der Wilk, and C. E. Rasmussen. "Understanding probabilistic sparse Gaussian process approximations" In: Advances in Neural Information Processing Systems (NeurIPS). Vol. 29. 2016 (cit. on pp. 2, 7)
[19] J. Wenger, G. Pleiss, M. Pförtner, P. Hennig, and J. P. Cunningham. "Posterior and Compu- tational Uncertainty in Gaussian processes". In: Advances in Neural Information Processing Systems (NeurIPS), 2022. DOI: 10. 08.558055rrr. /ar//vriiii.///..///eaa/eiaaaa.a 2205 15449 (cit. on pp. 2, 5, 7, 16, 17,
23).
[20] D. R. Burt, C. E. Rasmussen, and M. v. d. Wilk. "Convergence of Sparse Variational Inference in Gaussian Processes Regression". In: Journal of Machine Learning Research (Aug. 2020).
DOI: 10. 18550/arXivs 2008. 00323 (cit. on p. 3).
[21] N. Houlsby, F. Huszár, Z. Ghahramani, and M. Lengyel. "Bayesian active learning for classi- fication and preference learning". In: arXiv (2011). URL: https:///artti http//arxvvvoooraabaaa/ 1112.
5745 (cit. on pp. 5, 21).
[22] H. Zhu, C. K. Williams, R. Rohwer, and M. Morciniec. "Gaussian regression and optimal finite dimensional linear models". In: Neural Networks and Machine Learning. 1997 (cit. on
p. 5).
[23] C. Lanezos. An iteration method for the solution of the eigenvalue problem of linear differen- tial and integral operators. United States Government Press Office Los Angeles, CA, 1950 (cit. 5,
[24] M. R. Hestenes and E. Stiefel. "Methods of conjugate gradients for solving linear systems". In: Journal of Research of the National Bureau of Standards 49 (1952) (cit. on pp. 21).
[25] P. Hennig. M. A. Osborne, and M. Girolami. "Probabilistic numerics and uncertainty in com- putations". " In: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471.2179 (2015) (cit. on p. 7).
[26] J. Cockayne, C. Oates, T. Sullivan, and M. Girolami. "Bayesian probabilistic numerical meth- ods". In: SIAM Review 61.4 (2019), pp. 756-789 (cit. on p. 7)
[27] P. Hennig, M. A. Osborne, and H. P. Kersting. Probabilistic Numerics: Computation as Ma- chine Learning. Cambridge University Press, 2022. ISBN: 778-1-316-6814-11 DOI: 10. 0017/9781316681444 (cit. on p. 7).
(28) M. Pförtner, I. Steinwart, Hennig, and J. Wenger. Physics-Informel Gaussian Process Re- gression Generalizes Linear PDE Solvers. 2023. DOr: 10. 48550/arXiv 2212, 12474 (cit. on p. 7).
[29] L. Tatzel, J. Wenger, F. Schneider, and P. Hennig. Acceleruting Generalized Lineur Models by Trading off Computation for Uncertainty. 2024. DOI: 10. 48550 / arxiv. 2310 20285 (cit. on p. 7).
[30] M. Pförtner, J. Wenger, J. Cockayne, and P. Hennig. Computation-Aware Kalman Filtering and Smoothing. 2024. DOt: 10. 48550/ 48550/arxiv arxiv. 2405. 08971 (cit. p. 7).
[31] D. Hegde, M. Adil, and Cockayne. Calibrated Computation-Aware Gaussian Processes.
2024. DOI: 11. 4844488555Xii 2410. .87996 (cit. on p. 7).
[32] P. Hennig. "Probabilistic Interpretation of Linear Solvers". In: SIAM Journal on Optimization
25.1 (2015), pp. 234-260 (cit. on p. 7).
[33] J. Cockayne, C. Oates, I. C. Ipsen, and M. Girolami. "A Bayesian Conjugate Gradient Method", In: Bayesian Analysis 14.3 (2019), pp. 937-1012 (cit. on p.p.7 7).
[34] S. Bartels, J. Cockayne, I. C. Ipsen, and P. Hennig. "Probabilistic linear solvers: unifying view". In: Statistics and Computing 29.6 (2019), pp. 1249-1263 (cit. on p. 7).
12

---

[35] J. Wenger and P. Hennig. "Probabilistic Linear Solvers for Machine Learning". In: Advances in Neural Information Processing Systems (NeurIPS), 2020. DOI: 10. 00..44550///a0 http://XXvv....00 2010. 09691. URL: //atp///rxivv org/ abs 1200000....... 096991 (cit. on p. 7).
[36] H. Liu, Y.-S. Ong, X. Shen, and J. Cai. "When Gaussian process meets big data: A review of scalable GPs". In: Transactions on Neural Networks and Learning Systems 31.11 (2020),
Pp. 4405 4423 (cit. on p. 7).
[37] J. Hensman, A. Matthews, and Z. Ghahramani. "Scalable variational Gaussian process clas- sification". In: International Conference on Artificial Intelligence and Statistics (AISTATS). PMLR, 2015. pp. 351-360 (cit. on p. 7).
[38] H. Salimbeni, C.-A. Cheng, B. Boots, and M. Deisenroth. "Orthogonally Decoupled Vari- ational Gaussian Processes". In: Advances in Neural Information Processing Systems (NeurIPS). Vol. 31. 2018 (cit. on p. 7).
[39] L. Wu, G. Pleiss, and J. P. Cunningham. "Variational nearest neighbor Gaussian process". In: International Conference on Machine Learning (ICML). PMLR, 2022, pp. 24114-24130 (cit. on p. 7).
[40] J. Hensman, N. Durrande, and A. Solin. "Variational Fourier Features for Gaussian Pro- cesses". In: Journal of Machine Learning Research 18.151 (2018), pp. 1-52 (cit. on p. 7). [41] V. Dutordoir, N. Durrande, and J. Hensman. "Sparse Gaussian Processes with Spherical Har- monic Features". In: International Conference on Machine Learning (ICML.). Vol. 119. 2020, pp. 2793-2802 (cit. on p. 7).
[42] M. Van der Wilk, V. Dutordoir, S. John, A. Artemev, V. Adam, and J. Hensman. A framework for interdomain and multioutput Gaussian processes. 2020. DOE: 10. 48550/arXiv 2003. 01115 (cit. on p. 7)-
[43] A. Artemev, D. R. Burt, and M. van der Wilk. "Tighter Bounds on the Log Marginal Likeli- hood of Gaussian Process Regression Using Conjugate Gradients". In: International Confer- ence Machine Learning (ICML). 2021 (cit. on p. 7).
[44] M. Jankowiak and G. Pleiss. Scalable Cross Validation Losses Gaussian Process Models.
2022. DOI: 10. 48550/arKii. 2105. 11535 (cit. on p. 7).
[45] M. Jankowiak, G. Pleiss, and J. Gardner. "Parametric Gaussian Process Regressors". In: In- ternational Conference on Machine Leurning (ICML). Vol. 119. 2020, pp. 4702-4712 (cit. on
Pp. 7. 9).
[46] Y. Wei, R. Sheth, and R. Khardon. "Direct Loss Minimization for Sparse Gaussian Processes". In: International Conference on Artificial Intelligence and Statistics (AISTATS). Vol. 130.
2021, pp. 2566-2574 (cit. on p. 7).
[47] A. Tsanas and M. Little. Parkinsons Telemonitoring. UCI Machine Learning Repository. 2009. DOr: 10. 24432/ (cit. on p. 8).
[48] H. Fanace-T and J. Gama. Bike Sharing. UCI Machine Learning Repository. 2013. DOI: 10.
24432/C5W894 (cit. on p. 8).
[49] P. Rana. Physicochemical Properties of Protein Tertiary Structure. UCI Machine Learning Repository, 2013, DOI: 10. 13443/C55QW3 (cit. on p. 8).
[50] M. Naeem and S. Asghar. KEGG Metabolic Reaction Network (Undirected). UCI Machine Learning Repository. 2011. DO1: 10. 24432/C56609 (cit. p. 8).
[51] M. Kaul. 3D Road Network (North Jutland, Denmark). UCI Machine Learning Repository.
2013. DOI: 10. 24432/C5GGP1 (cit. on p. 8).
[52] G. Hebrail and A. Berard. Individual Household Electric Power Consumption. UCI Machine Learning Repository. 2006. DOI: 10. 24432/C58K54 (cit. on p. 8).
[53] M. Kelly, R. Longjohn, and K. Nottingham. The UCI Machine Learning Repository. 2017. URL: https: whh.vee...ve.e// ics uci edu (cit. on p. 7).
[54] D. P. Kingma and J. Ba. "Adam: A method for stochastic optimization". In: International
Conference on Learning Representations (ICLR) (2015) (cit. on p. 7).
[55] J. Nocedal. "Updating quasi-Newtor matrices with limited storage". In: Mathematics of Com-
putation 35.15 151 (1980). pp. 773-782 (cit. on p. 7).
13

---

[56] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury. G. Chanan, T. Killeen. Z. Lin, N. Gimelshein, L. Antiga, Desmaison, A. Köpf. Yang, Z. DeVito, M. Raison, A. Te- jani, S. Chilamkurthy, B. Steiner, I.. Fang, J. Bai, and S. Chintala. "PyTorch: An Impera- tive Style, High-Perffrmance Deep Learning Library". In: Advances Neural Information Processing Systems (NeurIPS). 2019. DOT: 10. 48550/ arXiv 1912 01703. URL: http:
/ http:///////b//////////// 01703 (cit. on p. 7).
[57] G. Pleiss, M. Jankowiak, D. Eriksson, A. Damle, and J. Gardner. "Fast matrix square roots with applications to Gaussian processes and Bayesian optimization". In: Advances in Neural Information Processing Systems (NewrIPS). 2020, pp. 22268-22281 (cit. on 9).
[58] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand. "Hierarchical nearest-neighhor Gaus- sian process models for large geostatistical datasets". In: Journal of the American Statistical Association 11 111.514 14 (2016), pp. 800-812 (cit. p. 9).
[59] M. Katzfuss and J. Guinness. "A General Framework for Vecchia Approximations of Gaus- sian Processes". In: Statistical Science 36.1 (2021). DOI: 10. 1214/19 (cit. on p- 9)-
[60] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur. Gaussian Processes and Kernel Methods: A Review on Connections and Equivalences. arXiv:1807.0258; [cs, stat]. 2018. DOI: 10. 48550/arXiv 1807 02582 (cit. on p. 16).
[61] V. Wild, M. Kanagawa, and D. Sejdinovic. "Connections and Equivalences between the Nyström Method and Sparse Variational Gaussian Processes". In: arXiv (2021), URL: http:
/ //arxiv.. org/ abb/2106......////////////////// 01121 (cit. pp. 16, 17).
[62] G. H. Golub and C. F. Van Loan. Matrix computations. John Hopkins University Press, 2012 (cit. on pp. 21, 25).
[63] G. Meurant and Z. Strakoš. "The Lanczos and conjugate gradient algorithms in finite pre- cision arithmetic".. en. In: Acta Numerica 15 (May 2006), pp. 471-542. DOE: 10. 1017 $0962492906260011 (cit. on 21).
14

---

Supplementary Material
This supplementary material contains additional results and in particular proofs for all theoretical statements. References referring to sections, equations or theorem-type environments within this document are prefixed with 'S', while references to, results from, the main paper are stated as is.
S1 Theoretical Results
15
$1.1
Altemaatiee Derivation CaGP Posterior
15
S1.2
Worst Case Error Interpretations of the Variance of Exact GPs, CaGPs and SVGPs
16
$1.3
CaGP's Variance Decreases Monotonically the Number of Iterations Increases.
17
S2 Training Losses
17
S2.1
Projected-Date Log-Marginal Likelihood
17
$2.2
Evidence Lower Bound (EIBO)
18
S2.3
Comparison of Training Losses
20
S3 Choice of Actions
20
$3.1
(Conjugate) Gradient/ Residual Policy
21
S3.2
thommation-theoretie Policy
21
S4 Algorithms
23
$4.1
Iterative and Batch Versions CaGP
23
S4.2
Implementation.
23
SS
Additional Experimental Results and Details
24
$5.1
Inducing Points Placement and Uncertainty Quantification of SVGP
..C.............................
24
S5.2
Grassman Distance Between Subspaces
25
S5.3
Generalization Experiment
25
S5.3.1 Impact of Learning Rate Generalization
25
S5.3.2
Evolution Of Hyperparameters During Training
27
S1 Theoretical Results
$1.1 Alternative Derivation of CaGP Posterior
Lemma S1 (CaGP Inference as Exact Inference Given Modified Observation Model)
Given Gaussian process prior f GP(u, K) and training data (X, the computation
GP posterior GP (ps, Ki) (see Equation (5)) with linearly independent and fixed actions equiv-
alent to exact butch GP posterior (f T) given data ÿ = S' observed according to the likelihood y f(X)~= N f(X), I where S' = S chol(ST S)-
Proof. First note that by definition 5' has orthonormal columns, since
T s' = (S cho S) T ch S)- = SS (T ((SS)-1 s>- chol(S" S)-T
L LL
Now by basic properties of Gaussian distributions, we have for arbitrary X. Ro xd that ( ((x))
( ((xx)) >>(((>(()))
( (S"
(5 K(X, S' + S' S k(X)
(xx))
(xX.X)s
K(Xo,Xo)
jointly Gaussian, where used that I= S'
Therefore we have that (f(X.) Ø) ~ N(m((X)), K, Xo)) with
(X,)) = ((X.)) + K(X,, X)S(S (K X) +02 StT p): x(Xo) + K K(X., () + (T n)
15

---

Ki(Xo, Xx) = K(X,, X.) K(X. (KX K((,) Xx)
o, X.) k(X., () x kT X.)
which is equivalent to the definition of the CaGP posterior in Equation (5). This proves the claim.
S1.2 Worst Case Error Interpretations of the Variance of Exact GPs, CaGPs and SVGPs
In order to understand the impact of approximation on the uncertainty quantification of both CaGP and SVGP, it is instructive to compare the theoretical guarantees they admit, when the latent function assumed to be in the RKHS of the kernel. In the context of model selection, this corresponds to the ideal case where the optimization has converged to the ground truth hyperparameters.
Let ʃ GP(O, K) be Gaussian process with kernel K and define the observed process v(.) f(+) (1)+oe(f ae()) where ~ GP(O, 5) is white noise with noise level o² > 0, i.e. 8(xr, x') = 1(x = x¹).
Consequently, the covariance kernel of the data-generating process y(+) is given by K (x, x')
K ((x, x²) + 020(2,20 x²) and we denote the corresponding RKHS as Hl
For exact Gaussian process inference, the pointwise (relative) worst-case squared error of the poste- rior mean is precisely given by the posterior predictive variance.
Theorem S1 (Worst Case Error Interpretation of GP Variance [60])
Given a set of training inputs Ii,..., In E X. the GP posterior GP(), K.) satisfies for any
+ Bj that
sup yE llllll go $1
(u(2))) p. (x=)))
= K, (2o, 2o) + o²
(S15)
posterior mean
prodictive variance
If of 0, then the above also holds for I.- = z.
Proof. See Proposition 3.8 of Kanagawa et al. [60].
CaGP admits precisely the same guarantee just with the approximate posterior mean and covariance function. The fact that the impact of the approximation on the posterior mean is exactly captured by the approximate predictive variance function is what is meant by the method being computation- aware.
Theorem S2 (Worst Case Error Interpretation of CaGP Variance [19])
Given a set of training inputs B1,..., Fn X, the CaGP posterior (p(ps Ki) satisfies for any x + 2j that
sup MEHge lllllll llo 2
(u(2.) n (x)))²
Kz xo) +
(S16)
appouximate posterior
approximate prodictive variance
02 = 0, then the above also holds for I. Zz.
Proof. See Theorem 2 of Wenger et al. [19].
While SVGP also admits decomposition of its approximate predictive variance into two (relative) worst-case errors, neither of these is the error we care about, namely the difference between the date-genng function y E H and the approximate posterior mean uisvcp(xt)) (x.). only involves a worst-case eтror term over the unit ball in the RKHS of the approximate kernel Q" = K".
Theorem S3 (Worst Case Error Interpretation of SVGP Variance [61])
Given a set of training inputs æ1,..., Fn E X and (fixed) inducing points Z E <Rmxd, the optimal variational posterior GP (J'vvaP: Kšvap) of SVGP is given by
USvGP ( = K(Eo, () Z) K(Z, 2)) (X) Kšvop (Xo, x²)) K(xo, 2.) Q( x')) + (xx, Z)(K Z) o 2 2K X)K(XX -))
where Q(x,, x²) = K(r,, 2)K (Z, Z, x²) the Nyström approximation of the covariance function K K(x x.) (see Egns. (25) and (26) of Wild et al. [61]). The optimized SVGP posterior
16

---

satisfies for any # 2 that
(u(2))
HSver (x:)) )3
+.
sup ME DEH a
Sl
ertor of exact posterior mean assuming v(-) is the RKHS of the approximate kernel Q"
(f(z)) K (20, Z Z) SZ))
sup /eHx, [ $1
(S17)
error of exact posterior mean given noise-frec observations inducing points
K5vGP (o, 2o) + o² approximate predictive variance
=
o 0, then the above also holds for Io Z-.
Proof. See Theorem of Wild et al. [61].
S1.3 CaGP's Variance Decreases Monotonically as the Number of Iterations Increases
Proposition S1 (CaGP's Variance Decreases Monotonically with the Number of Iterations) Given a training dataset of size n, let GP(M, K,) be the corresponding CaGP posterior defined in Equation (5) where i S n. denotes the downdate rank / number of iterations and assume the CaGP actions S [x are linearly independent. Then it holds for arbitrary ws E and i S j S n, that
K2 xo) 2 (z( Ro) 2 K2, (I,, (x xc)
(S18)
where (x,,,, 2o) is the variance of the exact GP posterior in Equation (1).
al. [19] originally defined the approximate precision matrix C, Ei=1 - dpd! Ei=1 d.d. as of rank-1 matrices and show that this definition is equiv- Te alent to the batch form C; = (S)) S.. we use in this work [see Lemma S1, Eqn. (S37) in 19]. Therefore we have that
Proof. Wenger
Ki(wo, 2o)
K (oo x.)
K( X)C,K(X, v.)
K (wo, vo) Σ K(wo; X)did] K ((X,..)
= =K - Σ
(K(xo X)di}
t=1
K(Fo, v.)
Σ
K (No, X)dj}
since i S
K K(wo, Øo)
Σ (K(wo; X)dr)²
K (o)
K (xo, X)C,
K. (. ato)
where the last equality follows from the fact that S, E nuxa has rank and therefore
C (S, -SST S k - (S,S,") K-
S2
Training Losses
S2.1
Projected-Daad -ata Log-Marginal Likelihood
Lemma S2 (Projected-Date Log-Margina Likelihood)
Under the assumptions of Lemma SI, the projected-data log-marginal likelihood given by
(NT (0) - -logplp 0) - log (j; S. p. S. ks:)
17

---

( (S)) (v 14) + log d K.S)) log dct(S) S,) + ilog(2x))
Proof.
By the same argument as in Lemma S1 obtain that
(u:j(jjs ( n. S k (proj (0) log 0) -og. N
1 ((5 y s u (S k KS)) (s) S 14) + log KSS)) + i'log(2m)
2
(vv - 1²) ( k 5 (y 1 + dS KS))) + (log((2))
Since = T where L T-T chol (5 S,) is the orthonormalizing matrix, cancels in the quadratic loss term, giving
1 2 (vv - Ts,
μi) + log det(s; KS)) + :lo((2x))
and finally we can decompose the log-determinant into a difference of log-determinants
1 2 (Iy s, s!! (y- p1) + dot KS,) log det [) + (2g(2)))
2 (y - s, KS.) s! ss (v (4) + log de! K.S)) log det [LT) log (2m))
which using LL = S, completes the proof.
S2.2 Evidence Lower Bound (ELBO)
Lemma S3 (Evidence Lower Bound Training Loss)
Define the variational family
Q = = (N::( K,(X, X)) S E Rnxi]
(S19)
then the evidence lower bound (ELBO) is given by
CCaGP (0)
log ply 0) (L(((f) p(f v))
= E.y p(y f)) + (L()()) p(f))
1
r² lly m4 (x)1
K(( (n ())o))og((3) nolog(27)
==1
+ ST "TSSKSö KSö; - tr(( KS) S KS) + log det( (KS) log det(sT 5)
where o; = -ST ST (1) are the "projected" representer weights.
Proof.
The ELBO is given by
ELBO CaGP
(0) = E. (logp(y) f)) - KI(((f)) p(f)).
We first compute the expected log-lik regliikllhoodd term.
1
Eg (ly f)) EE,
_
(y-
f) f) log ( + n log(2T)
2
1
1
Eg
((v
f) -)) + log(o +
2
02
Now using B Ax)
T(x)
T
A AE(æ)
tr(A (o)()))),, obtain
1
-
1 n.
12 ((( + log((? + n log(2T)
(x)
G2
18

---

2 の (1v ps (X)
1
Σ (i))),
log(o?) + n log (2+)
Since both qff) and the prior p(f) are Gaussian, the KL divergence term between them given by
KL(a(f) |I p(f) 1 -
-
x(x)) ((x))) (_(())(X))))))))))(((((())()))<) ((x(x)-((x)) + log (Xe((K:(X,>)) x))
det(K)
+ Ta(KK'',(X,,)) -n)
(KCi(y (KC))) KC:ly x(x)) log (X
+ +tr(Inxn--Ca - -CKK -n)
1
yy (XC-)(())CCCC ((X)) log det -
CC,K) +t(K)ttr(Ixxx))) tr -CK) -
2
1
S. KSū; logdet((nxx. C₂K) tr(Inxx CCK)- m) -n)
2
logdet(Inxx C₂K) tr( -S KS)
0
Next, we use the matrix determinant lemma det (A + AVV = det Im + V A A-
0. S.S. KSo: - log - 15 -
2
KSü; - logdet(( (S)) KS 5 KS)) (()
ST KS)
2
J.s.K KSU; logdet (S)) 8)) - -
2
8
S. KSU; - det ST KS) log det 0² ST "
- Ks)) - S S K5)
2
STSKSS - leg KS) iloo(o - log det ST 5) tr((St KS) ST KS) (TSK))
2
19

---

S2.3
Comparison of Training Losses
2
Train Loss
Test NLL
0
Test MSE
10-1
10-1
Epoch
Wallclock Time (s)
CaGP-CG Projected-Dana Data Log- Marginal Likelihood
CaGP-CG ELBO Training Loss
Figure S1: Comparison of two different training losses for CaGP. The naive choice of the projected- data log-marginal likelihood leads to increasingly worse generalization performance as measured by NLL. In comparison, the ELBO training loss leads much better performance.
CaGP-CG + Projected- Datt Log-Marginal Likelihood
CaGP-CG ELBO Training Loss
0--
-5-
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Test data
GP mean
GP uncertainty
Figure S2: CaGP predictive distributions with hyperparameters optimized using different losses. When optimizing hyperparameters with respect to the projected-data log- magginala likelihood, CaGP-CG completely overestimates the noise scale, which leads to increasingly worse general- ization performance. comparison, the ELBO training loss leads a much better overall fit.
S3 Choice of Actions
We begin by proving that the CaGP posterior in Equation (5) is uniquely defined by the space spanned by the columns of the actions colsp(S), rather than the specific choice of the matrix S.
Lemma S4 (The CaGP Posterior Is Uniquely Defined by the Column Space of the Actions) Let S, S' Raxi be two action matrices, each of which consists of non-zero and linearly indepen- dent action vectors, such that their column spaces identical, i.e.
colsp( S) colsp(S'),
(S20)
then the corresponding CaGP posteriors GP(4; K:) and GP(p,, K!) are equivalent.
Proof. By assumption (S20) there exists W iixi such that S' SW. Since action vectors are assumed to be linearly independent and non-zero, it holds that rank(5") = rank(SW) rank W), where the last equality follows from standard properties of the matrix rank. Therefore W is invertible, and we have that
C' = S KS - -
SW (WWWTSS (WTSTKSW) - ST S(S Kis C.
20

---

Since the CaGP posterior in Equation (5) is fully defined via the approximate precision matrix C, the desired result follows.
Corollary S1 (Action Order and Magnitude Does Not Change CaGP Posterior) The CaGP posterior in Equation (5) is invariant under permutation and rescaling of the actions.
Proof. This follows immediately by choosing permutation diagonal matrix W, respectively, such that S' SW in Lemma S4.
S3.1 (Conjugate) Gradient/ Residual Policy
Consider the following linear system
Kv. y μ
(S21)
with symmetric positive definite kernel matrix K K I observations y. prior mean evaluated at the data n ()) and representer weights vx..
Lanczos process [23] The Lanczos process is an iterative method, which computes approximate eigenvalues Å diaagAi...... 2s) E Rixi and approximate eigenvectors Ü = (ũ ....;,.))))) ú;) R
for a symmetric positive definite matrix K by repeated matrix-vector multiplication. Given an arbitrary starting vector 91 R², s.t. | | 1, it returns i orthonormal vectors Q = (qi q:)
Raxi and a tridiagonal matrix T Q KQ E Rixi, The eigenvalue approximations are given by an eigendecomposition of = AW where W E gixi orthonormal, and the eigenvector approximations are then given by Ú = QW E maxi [e.g. Sec. 10.1.4 of 62].
Conjugate Gradient Method [24] The conjugate gradient method is iterative method to solve linear systems with symmetric positive definite system matrix by repeated matrix-vector multiplica- tion. When applied to Equation (S21), it produces sequence of representer weights approximations
V; U, K (y- μ.). Its residuals
y μ K are proportional to the Lanczos vectors
ro i.e. Q = RD where E Rxxt is the matrix of for a Lanczos process initialized at = |rall"
residuals and D Rixi a diagonal matrix (e.g. [Alg. 11.3.2 in 62] or [See. Equ. (3.4) of 63]).
Therefore choosing actions defined by the residuals of CG CaGP-CG, i.c. S = R, is equivalent to choosing actions S' = Ú given by the eigenvector approximations computed by the Lanczos process initialized as above, since
colsp( (S) colsp(R) colsp RD) CQI colsp(QW) colsp(Ú) colsp(S')
and by Lemma S4 it holds that the corresponding CaGP posteriors with actions S and S' are equiv- alent.
S3.2 Information-theoreti Policy
In conommation-thhoredee formulations of active learning, new data is selected to minimize uncertainty about set of latent variables In other words, would aim to minimize the entropy of the
posterior (z) = log X)p(z X) dz as a function of the data X [21]. In analogy to active learning. in our setting we propose to perform computations y + S. to maximally reduce uncertainty about the latent function fXX) evaluated at the training data.
Lemma S5 (Information- coaematton--thoretiic Policy)
The actions S minimizing the entropy of the computation- conputation-awwar aware posterior (x) STy) at the training data, or equivalently the actions maximizing the mutual information between f(X) and the projected data S y, are given by
(S1,..., 8;) arg min ( (f()))
(S22)
SERext
arg max H(f(X)) - ((f STyy)
(S23)
SEROXi
=MI(f(X)):TTT)
where S1...., Sü are the top-i eigenvectors of K in descending order of the eigenvalue magnitude.
21

---

Proof. Let ÿ == S y and f = f(X). By assumption, we have that f ~ (u, K). Recall that
the entropy of a Gaussian random vector f ~ (m,,)) is given by H(f) = 1 log det + (og()) Now since the covariance function of the computation-aware posterior in Equation (5) does not depend on the targets neither does its entropy Hp(f)sT (f).
Therefore, by definition of the conditional entropy and using the law of the unconscious statistician, it holds that
H(f 9) = 1000--] l log log 1pp(f of
Epiin (Hrifls((f)
pl) (HInF)sTT))())
and since the covariance of a Gaussian conditioned on data doesn't depend on the data, we have
H p(nSst)
Therefore we can rewrite the mutual information in terms of prior and posterior entropy, such that H(f) Hf y) H(f) H))
log det (K) + n mlog(2re) - det(de KCCKK - log(22e)
1 log (det(K- - KS(S - det((K)))
Via the matrix determinant lemma dct(AAUWV)) det W - det (W) det(A), we obtain
1
log (det(- - dKk--K)de(((SST KS) det
det K 1
2
1
= log (det(S KS ST ds(
2
= log (ots
2
log det (S
log det ((5 log(o²))
log det LL T ilog(o2)
for L a square root of S Now can upper bound the above as follows
H(f STy) y)
(log (log det( ST K5) :log(o²))
max SER#x+
max SERax: 2
log det det(UU - -
1
2 log(A,(K)) - illo(o?))
where U are the orthonormal eigenvectors of K for the largest eigenvalues. Now choosing $ U achieves the upper bound since T = and therefore S, = U is a solution to the optimization problem.
Finally using the argument above and since H(f) does not depend on S, we have that
arg max H(f) - H(f (Ty) arg max H(f) H H) (f) = arg min H. Hp( (f). SeRex1
SERaxi
SeRaxi
This proves the claim.
22

---

S4
Algorithms
S4.1 Iterative and Batch Versions of CaGP
Algorithm S1: CaGP IterGP: Iterative formulation as in Wenger al. [19]
Input: GP prior (P, K). training data (X, y) Output: (comhined) GP posterior ()( Kj)
procedure CaGP(u, K, X, y, C. 0)
Time
Space
2 while not STOPPINGCRITERIONO do
3
8; <- POLICY()
Select action via policy.
4
< (y 1t) Ku;-1
Residual.
O(n²)
(m)
5
0: - s/r-11
Observation.
O(k)
O(1)
6
t Ks,
O(nk)
O(n)
7
d + Σ,-1 Ks, = 8: Ci-1 d:
Search direction.
O
O(n)
7. -- 8 KE Ka. = = =,dd
((n)
0(1)
9
C, < C:-1 ddd
Precision matrix approx. C₁ K
O(n)
O(ni
10
v. < 8:11 + 21 d:
Representer weights estimate.
O(n)
O(n)
11
E, - Eg C₁
Representer weights uncertainty.
12
na(-) p()) K)
Approximate posterior mean.
((non)
O(no)
13
K.(:.)) + k(,,)) -) K(; x)C,K(X
Combined covariance function.
(n,m)))
O
return (P(Mi,
Algorithm $2: CaGP: Batch Version
Input: GP prior (p, K), training data (X,y) Output: (combined) GP posterior GP(p,,
procedure CaGP(μ, K, x,
Time
Space
2
S, POLICY()
Select batch actions via policy. "Projected" data.
3
ÿ - s (y 14)
O
O
4
2, + KS.
O
O
5
L, 4 (SOLESKY(S Z,)
O(P(1+ R))
(
8; < L. L,
"Projected" representer weights.
((P)
O
7
Ks(:X)) ((6,))(Ki [+(K(x)),
O
O(n.t)
8
m²) + -1()) (Ks), X)ú,
O(noi)
O(n.)
9
Ki(,,) + K(.)) Ks(,X)LL X) kss(X,))
O(nof)
O
10
return GP(p;,
S4.2
Implementation
We provide an open-source implementation of CaGP-Opt as part of GPyTorch. To install the pack-
age via pip, execute the following in the command line:
pip install http:///ttteeegee........aaaaaaa http://./.ooreellleee///eeeeet www.riatrrrpprroorooooo ty pip install stt-https://githhubb http:///wooooooo/////oooooo/o/ee www.ccomtttaiiooooaaaarrrrroreea pip install pykeops
23

---

S5 Additional Experimental Results and Details
S5.1 Inducing Points Placement and Uncertainty Quantification SVGP
To better understand whether the overconfidence of SVGP at inducing points observed in the visu- alization Figure 1 holds also higher dimensions, we do the following experiment. For varying input dimension d E {1, (1.2..... 25). we generate synthetic training data by sampling inputs X uniformly at random with corresponding targets sampled from zero-mean Gaussian pro- cess ~ GP(O, K where (- ) = (-...) + o 5(-,-) is given by the sum of Matérn(3/2) and a white noise kernel with noise scale a. We optimize the kernel hyperparameters, variational
500
parameters and inducing points (m 64) jointly for 300 epochs using Adam with a linear learning rate scheduler. convergence we measure the average distance between inducing points and the nearest datapoint measured in lengthscale units, i.e.
dZz
minl|z; 23 [lliing(1=2)
(S24)
m
where E Rd is the vector of lengthscales (one per input dimension). We also compute the average ratio of the posterior variance to the predictive variance at the inducing points, i.e.
172
K pasterior (Zi, 2.)
(
(S25)
771
K Xposterior Zis 2i) + 02
The results of our experiments are shown in Figure S3. We find that expected the inducing points are optimized to closer to datapoints than points sampled uniformly at random. However, inducing points lie increasingly far away from the training data the dimension increases relative to the lengthscale that SVGP learns. Therefore this experiment suggests that the phenomenon observed in Figure 1, that SVGP can be overconfident at inducing points if they are far away from training datapoints, to be increasingly present as the input dimension increases. This is further substantiated by Figure S3(b) since the proportion of posterior variance to predictive variance at the inducing points is very small already in d = dimensions. This illustrates both SVGP's overconfidence at the inducing points (in particular in higher dimensions) and that its predictive variance is dominated by the learned observation noise, as we also saw in the illustrative Figure 1.
Distance to data di Z X
1.00
0.6
Random Inducing Points
Variance ratio (
SVGP's Optimized Inducing Points
0.75
CholeskyGP w/ true hyperparameter:
0.50
SVGP
0.2
0.25
0.0
0.00
12
16
20
24
12
16
20
24
Input Dimension
Input Dimension
(a) Average distance inducing points the nearest datapoint measured in lengthscale units.
(b) Average ratio posterior to predictive variance at SVGP's inducing point locations.
Figure S3: SVGP's inducing point placement and uncertainty in higher dimensions. (a) As the dimension increases, the inducing points SVGP learns lie increasingly far away from the data mea- sured in lengthscale units given a fixed training data size and number of inducing points. (b) SVGP's variance at the inducing points is dominated by the learned observational noise in higher dimensions, rather than by the posterior variance. The comparison to a CholeskyGP with the data- generating hyperparameters shows that SVGP compensates for lack of posterior variance at the in- ducing points by artificially inflating the observation noise. This illustrates both the overconfidence (in terms of posterior variance) of SVGP at the inducing points and its tendency to oversmooth.

---

S5.2 Grassman Distance Between Subspaces
In Figure 2 we compute the distance between the subspaces spanned by random vectors, the actions S of CaGP, and the space spanned by the top-i eigenvectors. The notion of subspace distance we use the Grassman distance, i.e. for subspaces spanned by the columns matrices A E gaxp and B E TR s.t. 2 q the Grassman subspace distance is defined by
A, B) = ||0|2
(S26)
where 0 E IR9 the vector of principal angles between the two spaces, which can computed via an SVD [e.g. Alg. 6.4.3 in 62].
S5.3
Generalization Experiment
Table S1: Detailed configuration of the generalization experiment in Section 5,
Posterior Approximation
Model Selection Training
Iters, i Ind. Points
Solver Tol.
Optimizer
Epochs
(Initial) Learning Rate
Batch Size
Procision
CholeskyGP
LBFGS
100
[1,10"
10 10"3 10
Floate4
SGPR
1024
Adam
1000
[1, .100
10
float32
LBFGS
..10
2.10
SVGP
1024
Adam
1000
11,10-1,0-2,4
1,100
1024
float32
CGGP
512
10
LBEGS
100
(1, 1.100".1072,10
float64
CaGP-CG
10-4
(1,10-1)
CaGP-Opt
512
Adam
1000
(1,10 10
float32
512
L.BFGS
100
(1, 10 10
float64
S5.3.1 Impact of Learning Rate on Generalization
To show the impact of different choices of learning rate on the GP approximations we consider, we show the test metries for the learning rate sweeps our main experiment in Figure S4. Note that not all choices of learning rate appear since small minority of runs fail outright, for example if the learning rate large.
25

---

Train Loss
Test NLI
Test RMSE
10"
Parkinsons
10-1
-
10-2
33
10"
3min
10a
10-1
10-2
4mii
Smin
4min
Smin
4min
8min
2.0
1.00
0.7
Protein
1.5
0.85
0.80
0.5
8min
16min
8min
16min
Smin
16min
0.5
0.5
2
0.4
KEGCCu
0.3
0.5
0.2
20min
20min
2
1.0
1.00
0.5
0.50
0.0
0.25
0
1h
2h
1h
2h
1h
2h
2
10a
0-
10-1_
0
1.5h
3/
4.5h
1.5h
4,5h
1.5b
3b
4.5h
Time
Time
Time
model
learning rate
CholeskyGiP
0.0001
SGPR
0.001
SVGP
0.01
COGP
CaGP-CG
CaGP-Opt
Figure $4: Effects of (initial) learning rate when using either LBFGS with Wolfe line search (CholeskyGP. SGPR) or Adam (SVGP, CaGP-CG, CaGP-Opt) for hyperparameter optimization.
26

---

S5.3.2 Evolution Of Hyperparameters During Training
To better understand how the kernel hyperparameters of each method evolve during training, show their trajectories in Figure S5 for each dataset. Note that we only show the first three length- scales per dataset (rather than up 26).
Outputscale
Observation Noise
Lengthscale;
Lengthscale: d=2
Lengthscale: d=3
Parkinsons
10
10
10a
100
10
10-1
10-3
10-1
100
1.7min
1.7min
1.7min
1.7miin
1.7min
10
10a
109
Bike
10-1
10-3
10
Amin Smin
4min 8min
Amin Smin
4min Smin
4min 8min
2 x 10"
4 x 10"
3 10"
Protein
6 10
10a
10
2 x 100
x 10
10-1
x 10-1
109
108
0
Erinn 16min
Smin 16min
Emin 16min
Smin 16min
8min 16min
10a
109
KBOGG
100
10
10
104
10-2
10
109
108
20min 40min
20min 40min
20min 40min
0
20min 40min
20min 40min
10a
10a
10
10-1
10-1
10-1
1.4h
10a
Power
10
10-1
10.
104
10-2
10:
10-3
109
10a
109
2.8h
2,8h
Time
Time
Time
CholeskyGP
SGPR
SVGP
CGGP
CaGP-CG
CaGP-Opt
Figure S5: Learned hyperparameters for different GP approximations UCI datasets. Showing only results for the best choice of learning rate per method.
27

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We give both illustrative and theoretical justification for the claims about our method in Section 3 and provide extensive empirical results in Section 5.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
We explicitly discuss the limitations of our method in the conclusion in dedicated paragraph.
Justification:
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the ap- proach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech- speech-to-text system might not be used reliably provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
28

---

Answer: [Yes]
Justification: We provide proofs to all theoretical statements in the supplementary material cite the appropriate reference for the result.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theo- rems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data provided not)?
Answer: [Yes]
Justification: We describe the experiments including datasets, number of repeats, hyperpa- rameters, method implementations and hardware in Section
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided not.
If the contribution is a dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate the research performed.
While NeurlPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature the contribution. For example
(a) the contribution primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results way produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case au- thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
29

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification: All data we use is publicly available in the UCI repository. We provide an open-source implementation of our method in Section S4.2.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https: ptttps://nips.ooo cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: /nips. cc/ wtc/public/aae popliccgguide suuidee/ooeeeee
ssi onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification: We give a complete description the choices we made for our benchmark experiments in Section 5.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
We repeat all experiments multiple times and report bootstrapped confidence intervals.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
30

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (c.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state The authors should prefer- ably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative emror rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In Section 5, we describe the specific GPUs use for each experiment and report wallclock time of all training runs.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't t make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: htttp::///neripsee wwwpehewwws...uuppppeeeeeececeee cttppccic/cc/ccccccccccccccccccc httpppppllctticieeeeeeeeep Eet//Ehhiiiiiiiii csGuidelinee"
Answer: [Yes]
Justification: Due to the methodological nature of this work, there potential harmful consequences of this work that we think need to be explicitly highlighted here.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work is foundational methodological research and does not have a spe- cific negative societal impact that we feel must be explicitly highlighted here.
Guidelines:
31

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negatíve societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This paper does not release data or models with high risk for misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by re- quiring that users adhere to usage guidelines or restrictions access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Ys]]
Justification:
We cite the appropriate source for all assets used our work. All datasets for our experiments are licensed under a Creative Commons Attribution 4.0 International (CC
BY 4.0) license.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
32

---

For scraped data from a particular source (e.g., website), the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms of use the pack- age should be provided. For popular datasets, paperswi snperswithcode www.coo//aa.... com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documenta- tion provided alongside the assets?
Answer: [NA]
Justification:
There are no new assets released with the paper.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
For crowdsourcing experiments and research with human subjects, does the pa- per include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
Justification:
This paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, cura- tion, or other labor should paid at least the minimum wage the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
33

---

Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
34


Reviews:
[{"review_id": "TszHKFL3NJ", "reviewer": "Reviewer_7YnY", "summary": "This paper presents a substantial extension to an existing class of models called IterGP. By introducing a novel training loss which combines both the hyperparameters and a sparse action matrix, the IterGP-Opt model offers linear-time scalability. Experiments over a range of UCI datasets demonstrate that IterGP-Opt matches or outperforms the benchmark SGPR and SVGP models.", "questions": "The final sentence of the abstract claims \"As a result of this work, Gaussian processes can be trained on large-scale datasets without compromising their ability to quantify uncertainty\" - Doesn't this imply we are conducting exact inference on these large-scale datasets? I'd recommend the authors consider a suitable rephrasing such as \"significantly compromising\".\n\n\"SGPR typically requires more gradient iterations during training as it introduces new parameters (the inducing point locations Z).\"\nI was surprised by this comment, because in my experience SGPR typically requires significantly fewer training iterations than SVGP. In part because training is unbatched, allowing it to exploit standard LBFGS optimisation as opposed to using Adam. And also because it is often not desirable to train the inducing point locations jointly with the hyperparameters.", "limitations": "The authors include a well written section on the limitations of IterGP-Opt.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and clearly structured.\n\nThe proposed IterGP-Opt model displays a high degree of novelty.\n\nThe results and methodology are for the most part well presented.\n\nScalable inference is an area of high impact, so potential improvements over the SVGP model are of high interest.", "weaknesses": "I have a couple of concerns regarding the models used in the experiments:\n\nThe paper cites Ober et al \"Recommendations for Baselines and Benchmarking Approximate Gaussian Processes\" as inspiration for including SGPR as a strong benchmark, but does not follow the guidance for SGPR training outlined in that paper. Both in terms of the choice of optimiser, and the approach of not jointly optimising the inducing points with the kernel hyperparameters. This is concerning as it could be having a significantly detrimental impact on the performance of the SGPR baseline - as is further evidenced by it underperforming SVGP on three of the four tasks. In any case, given the popularity of the datasets used, it ought to be straightforward to verify the SGPR performance is consistent with other publications in the literature.\n\nOne of the other methods presented in Table 1 is \"IterGP-CG\", but its status should be clarified for the reader. The first reference to IterGP-CG in the paper is in the caption to Figure 1, where it is introduced as \"IterGP-CG (ours), ...\", implying this model is a novelty of the paper. Yet later in the Background section, and also in the footnote on page 7, it is pointed out that the IterGP-CG model is not new. Given that IterGP-CG is not mentioned at all in the Conclusion section, it seems likely the \"ours\" was just a typo, but this is of course a very important point to clarify! Table 1 itself should also be updated to clarify which method(s) are claimed to be novel to the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Summary\n\nThank you for your time and effort in reviewing our paper and in particular for suggesting improvements to our benchmark experiments.\nWe have addressed your main concerns with a set of new experiments where we train SGPR using LBFGS in double precision. These changes close the gap between SGPR and SVGP, but IterGP-Opt still matches or outperforms SGPR on all datasets, except \"Bike\" (see the rebuttal PDF).\nYou can find our responses to your questions below. If anything remains unclear, we are happy to respond during the discussion period.\nShould the changes we've made to our paper alleviate the concerns you voiced, we would be grateful if you would update your score to reflect this. Thank you!\n\n### Weaknesses\n\n> The paper cites Ober et al [...] as inspiration for including SGPR as a strong benchmark, but does not follow the guidance for SGPR training outlined in that paper. [...] This is concerning as it could be having a significantly detrimental impact on the performance of the SGPR baseline - as is further evidenced by it underperforming SVGP on three of the four tasks. In any case, given the popularity of the datasets used, it ought to be straightforward to verify the SGPR performance is consistent with [...] the literature.\n\nTo ensure the SGPR baseline is well-tuned and to alleviate your concerns, we made the following improvements to our main experiment: \n\n- We trained SGPR on all datasets in double precision using L-BFGS with a Wolfe line search for 100 epochs in line with the recommendations of Ober et al. (2022). We performed a sweep over the initial learning rate $\\text{lr}_{\\text{SGPR}} \\in \\\\{10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}, 1\\\\}$ and repeated each run for five random seeds. Both in our original experiments using single precision and Adam and the new ones, we computed Cholesky decompositions with adaptively increasing kernel jitter as per Section 4.1 of Ober et al (2022). \nWe find that this indeed increases its performance. SGPR now outperforms SVGP on all small to medium datasets as expected, but IterGP-Opt still matches or outperforms SGPR on all datasets, except \"Bike\", as both Table R1 and Figure R1 in the attached rebuttal PDF show. Thank you for suggesting this improvement to the SGPR baseline!\n\n- Additionally, we checked our newly generated results against the corresponding results reported by Ober et al. (2022) in Tables 2 and 3 of Appendix H and found the newly reported performance of SGPR in our experiments to be largely consistent on all datasets in common between their and our work (\"Bike\", \"Protein\", \"KEGGu\") for $M=1000$ inducing points. However, one should be careful when comparing these numbers, since Ober et al. (2022) use a squared exponential kernel, while we use a Matern(3/2) kernel.\n\n- Finally, based on the suggestions of reviewer wkY4, we added an exact (Cholesky)GP baseline on the two smallest datasets and\nalso trained IterGP-Opt using L-BFGS. We found that this marginally improves the performance of IterGP-Opt. However, due to the higher memory consumption and slightly larger runtime, we would still recommend using Adam in practice.\n\n> One of the other methods presented in Table 1 is \"IterGP-CG\", but its status should be clarified for the reader. The first reference to IterGP-CG in the paper is in the caption to Figure 1, where it is introduced as \"IterGP-CG (ours), ...\", implying this model is a novelty of the paper. Yet later in the Background section, and also in the footnote on page 7, it is pointed out that the IterGP-CG model is not new. Given that IterGP-CG is not mentioned at all in the Conclusion section, it seems likely the \"ours\" was just a typo, but this is of course a very important point to clarify! Table 1 itself should also be updated to clarify which method(s) are claimed to be novel to the paper.\n\nWe apologize for any confusion this typo may have caused. \nAs we write in the introduction, IterGP as a class of GP *posterior approximation* methods was introduced by Wenger et al. (2022), one special case being IterGP-CG. \nWe did not intend in any way to claim credit for this method.\nThe work by Wenger et al. (2022), however, does *not* describe a way to perform *model selection*, which is one of the two main contributions of our work. In particular, we demonstrate in this work how to perform model selection for IterGP-CG. Our second contribution is that we introduce a new variant of IterGP (i.e. IterGP-Opt), which enables posterior approximation in *linear time*. We have modified the draft to more clearly indicate our contributions without inadvertently suggesting we proposed IterGP-CG.\n\n\n### Questions\n\n> [...] I'd recommend the authors consider a suitable rephrasing such as \"significantly compromising\" [in the abstract].\n\nThank you for this feedback. We will adjust the last sentence of the abstract accordingly.\n\n> SGPR typically requires more gradient iterations during training as it introduces new parameters (the inducing point locations $Z$).\" I was surprised by this comment, because in my experience SGPR typically requires significantly fewer training iterations than SVGP. [...]\n\nWe would like to politely point out that this is a misunderstanding of what we write. *We are comparing SGPR to exact GPs, not to SVGP* in the quoted sentence. We write -- *after* introducing exact GPs and *before* introducing SVGP -- in the paragraph on SGPR in lines 93-95 that \"While these complexities significantly improve upon $O(n^3)$ computation/$O(n^2)$ memory, SGPR typically requires more gradient iterations during training as it introduces new parameters (the inducing point locations $Z$).\"  We will clarify this in the final version and also add that the inducing point locations can be\nchosen in a way that does not require optimization, e.g. as outlined in Burt et al. (2020) and recommended by Ober et al. (2022).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "WvI5ygRTT0", "reviewer": "Reviewer_LXxc", "summary": "The authors cast the IterGP method of Wenger et al, which is guaranteed to give conservative estimates of uncertainty relative to the posterior, as a variational procedure. This allows for hyperparameter selection and selection of parameters controlling the quality of the approximation ($\\mathbf{S}$) through maximization of an evidence lower bound. The authors compare this method to existing variational approximations on several benchmark datasets.", "questions": "- The authors claim SVGP has “no latent uncertainty” where data is absent. Do the authors mean that the posterior variance is actually 0 or extremely close to it away from the data? I don’t see how this can be true since SVGP recovers the posterior distribution in certain limits. If what is meant is that uncertainty is underestimated away from the data, this is quite a bit weaker as a statement and should be made clear.\n- “Unlike these other methods, the IterGP posterior updates are defined by a linear combination of kernel functions on training data, which guarantees that posterior variance estimates always over estimate the exact GP uncertainty” — Could you clarify this sentence? First, what is meant by “posterior updates”? Posterior mean updates? Posterior variance updates? Posterior sample path updates? Also, I don’t think this is the property that guarantees that posterior variances are overestimated. It is easy to write down a form of posterior variance that is both a linear combination of kernel function and underestimates the posterior variance, for example 0. \n- How does the expression for the approximate predictive standard deviation (section 2) in terms of a worst case quantity over an RKHS ball compare to the analogous statement for low-rank variational methods (given in Theorem 6 of Wild et al 2023, https://arxiv.org/pdf/2106.01121)? I think a brief comparison (possibly in the appendix if space does not allow) would be useful.\n- When stating what is meant by computation aware, should there be a convergence condition? I.e. not only should the precision increase, but it should converge to the posterior precision as the amount of compute increases under some sensible restrictions.\n- I don’t understand what is plotted on the right of figure 2, or what my takeaway from it should be. What is meant by plotting the magnitude of the eigenvalues of the matrix for each $x_j$? On the plot on the right, what notion of distance between subspaces is used?\n- Should I interpret equation 12 as stating the $\\mathbf{S}$ is block diagonal? If so, is this said explicitly somewhere?\n- Could you explain in more detail why SGPR diverges in these experiments? Have you checked if this is resolved by double precision?\n\n\n### Other\n\n- Line 49, missing space after comma before “a”\n- Line 61, the authors claim that IterGP has “calibrated” confidence. Is the claim that the confidence intervals are calibrated, or that they are conservative in the sense of overestimating the posterior credible intervals?\n- Line 73, h is introduced here as the latent function and never appears again. I’m a bit unclear why it is needed.\n- Line 82, before discussing computation in model selection, it should be stated that you assume maximum (marginal) likelihood will be used. There are alternative approaches (e.g. sampling)", "limitations": "The authors point out limitations well; most importantly, the restriction to Gaussian likelihoods in the present work.", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "- The exposition is clear and well motivated.\n- Placing IterGP within a variational context both allows for a conceptual connection to SGPR/SVGP and (more importantly) practical progress in terms of selection of parameters.", "weaknesses": "- Datasets have not been properly cited. In particular, a general citation to UCI doesn’t give credit to the creators of each dataset. License information was also not stated.\n- I found the discussion of information gain not well-motivated. The definition of information gain was assumed to be known, and it wasn’t clear from the text why maximizing this quantity would lead to a good quality approximation. If the only goal is to suggest projection onto the largest $i$-eigenvalues, there are many results in the matrix approximation literature that also motivate this approach. Moreover, there is prior work in the Nystrom approximation literature considering the top $i$-eigenvalues as a subspace for performing inference, for example Ferrari-Trecate et al, https://proceedings.neurips.cc/paper/1998/hash/55c567fd4395ecef6d936cf77b8d5b2b-Abstract.html.\n- The phrasing in the caption of figure 4 conflates credible and confidence intervals. The accuracy of the credible interval should be with reference to the (exact) posterior and cannot be shown by frequentist coverage. To be clear, I think checking the empirical coverage of the methods is reasonable, but don’t think you can draw conclusions about the accuracy of the credible interval from this.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Summary\n\nThank you for your positive review of our paper and helpful input! \nWe responded to your questions below. In particular, we've reformulated the motivation behind choosing actions according to an information-theoretic objective and we've added a statement proving the monotonic decrease in marginal variance to the exact posterior marginal variance as a function of $i$.\nIf anything remains unclear, we are happy to elaborate further during the discussion period!\n\n### Weaknesses\n\n> - Datasets have not been properly cited. In particular, a general citation to UCI doesn’t give credit to the creators of each dataset. License information was also not stated.\n\nWe have added citations for each of the datasets we have used and added license information to the supplementary.\n\n> - I found the discussion of information gain not well-motivated. The definition of information gain was assumed to be known, and it wasn’t clear from the text why maximizing this quantity would lead to a good quality approximation. If the only goal is to suggest projection onto the largest $i$-eigenvalues, there are many results in the matrix approximation literature that also motivate this approach. Moreover, there is prior work in the Nyström approximation literature considering the top $i$-eigenvalues as a subspace for performing inference, for example Ferrari-Trecate et al, https://proceedings.neurips.cc/paper/1998/hash/55c567fd4395ecef6d936cf77b8d5b2b-Abstract.html.\n\nBased on your feedback, we've reformulated the information-theoretic motivation for projecting onto the top $i$-eigenvalue subspace as follows. We've replaced Lemma S2 with an arguably more intuitive result, namely that an \"optimal\" linear combination of the data should maximally reduce our uncertainty about the latent function at the training data or equivalently maximally increase the divergence of our updated belief from the prior. More formally, it holds that the actions maximizing the information gain, defined as the difference in entropy between prior and posterior, i.e.\n\n$$\nS_i = \\text{argmax}_{S \\in \\mathbb{R}^{n \\times i}} ( H(f(X)) - H(f(X) \\mid S^{\\top} y) ) = \\text{argmax}\\_{S \\in \\mathbb{R}^{n \\times i}} \\text{KL}(p(f \\mid S^\\top y)\\ || \\ p(f))\n$$\n\nare given by the top-$i$ eigenvectors of $\\hat{K}$, where $y \\in \\mathbb{R}^n$. This is different from the previous formulation in equation (10) in that here we are directly reasoning about the latent function of interest, $f$. \nIn Bayesian information-theoretic formulations of active learning, choosing observations that minimize posterior entropy (or a myopic approximation defined by the expected information gain) has a long history (e.g. MacKay (1992), Section 2 of Houlsby et al. (2011)). Here, rather than choosing observations directly, we generalize this idea to choosing *linear combinations of observations*.\nWe've included this improved result and its motivation in the final version.\nWe've also added a citation to the work by Ferrari-Trecate to make the connection to prior work explicit.\n\n- MacKay, D. J. C. (1992). Information-Based Objective Functions for Active Data Selection. Neural Computation, 4(4), 590–604. https://doi.org/10.1162/neco.1992.4.4.590\n- Houlsby, N., Huszár, F., Ghahramani, Z., & Lengyel, M. (2011). Bayesian Active Learning for Classification and Preference Learning. arXiv. https://doi.org/10.48550/arXiv.1112.5745\n\n> - The phrasing in the caption of Figure 4 conflates credible and confidence intervals. The accuracy of the credible interval should be with reference to the (exact) posterior and cannot be shown by frequentist coverage. To be clear, I think checking the empirical coverage of the methods is reasonable, but don’t think you can draw conclusions about the accuracy of the credible interval from this.\n\nThank you for pointing this out to us. We have modified the caption to describe that we are measuring the empirical coverage of the credible interval of IterGP-Opt and SVGP and have removed any claim about the accuracy of the credible intervals.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "2x9ZDcYvId", "reviewer": "Reviewer_wkY4", "summary": "The authors extend the work of Wenger et al. (2022) on \"computational uncertainty\" in Gaussian process regression models, which introduced the IterGP class of approximations. These approximations treat the limited computation in approximate GP methods as a source of uncertainty, resulting in uncertainty predictions that provably do not suffer from overconfidence. However, the original IterGP paper did not provide a straightforward method for model/hyperparameter selection: this paper addresses that by discussing two approaches for model selection, eventually settling on a new evidence lower bound using the variational family implied by the IterGP approximation, with different proposals for the \"action matrix.\" They compare against SVGP and SGPR, two popular baselines for approximate Gaussian process regression.", "questions": "In order of relative importance:\n1. Following from my main weaknesses, could you confirm whether you used the default GPyTorch settings in terms of faster computations? Apologies again if I missed this!\n2. Could you confirm why SGPR was not used for Road and Power? My understanding is that the memory cost, which would be the main limitation, should be roughly on the same order as for IterGP-Opt, or am I incorrect?\n3. In Fig. 1, could you confirm how you treated the inducing point locations? Although it shows your point nicely, I find it a bit misleading, as in my experience SVGP would never really optimize its inducing locations to be in a (relatively) large gap between datapoints, and so I was surprised to see an inducing location at about 0.6 on the x-axis.\n4. Out of interest, since the ELBO for IterGP-Opt is not stochastic, could L-BFGS-B be used there to speed up optimization?\n5. It's not clear to me how the statement in lines 111-113 would be implied by the equation above it - could you elaborate please?\n\nBelow I'll list some minor questions and typos:\n- Could you please change the colors used in the experiment plots and tables? I am unfortunately red-green colorblind, and I found it very difficult to disentangle the colors between SVGP and SGPR and again between IterGP-Opt and -CG. One example of colorblind-friendly colors in matplotlib can be `plt.style.use('tableau-colorblind10')`. Thanks!\n- Missing space in line 49 in \"introduced,a\"\n- Admittedly really nitpicky, but could you reframe Eq. 3 and so forth as a maximization problem, or be more clear that you're minimizing the negative ELBO? I think it makes more sense to keep the ELBO as a lower bound to the log marginal likelihood, and so maximizing makes more sense\n- line 99: \"we optimizer\" -> \"we optimize\"", "limitations": "In my opinion the authors have fairly assessed the limitations of their work, namely the inability to do minibatching and the method's limitation to non-conjugate likelihoods.\n\n-----\nPost-discussion update\n-----\nI would like to thank the authors for an engaging discussion. Overall, I was pleased with the quality of the responses and the additional experiments provided. While the work is perhaps not the most ''exciting'' in terms of its empirical results and it has some limitations (e.g., not being able to be minibatched, as pointed out by another reviewer), I don't really think it's useful to judge a paper along these lines, and even so I think this paper represents a solid contribution and has the potential to pave the way for important developments further down the line. Therefore, I am raising my score accordingly.\n\nI would also like the opportunity to address the last comments left by the authors. As I was not able to respond to these before the author-reviewer discussion period ended, I have not taken these into account for my final score. However, I hope that the authors will take them into consideration when they update their paper.\n\n> As we write in the caption of Table 1 in the original submission, the reported number of epochs is determined as follows: \"[...] We report the epoch where each method obtained the lowest NLL, and all performance metrics (NLL, RMSE, and wall-clock time) are measured at this epoch. [...]\". The rebuttal PDF did not contain this full caption due to the limited space available, but the final version will, of course. The total number of epochs was determined by resource constraints (1000 epochs for Adam, 100 for LBFGS, except for Power; see also lines 235-238 and our rebuttal response). As Figure R1 (rebuttal pdf) shows this results in roughly similar overall training time between for example SGPR trained with LBFGS and IterGP-Opt trained with Adam. We would recommend interpreting Table (R)1 and Figure (R)1 together to understand the generalization performance.\n\nThanks for your clarification here. Perhaps I have misunderstood which NLL you are looking at, but I am concerned that this amounts to \"cheating\" - as I understand it you are looking at test metrics to determine when to stop the optimization process, which of course you wouldn't have access to in a real setting. It would be better to either fix the number of epochs or determine some form of early stopping that only relies on training or validation data.\n\n> We would like to politely disagree with the characterization that the choice of optimizer and its parameters doesn't determine the final result (even in theory). The log-marginal likelihood (also the ELBO) as a function of the kernel hyperparameters is generally non-convex and often has multiple (local) minima (for an example see Figure 5.5 of Rasmussen and Williams (2008)). Which modes an optimizer converges to depends on the initialization, the choice of optimizer and the learning rate (schedule / line search). \n\nYes, I am in agreement here - I had meant my comment in a highly idealized sense where we could assume that different optimizers would reach the same optimum, but it was sloppy at best and likely unrealistic in many settings.\n\n> The reason we would still recommend using Adam over L-BFGS for IterGP-Opt is the increased resource usage. While the performance is slightly better on both datasets as you correctly point out, the runtime is also larger (e.g. Parkinsons ~=1.3x, Bike 1.8x longer) and the memory usage is higher. For IterGP-Opt one optimizes the kernel hyperparameters and the action entries, which in our case amounts to parameters. As, for example, the PyTorch documentation on L-BFGS points out, L-BFGS requires param_bytes * (history_size+1) amount of memory, which for a typical history size (e.g. PyTorch's default of 100) adds significant overhead, as compared to a single gradient, which only requires param_bytes.\n\nI'm ok with this logic favoring Adam over L-BFGS for memory reasons, as long as it is clearly stated. But I do think it's a bit hard to say that Adam is necessarily faster for IterOpt-GP - it seems fairly close, and in my experience (admittedly, with SGPR) L-BFGS has a tendency to optimize much quicker initially but then spend a lot of time searching for the optimum. It might be instructive to show the training curves as well for both optimizers, particularly as the early stopping I've mentioned above might confound these determinations.\n\nFinally, I did want to comment on something I didn't pick up on from a previous response you gave but caught during the reviewer discussion phase:\n> As per the recommendation in Section 4.1 of Ober et al (2022), in our original experiments we used a (single precision) Cholesky decomposition both for SGPR and SVGP\n\nWhere in Ober et al (2022) does it say they used single precision for SGPR and SVGP? I could not find any mention of single vs. double precision, and I believe it's fairly common for SGPR and SVGP to still be done in double precision, despite the default single precision in GPyTorch.", "rating": 7, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "The work presents a novel method for model selection and linear-time inference for the IterGP approximations originally proposed by Wenger et al. (2022). They discuss a number of viable options, and overall settle primarily on a newly derived ELBO with a sparse parameterization of an \"action matrix,\" although they do also provide experiments for conjugate gradient-based \"actions.\"\n\nIn my opinion, the quality of the work is largely high, is clearly written, and provides practical steps forward for the GP community. The paper is enjoyable to read and it was great to understand the various avenues that the authors explored to achieve practical model selection. Generally, the authors also seem relatively open about the limitations of the method and don't try to over-claim what they have achieved.", "weaknesses": "My main concerns with the paper concern its in its evaluation, particularly with regards to the sparse variational methods. I do not recall seeing these problems on this dataset, although I have admittedly found Parkinson's more problematic, but not in the way shown (I usually find that it just fails with a bad Cholesky decomposition instead of showing divergent behavior). \n\nIn my experience, this is because while GPyTorch is a state-of-the-art package for conjugate gradient-based methods, it is problematic when it comes to its implementations of sparse variational methods. I believe these instabilities are rather due to GPyTorch's use of Lanczos and conjugate gradient algorithms when computing covariance matrix decompositions and log probabilities, in addition to the use of float32. In my experience, stochastic variational methods are very sensitive to numerical accuracy and so should only be computed using the full Cholesky decomposition and float64. Therefore, since I did not see mention of this (although apologies if I have missed it!) I would recommend changing the GPyTorch defaults to use Cholesky and float64, or using a package that does not rely on approximate matrix decompositions such as GPflow or GPJax. \n\nMoreover, one of the key attributes of SGPR is that because the loss is not stochastic, it can be optimized (for both hyperparameters and inducing locations) using second-order optimizers such as L-BFGS-B, as is commonly done. This is typically much faster than Adam, and essentially allows optimization to be parameter-free (as the usual defaults generally work well). It would have been good to compare to this, as is common in the literature, although I do commend the parameter search that the authors did for Adam optimization.\n\nFinally, it would be interesting to see how the resultant ELBOs compare, as well as comparisons of ELBOs and hyperparameters to exact GPs on smaller datasets where feasible (I think that about 15-20k datapoints should be feasible given the compute the authors had access to). Indeed, a lot of the claims throughout the paper center around comparing IterGP to SVGP/SGPR with respect to the exact GP, and so it seems a bit strange that there aren't really any experimental comparisons with respect to exact GPR!", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "### Summary\n\nThank you for your feedback on our paper! We appreciate your detailed recommendations for improvement. Based on these we've added a new set of experiments where we trained SGPR and an exact (Cholesky)GP with LBFGS in double precision (see the rebuttal PDF).\nWe believe we have addressed all of your questions fully below, but should anything remain unclear, we are happy to respond during the discussion period.\nShould our response have influenced your initial evaluation of our work, we would be highly appreciative if you would consider updating your score.\nThank you!\n\n### Weaknesses\n\n> My main concerns with the paper concern its in its evaluation, particularly with regards to the sparse variational methods. I do not recall seeing these problems on this dataset, although I have admittedly found Parkinson's more problematic, but not in the way shown (I usually find that it just fails with a bad Cholesky decomposition instead of showing divergent behavior). In my experience, this is because while GPyTorch is a state-of-the-art package for conjugate gradient-based methods, it is problematic when it comes to its implementations of sparse variational methods. I believe these instabilities are rather due to GPyTorch's use of Lanczos and conjugate gradient algorithms when computing covariance matrix decompositions and log probabilities, in addition to the use of float32. In my experience, stochastic variational methods are very sensitive to numerical accuracy and so should only be computed using the full Cholesky decomposition and float64. Therefore, since I did not see mention of this (although apologies if I have missed it!) I would recommend changing the GPyTorch defaults to use Cholesky and `float64`, or using a package that does not rely on approximate matrix decompositions such as GPflow or GPJax. \n\nWe would like to clarify that we closely followed the recommendations of Ober et al. (2022) and computed any (covariance) matrix decompositions for SGPR and SVGP with Cholesky using an iteratively increasing additive jitter should the decomposition fail. We did *not* use Lanczos or CG.\nThe iteratively added jitter could explain the observed divergent behavior instead of an outright failure on Parkinson's.\nWe would like to politely point out that GPyTorch does not use Lanczos or CG for variational methods, but rather a Cholesky decomposition (see the source code for [SGPR](https://github.com/cornellius-gp/gpytorch/blob/283105b4f71b0e97fb2dbbf34c95396c7a88bd25/gpytorch/kernels/inducing_point_kernel.py#L65) and [SVGP](https://github.com/cornellius-gp/gpytorch/blob/283105b4f71b0e97fb2dbbf34c95396c7a88bd25/gpytorch/variational/variational_strategy.py#L203)).\nWe will expand our description of the experiments with this additional information.\nTo alleviate concerns about the use of single precision, we conducted new experiments in double precision as described in the answer to the next question.\n\n> Moreover, one of the key attributes of SGPR is that because the loss is not stochastic, it can be optimized (for both hyperparameters and inducing locations) using second-order optimizers such as L-BFGS-B, as is commonly done. This is typically much faster than Adam, and essentially allows optimization to be parameter-free (as the usual defaults generally work well). It would have been good to compare to this, as is common in the literature, although I do commend the parameter search that the authors did for Adam optimization.\n\nThank you for this suggestion! We've made the following improvements to our main experiment.\nWe trained SGPR on all datasets in double precision using L-BFGS with a Wolfe line search for 100 epochs, based on the recommendations of Ober et al. (2022). We performed a sweep over the initial learning rate $\\text{lr}_{\\text{SGPR}} \\in \\\\{10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}, 1\\\\}$ and repeated each run for five random seeds. Both in our original experiments using single precision and Adam and the new ones, we computed Cholesky decompositions with adaptively increasing kernel jitter as per Section 4.1 of Ober et al (2022). \nWe find that this indeed increases its performance as the updated Table R1 in the attached rebuttal PDF shows. SGPR now outperforms SVGP on all small to medium datasets as expected, but IterGP-Opt still matches or outperforms SGPR on all datasets, except \"Bike\", as both Table R1 and Figure R1 in the attached rebuttal PDF show.\n\n> Finally, it would be interesting to see how the resultant ELBOs compare, as well as comparisons of ELBOs and hyperparameters to exact GPs on smaller datasets where feasible (I think that about 15-20k datapoints should be feasible given the compute the authors had access to). Indeed, a lot of the claims throughout the paper center around comparing IterGP to SVGP/SGPR with respect to the exact GP, and so it seems a bit strange that there aren't really any experimental comparisons with respect to exact GPR!\n\nBased on your suggestion we've added results for exact (Cholesky)GPs to our experiments (see Table R1 of the rebuttal PDF). We trained a CholeskyGP on Parkinson's and Bike with LBFGS (using a Wolfe line search) in double precision for 100 epochs. We performed a parameter sweep over the initial learning rate $\\mathrm{lr}_{\\text{CholeskyGP}} \\in \\\\{10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}, 1\\\\}$ and repeated each run five times for different seeds.\nPerhaps unsurprisingly, the exact CholeskyGP outperforms all approximations on those small datasets.\nWe also added the CholeskyGP to the hyperparameter plot in Figure S4 of the supplementary of our paper. The exact GP tends to learn a larger outputscale and smaller observation noise than the approximate methods. Unfortunately, we could not add this to the rebuttal pdf due to the space constraint of one page only.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank all reviewers for their time and effort in reviewing our work!
We feel that your helpful suggestions enabled us to improve our paper, in particular, we've added the following additional experiments (see the attached PDF and our individual rebuttals below):

- **Generalization experiments** (Table R1, Figure R1)
    - SGPR trained with LBFGS in double precision on all small and medium sized datasets.
    - IterGP-Opt trained with LBFGS in double precision on all small and medium sized datasets.
    - CholeskyGP baseline trained with LBFGS in double precision on all datasets where this is possible given memory constraints.
- **Experiment on inducing point placement and posterior vs predictive uncertainty of SVGP in higher dimensions** (Figure R2).

We responded to all your questions in individual answers below. If anything remains unclear please let us know!

Meta-review: A summary of the strengths and weaknesses based on the reviews and rebuttal (including the follow-up discussion and that among the reviewers) is provided below:


**STRENGTHS**

1. This paper proposes a novel linear-time IterGP-Opt model by minimizing the (negative) ELBO training loss that optimizes the hyperparameters as well as parameters controlling the approximation quality via variational inference. Such models do not suffer from overconfidence in terms of the uncertainty of their predictions.

2. The paper is well-written.


**WEAKNESS**

1. As the authors have highlighted, the proposed IterGP-Opt is not amenable to stochastic minibatching and hence cannot achieve constant time/memory per iterative step. Nonetheless, the reviewers think that this paper can potentially inspire follow-up works with stochastic minibatching.

2. An initial concern was the unfair comparison with SGPR/DTC as little effort has been put into getting SGPR to perform competitively well. The authors' rebuttal has resolved this issue.

3. The empirical evaluation is lacking. In particular, the baselines are at least a decade old and more powerful derivatives of these baselines and others have since been developed. See, for example, the following work (and references therein) that can also perform hyperparameter optimization/model selection and inference, and the approximation quality can be improved by trading off time efficiency: A Distributed Variational Inference Framework for Unifying Parallel Sparse Gaussian Process Regression Models. ICML 2016.

If not empirically comparing with the above, the authors should at least give a qualitative comparison and its implications on the resulting performance difference.

Despite the weaknesses, we think that the novelty and significance of this work (specifically, point 1 in strength) are sufficiently strong and outweigh such weaknesses.
Decision: Accept (poster)