Computation-Aware Gaussian Processes: Model Selection And Linear-Time Inference

Authors: jonathan  wenger
kaiwen  wu
philipp  hennig
jacob r gardner
geoff  pleiss
john p cunningham
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 22:38:56

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 computation 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.

Main paper: 
 introduction: 
 Gaussian Processes (GPs) remain a popular probabilistic model class, despite the challenges in scaling them to large datasets. Since both computational and memory resources are limited in practice, approximations are necessary for both inference and model selection. Among the many approximation methods, perhaps the most common approach is to map the data to a lower-dimensional 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 a lower-dimensional latent inducing point space (e.g. SoR, DTC, FITC [3] , SGPR [4] , SVGP [5] ), or implicit-by using an iterative numerical method (e.g. CGGP [6] [7] [8] [9] [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).While effective and widely used in practice, the inevitable approximation error adversely impacts predictions, uncertainty quantification, and ultimately downstream decision-making. Many proposed methods come with theoretical error bounds [e.g. 2, [11] [12] [13] [14] , offering insights into the scaling and asymptotic properties of each method. However, theoretical bounds often require too many assumptions about the data-generating process to 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 (NeurIPS 2024). . Hyperparameters for each model were optimized using model selection strategies specific to each approximation. The posterior predictive given the datagenerating 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 data-generating 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 applications of GPs such as Bayesian optimization [e.g. variance starvation of RFF, 16], and manifests itself even in state-of-the-art variational methods like SVGP. SVGP, because it treats inducing variables as "virtual observations", can be overconfident at the locations of the inducing points if they are not in close proximity to training data, which becomes increasingly likely in higher dimensions. This phenomenon can be seen in a toy example in Figure 1 , 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 are a central issue in inference, but they are exacerbated in model selection, where errors compound and result in biased selections of hyperparameters [12, 17, 18] . Continuing the example, SVGP has been observed to overestimate the observation noise [18] , which can lead to oversmoothing. This issue can also be seen in Figure 1 , where the SVGP model produces a smoother posterior mean than the exact (Cholesky)GP and attributes most variation from the posterior 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 et al. [19] introduced computation-aware Gaussian processes (CaGP), a class of GP approximation methods which-for a fixed set of hyperparameters-provably does not suffer from overconfidence. Like SVGP and the other approximations mentioned above, CaGP also relies 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 approximation error due to 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.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 theoretical 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 linear-time inference, in a way that is amenable to hardware acceleration. We optimize these actions end-to-end alongside the hyperparameters with respect to a custom training loss, to optimally retain as much information from the data as possible given a limited computational budget. The resulting hyperparameters 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 is strongly competitive on large-scale data with state-of-the-art 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 a few hours on a single GPU without adversely impacting uncertainty quantification. 
background: 
 We aim to learn a latent function mapping fromX ⊆ R d to Y ⊆ R given a training dataset X = (x 1 , . . . , x n ) ∈ R n×d of n inputs x j ∈ R d and corresponding targets y = (y 1 , . . . , y n ) T ∈ R n .Gaussian Processes A Gaussian process f ∼ GP(µ, K θ ) is a stochastic process with mean function µ and kernelK θ such that f = f (X) = (f (x 1 ), . . . , f (x n )) T ∼ N (µ, K θ ) is jointly Gaussian with mean µ i = µ(x i ) and covariance K ij = K θ (x i , x j ).The kernel K θ depends on hyperparameters θ ∈ R p , which we omit in our notation. Assuming y | f (X) ∼ N f (X), σ 2 I , the posterior is a Gaussian process GP(µ ⋆ , K ⋆ ) where the mean and covariance functions evaluated at a test inputEQUATIONwhere K = K + σ 2 I and the representer weights are defined as v ⋆ = K-1 (y -µ).In model selection, the computational bottleneck when optimizing kernel hyperparameters θ is the repeated evaluation of the negative log-marginal likelihoodℓ NLL (θ) = -log p(y | θ) = 1 2 (y -µ) T K-1 (y -µ)quadratic loss: promotes data fitEQUATIONand its gradient. Computing (2) and its gradient via a Cholesky decomposition has time complexity O(n 3 ) and requires O(n 2 ) memory, which is prohibitive for large n.Sparse Gaussian Process Regression (SGPR) [4] Given a set of m ≪ n inducing points Z = (z 1 , . . . , z m ) T and defining u := f (Z) = (f (z 1 ), . . . , f (z m )) T , SGPR defines a variational approximation to the posterior through the factorizationp ⋆ (f (•) | y) ≈ q(f (•)) = p(f (•) | u) q(u)du, where q(u) is an m-dimensional multivariate Gaussian. The mean and covariance of q(u) (denoted as m := E q (u), Σ := Cov q (u)) are jointly optimized alongside the kernel hyperparameters θ using the evidence lower bound (ELBO) as an objective:EQUATIONEQUATIONThe inducing point locations Z can be either optimized as additional parameters during training or chosen a-priori, typically in a data-dependent way [see e.g. Sec. 7.2 of 20]. Following Titsias [4] , ELBO optimization and posterior inference both require O(nm 2 ) computation and O(nm) memory, a significant improvement over the costs of exact GPs.Stochastic Variational Gaussian Processes (SVGP) [5] SVGP extends SGPR to reduce complexity further to O(m 3 ) computation and O(m 2 ) memory. It accomplishes this reduction by replacing the first term in Equation ( 4) with an unbiased approximationE q(f ) (log p(y | f )) = E q(f ) ( n i=1 log p(y i | f (x i ))) ≈ n E q(f (xi)) (log p(y i | f (x i ))) . Following Hensman et al. [5] , we optimize m, Σ alongside θ, Z through joint gradient updates. Because the asymptotic complexities no longer depend on n, SVGP can scale to extremely large datasets that would not be able to fit into computer/GPU memory.whereC i = S i (S T i KS i ) -1 S T i ≈ K -1is a low-rank approximation of the precision matrix and v i = C i (y -µ) ≈ v ⋆ approximates the representer weights. Both converge to the corresponding exact quantity as the number of iterations, equivalently the downdate rank, i → n. Note that the CaGP posterior only depends on the space spanned by the columns of S i , not the actual matrix (Lemma S4). Finally, the CaGP posterior can be computed in O(n2 i) time and O(ni) 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 y ∈ H K σ in the RKHS H K σ defined by the kernel, the worst-case squared error of the corresponding approximate posterior mean µ y i is equal to the approximate predictive variance (see Theorem S2):EQUATIONThis guarantee is identical to one for the exact GP posterior mean and variance (see Theorem S1), except with the approximate quantities instead. 2 Additionally, it holds that CaGP's marginal (predictive) variance is always larger or equal to the (predictive) variance of the exact GP and monotonically decreasing, i.e.K i (x ⋄ , x ⋄ ) ≥ K j (x ⋄ , x ⋄ ) ≥ K n (x ⋄ , x ⋄ ) = K ⋆ (x ⋄ , x ⋄ ) for i ≤ j ≤ n (see Propo- 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 computation-aware, and we will extend the use of this object to model selection.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We explicitly discuss the limitations of our method in the conclusion in a dedicated paragraph.Guidelines:• The answer NA means that the paper has no 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 a 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 authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a 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 approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it 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 address 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 use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
model selection in computation-aware gaussian processes: 
 Model selection for GPs most commonly entails maximizing the evidence log p(y | θ) as a function of the kernel hyperparameters θ ∈ R p . 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-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 i ∈ R n×i , i.e.EQUATION) is the action matrix with orthonormalized columns. The corresponding likelihood is given byp( ỹ | f (X)) = N ỹ; S ′ i T f (X), σ 2 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 i (see Lemma S4), which is why Equation ( 5) can be written in terms of S i directly rather than using S ′ i . 3Projected-Data Log Marginal Likelihood The reinterpretation of the computation-aware posterior as exact Bayesian inference immediately suggests using evidence maximization for the projected data ỹ = S ′ i T y ∈ R i as the model selection criterion, leading to the following loss (see Lemma S2):EQUATION) Equation ( 9) involves a Gaussian random variable of dimension i ≪ n, and so all previously intractable 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 i without having to orthonormalize, which results in an additional term penalizing 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 i the data are projected onto. Specifically, the quadratic loss term only promotes fitting the projected data ỹ, not all observations y. See Figures S1 and S2 for experimental validation of this claim.ELBO Training Loss Motivated by this observation, we desire a tractable training loss that encourages maximizing the evidence for the entire set of targets y. Importantly though, we need to accomplish this evidence maximization without incurring prohibitive O(n 3 ) computational cost. We define a variational objective, using the computation-aware posterior q i (f | y, θ) in Equation (5) to define a variational family Q := q i (f ) = N (f ; µ i (X), K i (X, X)) | S ∈ R n×i parametrized by the action matrix S. We can then specify a (negative) evidence lower bound (ELBO) as follows:ℓ ELBO CaGP (θ) = ℓ NLL (θ)balances data fit and model complexityEQUATIONThis loss promotes learning the same hyperparameters as if we were to maximize the computationally intractable evidence log p(y | θ) while minimizing the error due to posterior approximation q i (f ) ≈ p ⋆ (f | y, θ). In the computation-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.EQUATIONwhere ṽi = (S T KS) -1 S T (y -µ). For a derivation of this expression of the loss see Lemma S3. If we compare the training loss ℓ ELBO CaGP in Equation ( 10) with the projected-data NLL 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 S2 show on synthetic data. In our experiments, this objective was critical to achieving state-of-the-art performance.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes]Justification: We provide proofs to all theoretical statements in the supplementary material or 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 crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems. • The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
choice of actions: 
 So far we have not yet specified the actions S ∈ R n×i 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 learning, where instead of just observing individual datapoints, we allow ourselves to observe linear combinations of the data y. Taking an information-theoretic perspective [21] , we would then aim to choose actions such that uncertainty about the latent function is minimized. In fact, we show in Lemma S5 that given a prior f ∼ 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 dataEQUATIONare the top-i eigenvectors s 1 , . . . , s i 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 minimize posterior entropy, we could try to do so approximately. The Lanczos algorithm [23] is an iterative method to approximate the eigenvalues 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 i = (y -µ)-Kv i-1 of the method of Conjugate Gradients (CG) [24] when used to iteratively compute an approximation v i ≈ v ⋆ = K-1 (y -µ) 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 Lanczos 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. 4 See Section S3.1 for details. 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 O(n 2 i).Learned Sparse Actions So far in our action choices we 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 to perform model selection, is recognizing that the actions 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 hyperparameters without spending unnecessary budget on computing approximately optimal actions for the current choice of hyperparameters. Specifically, the actions are chosen by optimizing ℓ ELBO CaGP as a function of the hyperparameters θ and the actions S, s.t.EQUATIONEQUATIONNaively this approach introduces an n × i dimensional optimization problem, which in general is computationally prohibitive.To keep the computational cost low and to optimally leverage hardware acceleration via GPUs, we impose a sparse block structure on the actions (see Eq. 14) where each block is a column vector s ′ j ∈ R k×1 with k = n/i entries such that the total number of trainable action parameters, i.e. non-zero entries nnz(S) = k • i = n, equals the number of training data. Due to the sparsity, these actions cannot perfectly match the maximum eigenvector actions. Nevertheless, we see in Figure 2 that optimizing these sparse actions in conjunction with hyperparameter optimization produces a nontrivial alignment with optimal action choice minimizing posterior entropy. Importantly, the sparsity constraint not only reduces the dimensionality 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.We give algorithms both for iteratively constructed dense actions (Algorithm S1), as used in CaGP-CG, and for sparse batch actions (Algorithm S2), as used for CaGP-Opt, in the supplementary material. 5 The time complexity is O(n 2 i) for dense actions and O(ni max(i, k)) for sparse actions, where k is the maximum number of non-zero entries per column of S i . Both have the same linear memory requirement: O(ni). Since the training loss ℓ ELBO CaGP only involves terms that are also present in the posterior predictive, both model selection and predictions incur the same complexity.Computational Uncertainty and Probabilistic Numerics All CaGP variants discussed in this paper fall into the category of probabilistic numerical methods [25] [26] [27] , which aim to quantify approximation error arising from limited computational resources via additional uncertainty about the quantity of interest [e.g. [28] [29] [30] [31] . Specifically, the iterative formulation of CaGP (i.e. Algorithm S1) originally proposed by Wenger et al. [19] employs a probabilistic linear solver [32] [33] [34] [35] .Scalable GPs with Lower-Bounded Log Marginal Likelihoods Numerous scalable GP approximations beyond those in Sections 1 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] [41] [42] learn a variational posterior approximation on a small set of linear functionals applied to the latent GP. However, unlike our method, their resulting approximate posterior is usually prone to underestimating uncertainty in 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] [17] [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] .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 conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We describe the experiments including datasets, number of repeats, hyperparameters, method implementations and hardware in Section 5. 
experiments: 
 We benchmark the generalization of computation-aware GPs with two different action choices, CaGP-Opt (ours) and CaGP-CG [19] , using our proposed training objective in Equation ( 10) on a range of UCI datasets for regression [53] . We compare against SVGP [5] , often considered to be state-of-the-art for large-scale GP regression. Per recommendations by Ober et al. [15] , we also include 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-CG 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 a (0.9, 0.1) split for five random seeds. We used a zero-mean GP prior and a Matérn(3/2) kernel with an outputscale o 2 and one lengthscale per input dimension l 2 j , as well as a scalar observation noise for the likelihood σ 2 , s.t. θ = (o, l 1 , . . . , l d , σ) ∈ R d+2 . We used the existing implementations 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 hyperparameters θ either with Adam [54] for a maximum of 1000 epochs in float32 or with LBFGS [55] for 100 epochs in float64, depending on 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 a batch size of 1024 throughout. We scheduled the learning rate via PyTorch's [56] LinearLR(end factor=0.1) scheduler for all methods and performed a hyperparameter sweep for the (initial) learning rate. All experiments were run on an NVIDIA Tesla V100-PCIE-32GB GPU, except for "Power", where we used an NVIDIA A100 80GB PCIe GPU to have sufficient memory for CaGP-Opt with i = 512. Our exact experiment configuration can be found in Table S1 .Evaluation Metrics 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 3 show generalization performance 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 (except on "Bike"). SGPR remains scompetitive for smaller datasets; however, it does not scale to the largest datasets. There are some datasets and metrics in which specific methods dominate, for example 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-Opt consistently performs either best or second best and scales to over a million datapoints. These results are quite remarkable for numerous reasons. First, CaGP is comparable in runtime to SVGP on large datasets despite the fact that it incurs a linear-time computational complexity while SVGP is constant time. 6 Second, while all of the methods we compare approximate the GP posterior with low-rank updates, CaGP-Opt (with i = 512) here uses half the rank of SGPR/SVGP m = 1024. Nevertheless, CaGP-Opt is able to substantially outperform SVGP even on spatial datasets like 3DRoad where low-rank posterior updates are often poor [57] . These results suggest that CaGP-Opt can be a 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 a "ramp-up" phase, compared to SVGP, but then improves or matches its generalization performance. This gap is particularly large on "Road", where CaGP-Opt is initially worse than SVGP but dominates in the second half of training. In Figure S5 we show that SVGP typically learns larger observation noise than other methods as suggested by previous work [18, 45] and hinted at by observations on synthetic data in Figure 1 and Figure S3 (b). Additionally on larger datasets SVGP also often learns large lengthscales, which in combination with a large observation noise can lead to an oversmoothing effect. In contrast, CaGP-Opt generally learns lower observational noise than SVGP. Of course, learning a 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 as noise.CaGP Improves Uncertainty Quantification Over SVGP A 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 α credible interval of the CaGP-Opt posterior; i.e. 4 compares the 95% coverage error 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.ε α coverage = |α -1 ntest ntest i=1 1(y i ∈ I α q(xi) )|. Figure 
conclusion: 
 In this work, we introduced strategies for model selection and posterior inference for computationaware Gaussian processes, which scale linearly with the number of training data rather than quadratically. The key technical innovations being a sparse projection of the data, which balances minimizing posterior entropy and computational cost, and a scalable way to optimize kernel hyperparameters, both of which are amenable to GPU acceleration. All together, these advances enable competitive or improved performance over previous approximate inference methods on large-scale datasets, in terms of generalization and uncertainty quantification. Remarkably, our method outperforms SVGP-often considered the de-facto GP approximation standard-even when compressing the data into a space of half the dimension of the variational parameters. Finally, we also demonstrate 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 as 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 is faster.) CaGP-Opt nevertheless requires larger GPUs than SVGP on datasets with more than a 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, is limited to GP regression with a conjugate observational noise model. We leave extensions to classification and other non-conjugate likelihoods as future work.Outlook and Future Work An immediate consequence of this work is the ability to apply computation-aware Gaussian processes to "real-world" problems, as our approach solves CaGP's open problems of model selection and scalability. Looking forward, an exciting future vision is a general framework for problems involving a Gaussian process model with a 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. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: In Section 5, we describe the specific GPUs we 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, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: Due to the methodological nature of this work, there are no 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.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
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 specific negative societal impact that we feel must be explicitly highlighted here.Guidelines:• 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 negative 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 specific 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 a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 to 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 or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation 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 ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper does not release data or models with a high risk for misuse.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to 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 not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: 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? 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 or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• Depending on the country in which research is conducted, IRB approval (or equivalent) 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. 
 
 Appendix: 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, or results from, the main paper are stated as is. Lemma S1 (CaGP Inference as Exact Inference Given a Modified Observation Model) Given a Gaussian process prior f ∼ GP(µ, K) and training data (X, y) the computation-aware GP posterior GP(µ i , K i ) (see Equation ( 5)) with linearly independent and fixed actions S is equivalent to an exact batch GP posterior (f | ỹ) given data ỹ = S ′ T y observed according to the likelihood ỹ | f (X) ∼ N S ′ T f (X), σ 2 I , where S ′ = S chol(S T S) -T .Proof. First note that by definition S ′ has orthonormal columns, sinceS ′ T S ′ = (S chol(S T S) -T ) T S chol(S T S) -T = chol(S T S) -1 S T S chol(S T S) -T = L -1 LL T L -T = I.Now by basic properties of Gaussian distributions, we have for arbitraryX ⋄ ∈ R n⋄×d that ỹ f (X ⋄ ) ∼ N S ′ T µ(X) µ(X ⋄ ) , S ′ T K(X, X)S ′ + σ 2 S ′ T S ′ S ′ T K(X, X ⋄ ) K(X ⋄ , X)S ′ K(X ⋄ , X ⋄ )is jointly Gaussian, where we used that I = S ′ T S ′ .Therefore we have that(f (X ⋄ ) | ỹ) ∼ N (µ i (X ⋄ ), K i (X ⋄ , X ⋄ )) with µ i (X ⋄ ) = µ(X ⋄ ) + K(X ⋄ , X)S ′ (S ′ T (K(X, X) + σ 2 I)S ′ ) -1 ( ỹ -S ′ T µ), = µ(X ⋄ ) + K(X ⋄ , X)S(S T (K(X, X) + σ 2 I)S) -1 S T (y -µ) K i (X ⋄ , X ⋄ ) = K(X ⋄ , X ⋄ ) -K(X ⋄ , X)S ′ (S ′ T (K(X, X) + σ 2 I)S ′ ) -1 S ′ T K(X, X ⋄ ) = K(X ⋄ , X ⋄ ) -K(X ⋄ , X)S(S T (K(X, X) + σ 2 I)S) -1 S T K(X, X ⋄ )which is equivalent to the definition of the CaGP posterior in Equation ( 5). This proves the claim.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 is 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 f ∼ GP(0, K) be a Gaussian process with kernel K and define the observed process y(•) = f (•) + σε(•) where ε ∼ GP(0, δ) is white noise with noise level σ 2 > 0, i.e. δ(x,x ′ ) = 1(x = x ′ ).Consequently, the covariance kernel of the data-generating process y(•) is given by K σ (x, x ′ ) := K(x, x ′ ) + σ 2 δ(x, x ′ ) and we denote the corresponding RKHS as H K σ . For exact Gaussian process inference, the pointwise (relative) worst-case squared error of the posterior 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 x 1 , . . . , x n ∈ X, the GP posterior GP(µ ⋆ , K ⋆ ) satisfies for anyx ̸ = x j that sup y∈H K σ ,∥y∥ H K σ ≤1 (y(x ⋄ ) -µ y ⋆ (x ⋄ )) 2 error of posterior mean = K ⋆ (x ⋄ , x ⋄ ) + σ 2 predictive variance (S15)If σ 2 = 0, then the above also holds for x ∈ x j .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 computationaware.Theorem S2 (Worst Case Error Interpretation of CaGP Variance [19] ) Given a set of training inputs x 1 , . . . , x n ∈ X, the CaGP posterior GP(µ i , K i ) satisfies for anyEQUATIONerror of approximate posterior mean= K i (x ⋄ , x ⋄ ) + σ 2 approximate predictive variance (S16)If σ 2 = 0, then the above also holds for x ∈ x j .Proof. See Theorem 2 of Wenger et al. [19] .While SVGP also admits a 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 data-generating function y ∈ H K σ and the approximate posterior mean µ SVGP (x ⋄ ). It only involves a worst-case error term over the unit ball in the RKHS of the approximate kernelQ σ ≈ K σ .Theorem S3 (Worst Case Error Interpretation of SVGP Variance [61] ) Given a set of training inputs x 1 , . . . , x n ∈ X and (fixed) inducing points Z ∈ R m×d , the optimal variational posterior GP(µ ⋆ SVGP , K ⋆ SVGP ) of SVGP is given byµ ⋆,y SVGP (x ⋄ ) = K(x ⋄ , Z)(σ 2 K(Z, Z) + K(Z, X)K(X, Z)) -1 K(Z, X)y(X) K ⋆ SVGP (x ⋄ , x ′ ⋄ ) = K(x ⋄ , x ′ ⋄ ) -Q(x ⋄ , x ′ ⋄ ) + K(x ⋄ , Z)(K(Z, Z) + σ -2 K(Z, X)K(X, Z)) -1 K(Z, x ′ ⋄ )where Q(x, x ′ ) = K(x, Z)K(Z, Z) -1 K(Z, x ′ ) is the Nyström approximation of the covariance function K(x, x ′ ) (see Eqns. ( 25) and ( 26) of Wild et al. [61] ). The optimized SVGP posterior satisfies for any x ̸ = x j thatsup y∈H Q σ ,∥h∥ H Q σ ≤1 y(x ⋄ ) -µ ⋆,y SVGP (x ⋄ )2 error of exact posterior mean assuming y(•) is in the RKHS of the approximate kernelQ σ + sup f ∈H K ,∥f ∥ H K ≤1 f (x ⋄ ) -K(x ⋄ , Z)K(Z, Z) -1 f (Z) 2error of exact posterior mean given noise-free observations at inducing points= K ⋆ SVGP (x ⋄ , x ⋄ ) + σ 2 approximate predictive variance (S17)If σ 2 = 0, then the above also holds for x ∈ x j .Proof. See Theorem 6 of Wild et al. [61] .Proposition S1 (CaGP's Variance Decreases Monotonically with the Number of Iterations) Given a training dataset of size n, let GP(µ i , K i ) be the corresponding CaGP posterior defined in Equation ( 5) where i ≤ n denotes the downdate rank / number of iterations and assume the CaGP actions S i ∈ R n×i are linearly independent. Then it holds for arbitrary x ⋄ ∈ X and i≤ j ≤ n, that K i (x ⋄ , x ⋄ ) ≥ K j (x ⋄ , x ⋄ ) ≥ K n (x ⋄ , x ⋄ ) = K ⋆ (x ⋄ , x ⋄ ) (S18)where K ⋆ (x ⋄ , x ⋄ ) is the variance of the exact GP posterior in Equation (1).Proof. Wenger et al. [19] originally defined the approximate precision matrixC i = i ℓ=1 1 η ℓ d ℓ d T ℓ = i ℓ=1dℓ dT ℓ as a sum of rank-1 matrices and show that this definition is equivalent to the batch form C i = S i (S T i KS i ) -1 S T i weK i (x ⋄ , x ⋄ ) = K(x ⋄ , x ⋄ ) -K(x ⋄ , X)C i K(X, x ⋄ ) = K(x ⋄ , x ⋄ ) - i ℓ=1 K(x ⋄ , X) dℓ dT ℓ K(X, x ⋄ ) = K(x ⋄ , x ⋄ ) - i ℓ=1 (K(x ⋄ , X) dℓ ) 2 ≥ K(x ⋄ , x ⋄ ) - j ℓ=1 (K(x ⋄ , X) dℓ ) 2 since i ≤ j ≥ K(x ⋄ , x ⋄ ) - n ℓ=1 (K(x ⋄ , X) dℓ ) 2 = K(x ⋄ , x ⋄ ) -K(x ⋄ , X)C n K(X, x ⋄ ) = K ⋆ (x ⋄ , x ⋄ )where the last equality follows from the fact that S n ∈ R n×n has rank n and thereforeC n = S n (S T n KS n ) -1 S T n = S n S -1 n K-1 (S n S -1 n ) T = K-1 . S2 Training Losses S2.1 Projected-Data Log-Marginal Likelihood Lemma S2 (Projected-Data Log-Marginal Likelihood)Under the assumptions of Lemma S1, the projected-data log-marginal likelihood is given byℓ NLL proj (θ) = -log p( ỹ | θ) = -log N ỹ; S ′ i T µ, S ′ i T KS ′ i = 1 2 (y -µ) T S i (S T i KS i ) -1 S T i (y -µ) + log det(S T i KS i ) -log det(S T i S i ) + i log(2π)Proof. By the same argument as in Lemma S1 we obtain thatℓ NLL proj (θ) = -log p( ỹ | θ) = -log N ỹ; S ′ i T µ, S ′ i T KS ′ i = 1 2 (S ′ i T y -S ′ i T µ) T (S ′ i T KS ′ i ) -1 (S ′ i T y -S ′ i T µ) + log det(S ′ i T KS ′ i ) + i log(2π) = 1 2 (y -µ) T S ′ i (S ′ i T KS ′ i ) -1 S ′ i T (y -µ) + log det(S ′ i T KS ′ i ) + i log(2π) Since S ′ i = S i L -T, where L -T = chol(S T i S i ) -T is the orthonormalizing matrix, L -T cancels in the quadratic loss term, giving= 1 2 (y -µ) T S i (S i T KS i ) -1 S T i (y -µ) + log det(S ′ i T KS ′ i ) + i log(2π)and finally we can decompose the log-determinant into a difference of log-determinants= 1 2 (y -µ) T S i (S i T KS i ) -1 S T i (y -µ) + log det(S T i KS i ) -2 log det(L) + i log(2π) = 1 2 (y -µ) T S i (S i T KS i ) -1 S T i (y -µ) + log det(S T i KS i ) -log det(LL T ) + i log(2π)which using LL T = S T i S i completes the proof.Lemma S3 (Evidence Lower Bound Training Loss) Define the variational familyQ := q(f ) = N (f ; µ i (X), K i (X, X)) | S ∈ R n×i (S19)then the evidence lower bound (ELBO) is given byℓ ELBO CaGP (θ) = -log p(y | θ) + KL(q(f ) ∥ p(f | y)) = -E q (log p(y | f )) + KL(q(f ) ∥ p(f )) = 1 2 1 σ 2 ∥y -µ i (X)∥ 2 2 + n j=1 K i (x j , x j ) + (n -i) log(σ 2 ) + n log(2π) + ṽT i S T KS ṽi -tr((S T KS) -1 S T KS) + log det(S T KS) -log det(S T S)where ṽi = (S T KS) -1 S T (y -µ) are the "projected" representer weights.Proof. The ELBO is given by-ℓ ELBO CaGP (θ) = E q (log p(y | f )) -KL(q(f ) ∥ p(f )). We first compute the expected log-likelihood term.E q (log p(y | f )) = E q - 1 2 1 σ 2 (y -f ) T (y -f ) + log det(σ 2 I n×n ) + n log(2π) = - 1 2 1 σ 2 E q (y -f ) T (y -f ) + n log(σ 2 ) + n log(2π) Now using E x T Ax = E(x) T A E(x) + tr(A Cov(x)), we obtain = - 1 2 1 σ 2 ∥y -µ i (X)∥ 2 2 + tr(K i (X, X)) + n log(σ 2 ) + n log(2π) = - 1 2   1 σ 2 ∥y -µ i (X)∥ 2 2 + n j=1 K i (x j , x j ) + n log(σ 2 ) + n log(2π)  Since both q(f ) and the prior p(f ) are Gaussian, the KL divergence term between them is given byKL(q(f ) ∥ p(f )) = 1 2 (µ i (X) -µ(X)) T K -1 (µ i (X) -µ(X)) + log det(K) det(K i (X, X)) + tr(K -1 K i (X, X)) -n = 1 2 (KC i (y -µ(X))) T K -1 KC i (y -µ(X)) -log det(K -1 K i (X, X)) + tr(I n×n -C i K) -n = 1 2 (y -µ(X)) T C i KC i (y -µ(X)) -log det(I n×n -C i K) + tr(I n×n -C i K) -n = 1 2 ṽT i S T KS ṽi -log det(I n×n -C i K) + tr(I n×n -C i K) -n = 1 2 ṽT i S T KS ṽi -log det(I n×n -C i K) -tr((S T KS) -1 S T KS)Next, we use the matrix determinant lemma When optimizing hyperparameters with respect to the projected-data log-marginal likelihood, CaGP-CG completely overestimates the noise scale, which leads to increasingly worse generalization performance. In comparison, the ELBO training loss leads to a much better overall fit.det(A + U V T ) = det(I m + V T A -1 U ) det(A): = 1 2 ṽT i S T KS ṽi -log det(I i×i -(S T KS) -1 S T KS) -tr((S T KS) -1 S T KS) = 1 2 ṽT i S T KS ṽi -log det((S T KS) -1 (S T KS -S T KS)) -tr((S T KS) -1 S T KS) = 1 2 ṽT i S T KSWe 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 ′ ∈ R n×i be two action matrices, each of which consists of non-zero and linearly independent action vectors, such that their column spaces are identical, i.e.EQUATIONthen the corresponding CaGP posteriors GP(µ i , K i ) and GP(µ ′ i , K ′ i ) are equivalent.Proof. By assumption (S20) there exists W ∈ R i×i such that S ′ = SW . Since action vectors are assumed to be linearly independent and non-zero, it holds that i = rank(S ′ ) = rank(SW ) = rank(W ), where the last equality follows from standard properties of the matrix rank. Therefore W is invertible, and we have thatC ′ = S ′ (S ′ T KS ′ ) -1 S ′ T = SW (W T S T KSW ) -1 W T S T = S(S T KS) -1 S T = C.Proof. Let ỹ := S T y and f := f (X). By assumption, we have that f ∼ N (µ, K). Recall that the entropy of a Gaussian random vector f ∼ N (m, S) is given by H(f ) = 1 2 log det(S) + n log(2πe) . Now since the covariance function of the computation-aware posterior in Equation ( 5) does not depend on the targets y, neither does its entropy H p(f |S T y) (f ).Therefore, by definition of the conditional entropy and using the law of the unconscious statistician, it holds thatH(f | ỹ) = - log p(f | ỹ)p(f | ỹ)p( ỹ) df d ỹ = E p( ỹ) H p(f | ỹ) (f ) = E p(y) H p(f |S T y) (f ) = H p(f |S T y) (f )Therefore we can rewrite the mutual information in terms of prior and posterior entropy, such thatH(f ) -H f | S T y = H(f ) -H p(f |S T y) (f ) = 1 2 log det(K) + n log(2πe) -log det(K -KC i K) -n log(2πe) = - 1 2 log det(K -KS(S T KS) -1 S T K) det(K -1 )Via the matrix determinant lemmadet(A + U W V T ) = det(W -1 + V T A -1 U ) det(W ) det(A), we obtain = - 1 2 log det(-S T KS + SKK -1 K) det(-(S T KS) -1 ) det(K) det(K -1 ) = - 1 2 log det(S T KS -S T KS) det(-(S T KS) -1 ) = - 1 2 log det(σ 2 S T S(S T KS) -1 ) = 1 2 log det(σ -2 (S T S) -1 S T KS) = 1 2 log det((S T S) -1 S T KS) -i log(σ 2 ) = 1 2 log det(L -T S T KSL -1 ) -i log(σ 2 )for L a square root of S T S. Now we can upper bound the above as followsmax S∈R n×i H(f ) -H f | S T y ≤ max S∈R n×i 1 2 log det( ST K S) -i log(σ 2 ) = 1 2 log det(U KU T ) -i log(σ 2 )) = 1 2 i j=1 log(λ j ( K)) -i log(σ 2 )where U are the orthonormal eigenvectors of K for the largest i eigenvalues. Now choosing S = U achieves the upper bound since U T U = I and therefore S i = U is a solution to the optimization problem.Finally using the argument above and since H(f ) does not depend on S, we have thatarg max S∈R n×i H(f ) -H f | S T y = arg max S∈R n×i H(f ) -H p(f |S T y) (f ) = arg min S∈R n×i H p(f |S T y) (f ) .This proves the claim.Algorithm S1: CaGP = IterGP: Iterative formulation as in Wenger et al. [19] Input: GP prior GP(µ, K), training data (X, y) Output: (combined) GP posterior GP(µi, Ki)1 procedure CAGP(µ, K, X, y, C0 = 0) Time Space 2 while not STOPPINGCRITERION() do 3 si ← POLICY() Select action via policy. 4 ri-1 ← (y -µ) -Kvi-1 Residual. O(n 2 ) O(n) 5 αi ← s T i ri-1 Observation. O(k) O(1) 6 zi ← Ksi O(nk) O(n) 7 di ← Σi-1 Ksi = si -Ci-1zi Search direction. O(ni) O(n) 8 ηi ← s T i KΣi-1 Ksi = z T i di O(n) O(1) 9 Ci ← Ci-1 + 1 η i did T i Precision matrix approx. Ci ≈ K-1 . O(n) O(ni) 10 vi ← vi-1 + α i η i di Representer weights estimate. O(n) O(n) 11 Σi ← Σ0 -Ci Representer weights uncertainty. 12 µi(•) ← µ(•) + K(•, X)vi Approximate posterior mean. O(n⋄n) O(n⋄) 13 Ki(•, •) ← K(•, •) -K(•, X)CiK(X, •) Combined covariance function. O(n⋄ni) O(n 2 ⋄ ) 14 return GP(µi, Ki)Algorithm S2: CaGP: Batch VersionInput: GP prior GP(µ, K), training data (X, y) Output: (combined) GP posterior GP(µi, Ki) 1 procedure CAGP(µ, K, X, y) Time Space 2Si ← POLICY() Select batch of actions via policy.EQUATIONreturn GP(µi, Ki)We provide an open-source implementation of CaGP-Opt as part of GPyTorch. To install the package via pip, execute the following in the command line:pip install git+https://github.com/cornellius-gp/linear_operator.git@sparsity pip install git+https://github.com/cornellius-gp/gpytorch.git@computation-aware-gps-v2 pip install pykeopsTo better understand whether the overconfidence of SVGP at inducing points observed in the visualization in Figure 1 holds also in higher dimensions, we do the following experiment. For varying input dimension d ∈ {1, 2, . . . , 25}, we generate synthetic training data by sampling n = 500 inputs X uniformly at random with corresponding targets sampled from a zero-mean Gaussian process y ∼ GP(0, K σ ), whereK σ (•, •) = K(•, •) + σ 2 δ(•,•) is given by the sum of a Matérn(3/2) and a white noise kernel with noise scale σ. We optimize the kernel hyperparameters, variational parameters and inducing points (m = 64) jointly for 300 epochs using Adam with a linear learning rate scheduler. At convergence we measure the average distance between inducing points and the nearest datapoint measured in lengthscale units, i.e.dl (Z, X) = 1 m m i=1 min j ∥z i -x j ∥ diag(l -2 ) (S25)where l ∈ R d 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.EQUATION)The results of our experiments are shown in Figure S3 . We find that as expected the inducing points are optimized to lie closer to datapoints than points sampled uniformly at random. However, the inducing points lie increasingly far away from the training data as 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 = 4 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 . To better understand how the kernel hyperparameters of each method evolve during training, we show their trajectories in Figure S5 for each dataset. Note that we only show the first three lengthscales per dataset (rather than up to d = 26). 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 in 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 as motivation as long as it is clear that these goals are not attained by the paper.• 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 or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to 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, it 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 to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , 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 authors 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 be possible for other researchers to have some path to reproducing or verifying the results.• 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 (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 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 error rates). • If error bars are 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.Wenger et al. [19] named their computation-aware inference algorithm "IterGP", to emphasize its iterative nature (see Algorithm S1). We adopt the naming "CaGP" instead, since if the actions S are not chosen sequentially, it is more efficient to compute the computation-aware posterior non-iteratively (see Algorithm S2).At first glance, SVGP satisfies a similar result (Theorem S3). However, this statement does not express the variance in terms of the worst-case squared error to the "true" latent function. See Section S1.2 for details.Given this observation, we sometimes abuse terminology and refer to the actions as "projecting" the data to a lower-dimensional space, although Si does not need to have orthonormal columns.Wenger et al. [Sec. 2.1 of 19] showed that CaGP using CG residuals as actions recovers Conjugate-Gradient-based GPs (CGGP)[7,9,10] in its posterior mean, extending this observation to CGGP.For a detailed analysis see Algorithms S1 and S2 in the supplementary material, which contain line-by-line time complexity and memory analyses.While SVGP is arguably linear time since it will eventually loop through all training data, each computation of the ELBO uses a constant time stochastic approximation. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: I have a couple of concerns regarding the models used in the experiments:

The 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.

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.
questions: 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".

"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. 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: limitations: The authors include a well written section on the limitations of IterGP-Opt.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- 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.
- 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.
questions: 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.
- “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. 
- 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.
- 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.
- 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?
- Should I interpret equation 12 as stating the $\mathbf{S}$ is block diagonal? If so, is this said explicitly somewhere?
- Could you explain in more detail why SGPR diverges in these experiments? Have you checked if this is resolved by double precision?


### Other

- Line 49, missing space after comma before “a”
- 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?
- Line 73, h is introduced here as the latent function and never appears again. I’m a bit unclear why it is needed.
- 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: limitations: The authors point out limitations well; most importantly, the restriction to Gaussian likelihoods in the present work.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

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.

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!
questions: questions: In order of relative importance:
1. 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!
2. 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?
3. 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.
4. Out of interest, since the ELBO for IterGP-Opt is not stochastic, could L-BFGS-B be used there to speed up optimization?
5. 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?

Below I'll list some minor questions and typos:
- 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!
- Missing space in line 49 in "introduced,a"
- 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
- line 99: "we optimizer" -> "we optimize"
limitations: 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.

-----
Post-discussion update
-----
I 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.

I 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.

> 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.

Thanks 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.

> 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). 

Yes, 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.

> 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.

I'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.

Finally, 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:
> 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

Where 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.

Comments: title: Response to rebuttal
Answers to Questions [continued]
Answers to Questions
Answers to Questions [continued]
Answers to Questions
comment: Thank you for your suggestions and continuous engagement during the discussion phase!
Thank you for providing the additional experiments, I feel it offers important context, as the reader can now gauge the impact of the contribution (i.e. compared to using sub-optimal hyperparameters for IterGP-CG). 
I'm happy to update my score accordingly, and looking forward to the discussion phase.
Thank you for carefully reading our responses and of course for your detailed and actionable feedback!

> I am somewhat confused by how the number of epochs is chosen for Adam in the updated results. [...] How did the authors determine when to stop Adam optimization (looking at the paper I couldn't find where this was mentioned, but apologies if I've missed it)?

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.

> At the end of the day, the results should be the same (theoretically speaking, of course) regardless of how you choose the optimizer, with the only difference really being time taken and memory allocation, [...]

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)](http://gaussianprocess.org/gpml/chapters/RW.pdf)). Which modes an optimizer converges to depends on the initialization, the choice of optimizer and the learning rate (schedule / line search). 

> [...] but this is not reflected in these updated results. This is especially striking in Parkinsons and Bike for both SGPR and IterGP-Opt, which do not seem to reflect the authors' recommendation for Adam over L-BFGS (as stated in a response to reviewer 7YnY).

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 $d+2 + n$ parameters. As, for example, the [PyTorch documentation on L-BFGS](https://pytorch.org/docs/stable/generated/torch.optim.LBFGS.html) 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`.

---

> I appreciate that the authors have agreed to include an analysis of the hyperparameter values of CholeskyGP/exact GPR in their manuscript, but I still think that explicitly considering the ELBOs/LML estimates could be somewhat instructive, as indicated in my original review. Firstly, because solely considering predictive metrics can be misleading when the model is mis-specified, as it most likely is here on these UCI datasets, but also in that it might shed light on some of the discussions being had about optimization.

We apologize for overlooking this suggestion for improvement in your original response! We logged training losses (i.e. ELBOs / (approximate) LMLs) for each epoch in all our runs and will gladly add these to the manuscript (e.g. as an additional column in Table 1, and as an additional row in Figure 1). 

---

We hope our additional answers were helpful in forming your final opinion about our work! If there are any outstanding questions, please do not hesitate to voice them.
Thank you for your rebuttal and detailed response to my review and others. I can confirm that I have carefully read the rebuttals as well as the discussion (up to this point). In general, I think the updated results represent a significant improvement to the quality of the paper. I would also like to thank the authors for correcting some of the confusion I had regarding their paper as well as GPyTorch. However, I generally like to see how the discussion period plays out, and especially if the other reviewers are satisfied with the authors' responses, before committing to change my score one way or another.

A couple of minor points still come to mind:
- I am somewhat confused by how the number of epochs is chosen for Adam in the updated results. At the end of the day, the results should be the same (theoretically speaking, of course) regardless of how you choose the optimizer, with the only difference really being time taken and memory allocation, but this is not reflected in these updated results. This is especially striking in Parkinsons and Bike for both SGPR and IterGP-Opt, which do not seem to reflect the authors' recommendation for Adam over L-BFGS (as stated in a response to reviewer 7YnY). How did the authors determine when to stop Adam optimization (looking at the paper I couldn't find where this was mentioned, but apologies if I've missed it)?
- I appreciate that the authors have agreed to include an analysis of the hyperparameter values of CholeskyGP/exact GPR in their manuscript, but I still think that explicitly considering the ELBOs/LML estimates could be somewhat instructive, as indicated in my original review. Firstly, because solely considering predictive metrics can be misleading when the model is mis-specified, as it most likely is here on these UCI datasets, but also in that it might shed light on some of the discussions being had about optimization.

Overall though, while I continue to reserve judgement about the final score until later, I do think that in my view this paper has already been significantly improved by the responses the authors have made.
Based on your request, we have added the following **experiment comparing CGGP trained via the procedure given in [Wenger et al. (2022b)](https://arxiv.org/abs/2107.00243) and our proposed model selection procedure for IterGP-CG**:

We trained CGGP on all small to medium-sized datasets via L-BFGS for 100 epochs in double precision with a sweep over the initial learning rate. Due to the very limited time until the end of the rebuttal period we ran this experiment only for one seed, we will add the remaining runs for the camera-ready version.
As the table below shows, *IterGP-CG trained via our proposed objective outperforms or matches CGGP on all datasets (i.e difference in performance larger than one standard deviation), except in terms of RMSE on Protein.* 
The significant lower runtime on KEGGu is explained by the fact that after epoch 9 CGGP ceases to improve.

| Dataset    | Method    | Optim. | LR  | Epoch | Test NLL | Test RMSE | Avg. Runtime |
|------------|-----------|--------|-----|-------|----------|-----------|--------------|
| Parkinsons | CGGP      | L-BFGS | 1.0 | 68    | -2.734   | 0.011     | 1min 23s     |
|            | IterGP-CG | Adam   | 1.0 | 250   | -2.936   | 0.009     | 1min 44s     |
| Bike       | CGGP      | L-BFGS | 1.0 | 15    | -2.053   | 0.021     | 2min 14s     |
|            | IterGP-CG | Adam   | 1.0 | 250   | -2.042   | 0.024     | 5min 17s     |
| Protein    | CGGP      | L-BFGS | 0.1 | 28    | 0.852    | 0.511     | 16min 13s    |
|            | IterGP-CG | Adam   | 1.0 | 27    | 0.820    | 0.542     | 1min 26s     |
| KEGGu      | CGGP      | L-BFGS | 1.0 | 9     | -0.510   | 0.127     | 7min 29s     |
|            | IterGP-CG | Adam   | 1.0 | 229   | -0.699   | 0.120     | 39min 5s     |


Now since CGGP and IterGP-CG share the same posterior mean (see Section 2.1 of [Wenger et al (2022)](http://arxiv.org/abs/2205.15449)), the Test RMSE for IterGP-CG is guaranteed to be identical when using the same hyperparameters. Empirically as Wenger et al (2022) show in Figure 4, the Test NLL is also the same. Therefore the reported numbers here will be the same if computing the approximate posteriors with IterGP-CG for the hyperparameters given by CGGP as done in Wenger et al. (2022). (Note that the reported 2x speedup of IterGP-CG vs CGGP by Wenger et al. 2022 only applies to inference, not the CGGP training procedure).

Thank you for your suggestion to improve our paper. We would kindly ask you to take this newly added experiment into account in your final evaluation.
The authors have adressed the concerns I raised in my review. I maintain my initial score and assessment that the paper should be accepted.
The hyperparameters were chosen with the training procedure for CGGP as proposed in Wenger et al. [29] to favor CGGP. They do the same in the follow-up experiment in Figure 5 where they choose the hyperparameters given by SVGP for their variant IterGP-PI.

Would you be satisfied if we added CGGP training runs and then reported the test performance for the approximate posterior computed by IterGP-CG (as done by Wenger et al 2022)? This would be possible.
Yes I understand the significance of Figure 4, that during optimisation the hyperparameters were fixed, but there is still a procedure on how to obtain those fixed hyperparameters as cited in the paper:
> For all datasets, we select hyperparameters using the training procedure of Wenger et al. [29].

The point is, it is important to demonstrate explicitly what improvement is brought to IterGP-CG compared to what was attainable with IterGP-CG before.
> The CholeskyGP model is chosen to be light grey, which when bolded barely stands out. Please adjust accordingly, simply making it black should suffice (or alternatively, don't set it to bold at all and only highlight the best of the approximate methods?).

We are happy to make the CholeskyGP stand out more.

> There seem to be learning rates given for the new BFGS-optimised methods, is this just a typo?

This is not a typo. As we write in the response above this is the *initial* learning rate, which then gets updated by the line search. See also the [corresponding lines](https://github.com/pytorch/pytorch/blob/26b0a0c2f37a8ad376f261df7bb4fee65ff2f230/torch/optim/lbfgs.py#L419-L423) in the PyTorch implementation of L-BFGS and its [`lr` argument](https://pytorch.org/docs/stable/generated/torch.optim.LBFGS.html).

> Regarding the presentation of the results for IterGP-CG: in order to make this contribution much clearer, might it be helpful to explicitly compare these results against the previous training procedure for IterGP-CG, as was used in e.g. Figure 4 of arXiv:2205.15449. This would aid the reader to gauge the significance of the contribution.

We believe there might be a misunderstanding here as to what the contribution of Wenger et al. (2022) is. In Figure 4 they show the loss as a function of the solver iterations when computing an *increasingly better approximate posterior with fixed hyperparameters*. We propose a *model selection / hyperparameter optimization procedure* for this approximate posterior computed with a fixed number of solver iterations (i=512). Wenger et al. (2022) *do not consider model selection / propose a training procedure*. It is therefore not possible to compare their experiment to ours, since the two papers consider two different problems.

Does this clarify our contribution? We are happy to give more detail if that should alleviate your concerns.
Thank you for responding to my concerns, and updating the experimental table accordingly. A few suggestions regarding the table: 

- The CholeskyGP model is chosen to be light grey, which when bolded barely stands out. Please adjust accordingly, simply making it black should suffice (or alternatively, don't set it to bold at all and only highlight the best of the approximate methods?).

- There seem to be learning rates given for the new BFGS-optimised methods, is this just a typo?

>  In particular, we demonstrate in this work how to perform model selection for IterGP-CG. 
- Regarding the presentation of the results for IterGP-CG: in order to make this contribution much clearer, might it be helpful to explicitly compare these results against the previous training procedure for IterGP-CG, as was used in e.g. Figure 4 of arXiv:2205.15449. This would aid the reader to gauge the significance of the contribution.
We would like to issue a brief correction here to one equation in our response above. The second equality in the definition of the information-optimal actions should be the entropy of the computation-aware posterior, rather than the KL divergence to the prior:

$$
S_i = \dots = \textrm{argmin}_{S \in \mathbb{R}^{n \times i}} \mathbb{E}\_{p(f(X) \mid S^\top y)}[- \log p(f(X) \mid S^\top y)]
$$
> 5. 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?

We intended to reference the way the approximate posterior is defined (i.e. Equation (5) and line 105) rather than the equation above line 111 in our statement. Equation (5) and line 105 imply a monotonic decrease in the marginal variance as a function of $i$. Observe that the precision matrix approximation $C_i$ (defined in line 105) has rank $i$, assuming the actions, i.e. columns of $S_i$, are linearly independent and non-zero (as assumed in Wenger et al. (2022)). Therefore we can rewrite 
$C_i = \sum_{\ell=1}^i \tilde{d}\_\ell \tilde{d}_\ell^\top $ as a sum of rank 1 matrices. Note that this is actually how Wenger et al. (2022) compute $C_i$ in Algorithm 1 in their paper. Therefore the marginal variance at iteration $i$ for $i \leq j$ is given by 

$$K_{i}(x, x) = K(x, x) - \sum_{\ell=1}^i K(x, X) \tilde{d}\_\ell \tilde{d}\_\ell^\top K(X, x)=K(x, x) - \sum_{\ell=1}^i (K(x, X) \tilde{d}\_\ell)^2 \geq K(x, x) - \sum_{\ell=1}^j (K(x, X) \tilde{d}\_\ell)^2= K_{j}(x,x)$$

Now for $i=n$, since the actions are assumed to be linearly independent and non-zero (see Wenger et al. (2022)), $S_i$ has rank $n$ and thus $C_i = S_i(S_i^\top \hat{K}S_i)^{-1}S_i^\top = \hat{K}^{-1}$. Therefore we have $K_{i}(x, x) \geq K_{j}(x, x) \geq K_{\star}(x, x)$.

> 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!

We apologize for not taking this into account. We have swapped out the colormap for [Seaborn's colorblind palette `sns.color_palette("colorblind", 5)`](https://seaborn.pydata.org/generated/seaborn.color_palette.html).
### Questions

> 1. 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!

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 to compute solves with $K(Z, Z)$ with an iteratively increasing jitter should the decomposition fail. In other words, we did not rely on iterative methods for accelerated solves for SGPR and SVGP. We will make this more explicit in the description of the experimental setup in the final version of the paper.

> 2. 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?

The asymptotic memory complexity $O(nm)$ is indeed the same for $m=i$, but since we are using $m=1024$ inducing points versus $i=512$ for IterGP-Opt the memory cost of SGPR is two times higher, ignoring any constant factors arising from implementation differences.
When training SGPR in double precision (with LBFGS), the memory required increases by another factor of two, giving an overall memory requirement that is at least four times higher. These differences lead to prohibitive memory consumption when we attempt to train SGPR on Road on the GPU we used in our experiments.

> 3. 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.

We optimized the hyperparameters and inducing points of SVGP in Figure 1 using Adam with an initial learning rate of 0.1 for 200 epochs and a linear learning rate decay. Your experience matches ours that in one dimension the depicted scenario does not occur frequently -- it depends on the initialization and choice of optimization hyperparameters. However, we believe *it faithfully reflects SVGP's behavior in higher dimensions*. 

To substantiate this claim, we ran an experiment on synthetic data based on a latent function drawn from a GP with increasing dimension of the input space. We find that the average Mahanalobis distance (induced by the lengthscales) across inducing points to the nearest datapoint increases with dimension (see Figure R2 of author response PDF). We compare this to randomly placed inducing points. While SVGP does move inducing points closer to data points on average, they are still becoming increasingly distant from the data measured in lengthscale units. Further, we also provide more evidence for our claim that SVGP's latent variance often is very small (at inducing points) relative to the total variance, meaning almost all deviation from the posterior mean is considered to be observational noise. The right plot in Figure R2 of the attached rebuttal PDF shows this across dimensions. For $d\geq 4$ approximately 95% of the predictive variance is due to observational noise.

> 4. Out of interest, since the ELBO for IterGP-Opt is not stochastic, could L-BFGS-B be used there to speed up optimization?

Thank you for this suggestion! We ran an additional set of experiments on all but the two largest datasets for IterGP-Opt. We trained using L-BFGS in double precision with a Wolfe line search run for a maximum of 100 epochs and averaged across five seeds. As the results in Table R1 of the attached rebuttal PDF show, performance does increase slightly at comparable runtime (only for Bike the runtime increased considerably).
> - 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.

In an earlier answer, we demonstrated that the marginal variance decreases monotonically as a function of $i$ (albeit not necessarily strictly), given linearly independent, non-zero actions. We immediately obtain that the marginal precision increases monotonically up to the marginal precision of the posterior as the computational budget increases, i.e. $i \to 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?

One way to understand how different instances of IterGP differ based on different choices of actions, i.e. the columns of the matrix $S_i$, is by interpreting the magnitude of the action entries as how the computational budget is distributed across the $n$ observations (see also Figure 2 of Wenger et al. (2022)). This interpretation arises from the formulation of IterGP conditioning on $S_i^\top y$ (see also lines 122-125). For example, unit vector actions correspond to spending all budget on a single data point, sparse actions target only a subset of data points, and CG and eigenvector actions distribute the budget across all data points with different weights. The top row in Figure 2 corresponds to eigenvector actions, hence we plot the magnitude of the eigen*vector* entries. The rows below illustrate the two variants of IterGP which we consider in our paper. In the plot on the right in Figure 2, we show the Grassman distance, which is the Euclidean norm of the vector of principal angles between the two subspaces. We've added its definition to the paper.

> - Should I interpret equation 12 as stating the $S$ is block diagonal? If so, is this said explicitly somewhere?

The matrix $S$ consists of $i$ blocks (*), which are of size $k \times 1$. We state this implicitly in lines 183-184, but we've made this more explicit in the final draft.

> - Could you explain in more detail why SGPR diverges in these experiments? Have you checked if this is resolved by double precision?

For SGPR, we computed Cholesky decompositions in single precision with an iteratively increasing additive jitter should the decomposition fail as per Section 4.1 of Ober et al (2022).
The iteratively added jitter could explain the observed divergent behavior instead of an outright failure in single precision.
During the rebuttal, we added results for SGPR trained in double precision with LBFGS to our experiments, which does not show any divergence during training as Figure R1 in the attached rebuttal PDF shows. Training in this way improves the performance of SGPR over SVGP, but IterGP-Opt still matches or outperforms SGPR, except on "Bike".

> - 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?

The claim is that they are conservative. We've corrected this in the final version.
### 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.

We would like to clarify that we claim this in reference to Figure 1 as an illustration of a more general phenomenon that can occur when using SVGP. Specifically, we claim that "SVGP, [...] is overconfident at the locations of the inducing points if they are not in close proximity to training data, which becomes increasingly likely in higher dimensions. This phenomenon can be seen in a toy example in Figure 1 (middle left), where SVGP has no latent uncertainty where data is absent (lower row, cf. other methods)." 
We've modified this claim to "has near zero latent uncertainty at the inducing point away from the data".
Additionally in the caption of Figure 1 we write referring to the same figure, that "SVGP expresses almost no posterior variance in the middle region, and thus almost all deviation from the posterior mean is considered to be observational noise."

To back up these two claims we added a new experiment, which measures the distance of inducing points to the data (measured in lengthscale units) as a function of input dimension and the ratio between posterior and predictive variance at the inducing points (see Figure R2 of rebuttal PDF). We find that inducing points are optimized to lie increasingly farther away from the data as the dimension of the input space increases. Further, we observe that in higher dimensions ($d\geq 4$) the latent variance can be as little as 5% of the predictive variance.
Note that all our statements and observations are made under the assumption of a fixed number of data points $n > m$ and inducing points $m$ and therefore do not conflict with statements about the convergence of SVGP to the exact posterior. 

> - “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.

We refer to the posterior variance updates and specifically *the form of the downdate*. We've clarified this in the draft. To see why this guarantees the variances are always overestimated (they monotonically decrease as a function of $i$), consider the following. Observe that the precision matrix approximation $C_i$ (defined in line 105) has rank $i$, assuming the actions, i.e. columns of $S_i$, are linearly independent and non-zero (as assumed in Wenger et al. (2022)). Therefore we can rewrite 
$C_i = \sum_{\ell=1}^i \tilde{d}\_\ell \tilde{d}_\ell^\top $ as a sum of rank 1 matrices. Note that this is actually how Wenger et al. (2022) compute $C_i$ in Algorithm 1 in their paper. Therefore the marginal variance at iteration $i$ for $i \leq j$ is given by 

$$K_{i}(x, x) = K(x, x) - \sum_{\ell=1}^i K(x, X) \tilde{d}\_\ell \tilde{d}\_\ell^\top K(X, x)=K(x, x) - \sum_{\ell=1}^i (K(x, X) \tilde{d}\_\ell)^2 \geq K(x, x) - \sum_{\ell=1}^j (K(x, X) \tilde{d}\_\ell)^2= K_{j}(x,x)$$
Now for $i=n$, since the actions are assumed to be linearly independent and non-zero, $S_i$ has rank $n$ and thus $C_i = S_i(S_i^\top \hat{K}S_i)^{-1}S_i^\top = \hat{K}^{-1}$. Therefore we have $K_{i}(x, x) \geq K_{j}(x, x) \geq K_{\star}(x, x)$.