Posterior and Computational Uncertainty in Gaussian Processes
Authors: Jonathan Wenger
Geoff Pleiss
Marvin Pförtner
Philipp Hennig
John Patrick Cunningham
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian processes
computational uncertainty
numerical methods
probabilistic numerics
probabilistic linear solvers

Abstract: Gaussian processes scale prohibitively with the size of the dataset. In response, many approximation methods have been developed, which inevitably introduce approximation error. This additional source of uncertainty, due to limited computation, is entirely ignored when using the approximate posterior. Therefore in practice, GP models are often as much about the approximation method as they are about the data. Here, we develop a new class of methods that provides consistent estimation of the combined uncertainty arising from both the finite number of data observed and the finite amount of computation expended. The most common GP approximations map to an instance in this class, such as methods based on the Cholesky factorization, conjugate gradients, and inducing points. For any method in this class, we prove (i) convergence of its posterior mean in the associated RKHS, (ii) decomposability of its combined posterior covariance into mathematical and computational covariances, and (iii) that the combined variance is a tight worst-case bound for the squared error between the method's posterior mean and the latent function. Finally, we empirically demonstrate the consequences of ignoring computational uncertainty and show how implicitly modeling it improves generalization performance on benchmark datasets.

Posterior and Computational Uncertainty in Gaussian Processes
Jonathan Wenger
Geoff Pleiss?
Marvin Pförtner
Philipp Hennig
John P. Cunningham
University of Tübingen
Columbia University
Max Planck Institute for Intelligent Systems, Tübingen
Abstract
Gaussian processes scale prohibitively with the size of the dataset. In response, many approximation methods have been developed, which inevitably introduce approximation error. This additional source of uncertainty, due to limited com- putation, entirely ignored when using the approximate posterior. Therefore in practice, GP models are often much about the approximation method as they are about the data. Here, we develop a new class of methods that provides consistent estimation of the combined uncertainty arising from both the finite number of data observed and the finite amount computation expended. The most common GP approximations map to an instance in this class, such methods based on the Cholesky factorization, conjugate gradients, and inducing points. For any method in this class, we prove (i) convergence of its posterior mean in the associated RKHS, (ii dccomposability of its combined posterior covariance into mathematical and computational covariances, and (iii) that the combined variance is tight worst-case bound for the squared error between the method's posterior mean and the latent function. Finally, empirically demonstrate the consequences of ignoring compu- tational uncertainty and show how implicitly modeling it improves generalization performance benchmark datasets.
Introduction
Gaussian processes (GPs) are an expressive probabilistic model class, but their prohibitive scaling necessitates approximation [1]. A range of approximations based kernel [2-10] or precision matrix [11-14] estimates, inducing point methods [15-22), and iterative solvers [23-29] have been proposed. These methods all mathematical posterior, which exists theoretically but cannot be accessed given limited computational resources. The approximate posterior is then used as a direct replacement of the mathematical posterior in downstream applications. Doing so, however, completely ignores the fact that we only expended limited amount of compute. By analogy the typical GP operation, where limited data induces modeling етrоr captured by mathematical uncertainty, our work is motivated by the fact that limited computation induces approximation error that must be captured by computational uncertainty.
affordable amount of computation obtain approximation of the
Here, we introduce IterGP, a class of methods which return a combined uncertainty that is the sum of mathematical and computational uncertainty. Figure 1 illustrates the difference between ignoring computational uncertainty and explicitly modeling it. We perform GP regression using Matérn ) kernel on toy dataset and compare SVGP [22] to its analog in framework, IterGP-PI O) for fixed set of inducing points. The computational shortcuts of inducing point methods can lead unavoidable biases in their posterior mean and covariance [30, 31]. As Figure illustrates, SVGP
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

SVGP
HeerGi.PP (ours)
Aggrss, Appros, Unominty
Compuratoral Uncensiity Uscarderr
Uncertainty
Inducing
Figure 1: Modeling computational uncertainty improves GP approximation.
may underestimate the marginal variance where inducing points do not coincide with datapoints. In contrast, IterGP is guaranteed to overestimate the mathematical uncertainty with the difference precisely given by the computational uncertainty (O). Additionally, the computational uncertainty is a worst-case bound (-) on the error the approximate posterior mean.
To be clear, this overestimation is desirable: IterGP is not typical approximation in the sense that its combined posterior attempts to approximate the mathematical posterior. Rather, IterGP recognizes that the mathematical posterior exists, but we do not have access to it, given computational constraints. Finite compute is as true a source of posterior uncertainty as finite data. Taking this view seriously, the true goal of GPs in the limited compute regime should in fact be to track combined uncertainty. This intuition motivates IterGP and is formally a feature of our results. We show that, you update your GP via computation, specifically matrix-vector multiplication, then the combined uncertainty of the IterGP algorithm is precisely the correct object to capture your belief (Theorem 2) - in the same way the mathematical posterior is the correct object given finite data and unlimited computation.
Formally, IterGP probabilistic numerical method [32-35]. It treats the (unknown) representer weights as a latent variable with prior belief that, when marginalized out, corresponds to GP prior conditioned on no data. We then use computational primitive (matrix-vector multiplication) that corresoonds to tractable Bayesian updates on the representer weight distribution; that is, conditioning on computations performed on the data. The resulting belief can then be marginalized out to obtain a closed form, tractable expression for the combined mathematical plus computational uncertainty. This uncertainty quantification can done exactly in quadratic time and linear space complexity.
Our framework admits three key theoretical properties. First, common GP approximations such as the partial Cholesky, the method of conjugate gradients and inducing point methods (e.g. SVGP) map to corresponding IterGP instance. Therefore, these approaches can either be directly extended or modified to properly account for computational uncertainty. Second, the approximate posterior mean of any method in our proposed class converges to the mathematical posterior mean in RKHS norm in most n steps, where the convergence rate determined by the choice method (Theorem 1). Third, the combined uncertainty is tight worst case bound on the relative error between the approximate posterior mean and the latent function (Theorem 2). To the best of our knowledge no existing GP approximation has this last property; analoguous guarantee only holds for exact GPs [36, Sec. 3.4].
Contribution This work introduces IterGP, which defines new class GP approximations that accounts for computational uncertainty arising from limited computation. Some IterGP instances extend classic methods with improved uncertainty quantification (Table 1). For any method this class, prove that the approximate posterior mean converges to the mathematical posterior mean (Theorem 1) and that the combined uncertainty is a tight worst-case bound on the relative distance to the latent function one is trying to learn (Theorem 2, Corollary 1). We demonstrate empirically that modeling computational uncertainty can either save computation or improve generalization on a set show that it is possible exactly quantify the
of regression benchmark datasets. In conclusion, inevitable error of GP approximations at quadratic cost by propagating said error to the posterior in the form of computational uncertainty.

---

2
Computation-Aware Gaussian Process Inference
We aim to learn latent function h: X y from X' S F. to y C R given training dataset X (ZI,..., In) E Randd of n. inputs z; E Rd and corresponding outputs y = (...) E IR..
Gaussian Processes A stochastic process f GP(H, k) with mean function μ.: : R - R and kernel : Rd x F- + is called Gaussian process (GP) if any collection of function values f =
~ (N, K) is jointly Gaussian with μi; u(æ;)) and K = k( 2z). Assuming observation noise | f ~~ [ff, o the posterior distribution at test inputs X, is given by fo ~ (N,,((X)) k. (X:: X))) where the posterior mean and covariance functions given by
(-(+) = u(-) + k(
(y ()) and k..,)) = - X) -() ((x,))))))))
(I)
where K = K + 1 E I=x, Computing the representer weights v. = K - 14) exactly (as well as the posterior variance) is prohibitive given our limited computational budget.
Learning Representer Weights Consider the conditional distribution of the latent GP given its representer weights:
p(fo v.) (N((X))) + ,(x,, X)e,, (X.,
When v, is known exactly, we recover eq. (1). However, if we instead treat v. as a random variable with the prior p(u.) (e.,; 0. -- then the resulting marginal J p(f,, | v. v.)p(e,)) dv, the GP prior N ((((X),, k(X., X.)). Our goal is to update this prior by iteratively applying the tractable computational primitive (i.e. matrix-vector multiplies). More specifically, each iteration conditions the current belief distribution p(u.) (Nv,; 1; E;- on one-dimensional projec- tion of the current residual r = (y) 11) KV where the projection is defined by an arbitrary vector s,:
dv, recovers
a, == s = s! ((y p kv, = k(, - -v,-1)
(3)
The choice of actions si, which intuitively weight the approximation error of selected datapoints, defines different instances of our IterGP framework. Computing eq. (3) requires single matrix vector multiplication. After computing Os, we can perform an exact Bayesian update of p(v.) via linear Gaussian identities. The updated p(v.) (conditioned on a;;) is N(v. v., E,), with
8 U; + E K (8. KE Ks,)
K K(u. U.-1) C,(y 11)
(4)
d;
E E, E (s K KE Ks,
KE
K-1 C..
(5)
=d,
where C, -- Eド
1 d,d. = (.S ST a rank-i matrix (see Proposition S3 for details).
nj With each computation, the uncertainty about the representer weights contracts C, - K Eg as + n. After n iterations, Cn K meaning we fully recovered the representer weights with zero uncertainty. The consistent estimate for the representer weights is consequently u, Cy 14).
Combining Mathematical and Computational Uncertainty We now have a belief p(v.) = Ve,,;,, Σ,) about the representer weights reflecting the expended computation. To account for this computational uncertainty, we treat the representer weights as a latent variable of the mathematical posterior by reparameterizing p(f. y) =f(fo v.) and then marginalizing. The resulting marginal considers all possible representer weights which would have resulted the same computational observations and therefore implicitly adds the uncertainty coming from the computation itself. Since the posterior mean of a GP is a linear function of the representer weights, the marginal distribution is given by p(fo)
pff. v. v.)p(v.) du. (fo: us (X.). Ks (X., X.)), where
(.)()) = p(-) + k., X)u; k.(,,.) k(,,) kC,, X k(X, ) + K(, X)E,,(X
XXCC (6)
mathematical uncertainty @
computational uncertainty
combined uncertainty
since E, = K
C,,' -C.'' As we perform more computation, the computational uncertainty reduces and we approach the mathematical uncertainty. We note that, while the individual terms are com- putationally prohibitive, the combined uncertainty can be evaluated cheaply since the approximate
'While we derive the combined posterior from probabilistic numeries perspective, we can alternatively interpret eq. (6) the GP prior conditioned on linearly transformed data S.y N o S.).)
3

---

Combined Uncertainty
Mathematical Uncertiitty
Computational Uncertuinty
Figure 2: Decomposition of the combined uncertainty. The combined uncertainty output by IterGP decomposes into the mathematical uncertainty and computational uncertainty (0). After i == iterations Algorithm computational uncertainty is small in parts of the input space where there either is data or computation was "targeted" (0). Which datapoints are targeted in each iteration and to what degree is defined by the magnitude of the action vector elements (si))- Different instances of IterGP either reduce computational uncertainty locally (e.g. ttrGP--Coll IterGP-PP PPPP or globally (e.g. IterGP-CG). After n iterations the mathematical uncertainty is recovered.
precision matrix C, of low rank. Figure 2 illustrates that computational uncertainty large where there are data and we have not targeted computation yet. Different methods from our proposed class target computation in different parts of the input space. Where there is no data the prior good approximation of the posterior and therefore computational uncertainty is low.
Algorithm 1 computes an estimate of the representer weights v, and the rank-i precision matrix approximation C. A specific instance of IterGP is defined by a sequence of actions 8,. To gain an intuition for how Algorithm operates, it helps to interpret it as targeting given computational budget towards certain datapoints defined by Near datapoints z that are not targeted, i.e. (s:), = O, computational uncertainty remains unchanged. In fact, datapoints (I,, Myy that are never targeted up to iteration i are not needed to compute GP(us, k;), meaning that Algorithm I is inherently online and we can observe data sequentially without having to restart the algorithm (see Theorem S7).
Algorithm 1: A Class of Computation-Awar Iterative Methods for GP Approximation
Input: prior mean function μ., prior covariance function/ kernel k training inputs X, labels y Output: (combined) GP posterior <P(j,kk)
procedure ITERGP(HA k X.y)
2 (Nu, ku) < (H,
Initialize mean and covariance function with prior. Prior predictive mean.
3 J4 + (XX)
K + k(x,X) I21
Prior predictive kernel matrix. Stopping criterion.
5
while not STOPPINCRITERRIONO do
#; + POLICY()
ts Select action via policy (see Table for examples). Predictive residual.
7
Ti-1 + (y 11) Kv-1
8
ai + s!ri-1
Observation via information operator. Search direction.
9
d, - E., = -1- C
10
H < 8 K
Kdd.
Normalization constant.
C. + C.-1 dd,d
Precision matrix approximation C, x R-1
12
Q + Qi-1 kd.d K
Kernel matrix approximation Q: R K.
13
+ 0: Si di
Representer weights estimate.
n
14
Σ; Eo C:
Computational representer weights uncertainty.
15
P(v.) < N(v.; Σ,)
Belief about representer weights. > Approximate posterior function. Combined uncertainty.
16
(i()) + u() + k(; jo
17
k6, - k(,, X),(()
return (P(,,,
Greyed out quantities not needed to compute the combined posterior and are only included for clarity of exposition.

---

Table 1: Instances of Algorithm which map to commonty used GP approximations.
Method
Actions Si
Classic Analog
Reference
IterGP-Chol
E, ev, ( PCG P-1r; x.)
(partial) Cholesky
Theorem S3
IterGP-PBR
(partial) EVD SVD
Theorem S4
IterGP-CG IterGP-PI
(preconditioned) CG
Theorem S5 and Corollary S2
Nyström (SoR, DTC), SVGP
Section 2.1 and Theorem
2.1 Connection to Other GP Approximation Methods
IterGP extends the most commonly used GP approximations include computational uncertainty, with at most quadratic cost (see Table for summary and Figure 2, Figure S3 for illustration).
Cholesky Decomposition The (partial) Cholesky decomposition iteratively chooses datapoints or pivots 2; based on a given ordering. The resulting Cholesky factor is lower triangular and increases in rank each iteration, and well-chosen ordering achieves fast convergence (ef. [37, Thm. 2.3]). If one chooses standard unit vectors e; as actions corresponding to the selected datapoint per iteration, then Algorithm 1 recovers the partial Cholesky factorization exactly (Theorem S3).
Conjugate Gradients CG [38] with preconditioning for GP inference has become increasingly popular [24-29, 39, 40J. Algorithm recovers preconditioned CG exactly, if we choose either preconditioned conjugate gradients or residuals as actions (see Theorem S5 and Corollary S2). In fact, Algorithm I can even construct its own Siagoonal-plu-lowwww a few iterations with an arbitrary policy and then using the byproducts of these iterations for the preconditioner. For example, if we run IterGP-Chol initially, can construct an incomplete Cholesky preconditioner for subsequent CG iterations.
preconditioner by first running
Inducing Point Methods Inducing point methods, such as variants of the Nyström approximation [16), i.e. subset of regressors (SoR) [15, 41] and deterministic training conditional (DTC) [18, 42], as well as SVGP [22] share posterior mean, which by Theorem S6 takes the form
(SVGP()) = a(- x XxxxxxXxxxXXXx (7(
x) +
1) 2
Kzx -4)
(7)
where Z Raxi is set of inducing points and g(:: k( Z) 22 k( ZZ -). These approximations
also have very closely related posterior covariance functions [20, 43]. If we choose actions &, =
k(X, 2;), by Proposition S3, Algorithm returns posterior mean given by
p((-) = kc., X) K Zx (k(X, X X) + o² 1 KZx - pt).
(8)
Gram matrix K KS:
Choosing such actions, given by kernel functions k(, z;) centered at inducing points 2;, reduces computational uncertainty in regions close to inducing points (see IterGP-PI in Figure 2). where closeness is determined by the kernel. Comparing SVGP's and IterGP-PI's posterior mean provides probabilistic numerical perspective on why even for small KL-divergence between the approximating distribution of SVGP and the true posterior, the mean estimate can be far from the true mean [31, Prop. 3. As outlined Section 2, eq. (8) is a Bayesian update on the initially unknown representer weights v. K 14). The Gram matrix in eq. (8) describes how surprising the computational observations Kzx (y (1) ST v1) ST K v. of the representer weights should be, given the prior uncertainty Eu about them. SVGP uses similar form for the posterior mean (c.f. (7) and (8)), but the Gram matrix is "smaller" since (X, 3 k(X,XX X). This can be interpreted as inducing point methods being overconfident in their update of the representer weight estimates to achieve linear time complexity. As the inducing points approach the data points the two posterior mean functions jSvvr and μ become closer and are equivalent if the inducing points equal the training data.
2.2 The Cost of Computational Uncertainty
Quantifying combined uncertainty has greater cost than linear time GP approximations such as inducing point methods, due to its use of metrix-vector multiplication as the computational operation condition on the data. Algorithm in its most general form performs three matrix-vector products per iteration resulting in a quadratic time complexity (n?i) overall for i iterations. In this sense,
5

---

Algorithm represents middle ground between the mathematical posterior-whie incurs a cubic time complexity- and (nni) approximations wwhich can only estimate their computational error through potentially loose theoretical bounds which may [e.g. 21, 22, 44] or may not be computable in less than (n [4, 37]. At any point during a run of Algorithm 1, computing the predictive mean on no new data points has cost (n,,),, while the marginal predictive (o-)variance can be evaluated in ni) since C; is of rank i. Additionally, using Matheron's 's rule [45-47], sampling from the approximate posterior at no evaluation points also only requires (n, computation (assuming we can sample from the prior-see Section S3.3). The objects required to make predictions and draw samples are the vector v; and low rank matrix C₂ which both require O(n)) memory. Finally, the memory cost of Algorithm is only linear in n, since matrix multiplication v + Kr can be computed without explicitly forming K [48].
2.3 Related Work
eectoor vector multiplies, particularly CG [38], has become popular recently [5, 24-29, 39). Advances in specialized hardware has boosted their scalability without excessive memory footprint (27, 48]. Such iterative methods typically rely on preconditioning, which has been shown to significantly improve their performance [25, 26, 29]. Our method generalizes CG in this setting and thus retains the same benefits. At its core Algorithm 1 employs a (Bayesian) probabilistic numerical method [32-35], more specifically probabilistic linear solver (PLS) [49-54] applied to the linear system KU, = y. The fact that a PLS using CG actions can recover CG in its posterior mean was observed previously [49, 53). Here, extend this result to residual actions and preconditioning. Further, we also demonstrate the connection to the Cholesky and singular value decompositions. For randomized actions, the PLS as part of Algorithm also recovers the randomized Kaczmarz method in its posterior mean [55-58]. Employing PLS for GP approximation by updating beliefs over the kernel and precision matrix was suggested previously [53, 59]. Our work differs in that it updates a belief over the representer weights, as opposed to the kernel function matrix, considers more general projections than just conjugate residuals, and, most importantly, provides a theoretically motivated combined posterior which can be computed exactly.
GP inference based on matrix-
3 Theoretical Analysis
The main goals of our theoretical analysis will be to prove
(a) convergence of IterGP's posterior mean in norm (Theorem 1) and pointwise (Corollary 1) and to provide rigorous justification for the combined and computational uncertainty. Importantly, the combined uncertainty a right worst-case bound on the relative distance all potential latent functions consistent with our (computational) observations (Theorem 2).
We will demonstrate similar interpretation of the computational uncertainty as bound on the relative error to the mathematical posterior mean (see eqs. (14) and (16)).
3.1 Estimation of Representer Weights
At the heart of Algorithm 1 is probabilistic linear solver [49-51, 53] iteratively updating belief about the representer weights. It constructs an expanding subspace s 8;} = span(di..... d spanned by the actions in which the inverse K is perfectly identified. Each step d, expanding this explored subspace K Koorthogonal to the previous ones.
Proposition 1 (Conjugate Direction Method)
Let the actions of Algorithm I be linearly independent. Then Algorithum I is a conjugate direction method, i.e. holds that di kd, = 0 for all i $ j.
Proof. Without loss of generality assume i > j. Then the result follows directly from Lemma $1. Li
Geometrically, Algorithm 1 iteratively projects the representer weights onto the expanding subspace span( (S.} with respect to (.. We can use this intuition to understand the convergence of the representer weights estimate. The relative error p(i) at iteration i given by how small the "angle" between this subspace and the representer weights vector is.

---

Proposition 2 (Relative Error Bound for the Representer Weights) For any choice of actions relative error bound p(i), s.f. I|u. -ill - p()) v. ll is given by
( = (ج) (I (T C) u.)² Amx - 1
(9)
projection onto span(S,)
where 8. U./Vlu,11 k. If the actions {8: are linearly independent, then p(i) S 8n-i.
Proof. See Section S2.2.
Proposition 2 guarantees convergence in at most iterations, if the actions are chosen to be linearly independent, since C, K a K -orthogonna projection onto (pan(S.) (see Lemma S1). Therefore, if our finite computational budget is large enough, we eventually recover the mathematical posterior. This is reflected by the contraction of the posterior over the representer weights (see Proposition S4). The bound in Proposition 2 is tight without further assumptions an adverserial sequence actions such that the first (n - 1) are span J- Then the inverse is perfectly identified in that subspace, but tighter convergence bounds for specific sequences actions. For example, for randomized actions the bound depends on their distribution [56, 57]. If residuals T, chosen actions, we obtain
the actions, since there exists
C.y Ku 0. In practice, one can derive
p(i)
V Vx+1
p(i)
Ar An-
(10)
since then Algorithm I's estimate of the representer weights equals that of CG (Corollary S2). Here is the condition number and > the eigenvalues of either (i) the kernel matrix K if 8; rú (ii)
the preconditioned kernel matrix P KKP if P r
3.2 Convergence RKHS Norm of the Posterior Mean
Having established convergence of the representer weights estimate, we can use this result to prove convergence norm of IterGP's posterior mean to the mathematical posterior at the same rate.
Theorem (Convergence in RKHS Norm of the Posterior Mean Approximation)
Let Hu be the RKIS associated with kernel k(-:-), o² > 0 and let μaa E Hr be the unique solution to the regularized empirical risk minimization problem
arg fEH. 1 (f(x)) y₁ + +
(11)
which is equivalent to the mathematical posterior mean to shift by the prior μ le.g. 1, Sec. 6.2). Then for E (O....,n] n) the posterior mean pi(-) computed by Algorithm I satisfies
I|w. millm, S p(i)c(o?))lu4 poll
(12)
where No μ. the prior mean and the constant c(o2)= /1 1 1+
0²
a 0² O.
Xain
Proof. See Section S2.3.
Theorem 1 gives a bound on the RKHS -nom error between the posterior mean 14 of IterGP and the mathematical posterior mean M4.. If for the given prior kernel a bound on the RKHS- norm error || /.,1 Hy between the latent function and the mathematical posterior mean μ4, is known, Theorem can be directly used to bound the RKHS-norm enror between IterGP's posterior mean and
the latent function via the triangle inequality: || p4 pillz,
< |/ Mt 14l12,, + ||,.. M4 Mill ll
+0 as n-000
+0 as 1-30
3.3 Combined and Computational Uncertainty as Worst Case Errors
While Theorem 1 shows convergence in norm for IterGP's posterior mean, the convergence rate p(i) may contain expressions which cannot be evaluated runtime with the limited computation at our disposal. For example, for residual actions evaluating eq. (10) requires computation of the kernel matrix spectrum. However, the combined uncertainty of IterGP a tight bound on the pointwise relative error to all possible latent functions which would have resulted in the same computations.
7

---

Varismce
X
Latent Function
Mathematical Posterior Mean
Mathematical Uncertainty
Combined Uncertainty
Data
Approximate Posterior Meun
Computational Uncertainty
Figure 3: Computational and combined uncertainty of IterGP as worsi-case bounds."
Theorem (Combined and Computational Uncertainty as Worst Case Errors)
Let o² 2 O and let -)(-.-) )) k.(-.-) K kammp (..-) be the combined uncertainty computed by Algo- rithm 1. Then, for any æ E X (assuming & X if > 0) we have
of approximate posterior
sup gE :1gl Rg S1
g(x)
(x)
(x) - (x)
v (x, x:) o2
and (13)
of math, post. mean
computational
sup
[((x)) - (x)
k k:omp (xx, x)
(14)
9 Ha flll Hge
computational crror
where p-(()) kk., -(-
(X)) is the mathematical and p.) = k(., ( IteGP's
posterior mean for the latent function g Hyg. If o² then the above also holds for X.
Proof. See Section S2.4.
Theorem rigorously explains why the combined (mathematical computational) uncertainty k; the correct object characterizing our belief about the latent function h. given that we are in the limited compute regime. the same way that the mathematical uncertainty a tight bound on the distance to all functions g which could have produced the data (see [36, Prop. 3.8]), the combined uncertainty a tight bound on all functions g which would have produced the same computations.
3.4 Pointwise Convergence of the Posterior Mean
In particular, Corollary 1 shows and Figure illustrates, the computational uncertainty (0) is pointwise bound on the relative distance the mathematical posterior mean (16) and the combined uncertainty (O + )) is pointwise bound on the relative distance to the true latent function (15).
Corollary (Pointwise Convergence the Posterior Mean)
Assume the conditions of Theorem hold and assume the latent function h E Hye. Let p. be the corresponding mathematical posterior mean and μs the posterior mean computed by Algorithm 1. Then holds that
[h(x)) - -x(30)] S (25 x)) + 02,
and
(15)
|1 Haa
1p4. (xx) - Ms(x))]
comp
(x. x)).
(16)
||Al|
Proof. This follows immediately from Theorem 2 by recognizing that h/|lh has unit norm.
"The combined (co-)vand cco-)variance decomposes into mathematical and computational covariances, opposed to the combined standard deviation since - v + v The bottom panel thus illustrates the variance decomposition. However, to better illustrate Theorem 2, in the upper panel we plot the combined standard deviation / + and computational standard deviation / within it, in line with standard GP plotting practice.

---

It is worth noting that Theorem 2 and Corollary 1 generally do not hold for other GP approximations. They explicitly rely on K being the K Koorthogonal projection onto the space spanned by the actions (see Lemma S1). Since orthogonal projections are unique, if another GP approximation is such a projection and therefore satisfies Theorem 2, fact an instance of IterGP.
Experiments
To demonstrate the effects of quantifying computational uncertainty we perform GP regression on synthetic and benchmark datasets for the two most common GP approximations in the large- scale setting, SVGP [22] and CGGP [26], and their direct analogs from class of methods. An implementation of Algorithm 1, based on KeOps [48] and ProbNum [60], is available at: https: /gi thub. com/ JonathaniWegger/ttergg.
Experimental Setup We consider a synthetic dataset of iid uniformly sampled inputs 2j E [-1, 1j with y(x) /T: sin(TxTTT
+ E, where E ~ (O, d²), as well as range of UCI datasets [61] with training set sizes = 5, 287 to 57, 247, dimensions d = 9 to 26 and standardized features. All experiments were run on an NVIDIA GeForce RTX 2080 graphics card. We perform GP regression using a zero mean prior and a Matérn()) kernel (for other kernels see Section S4). All experiments were run times with randomly sampled training and test splits of 90/10 and we report average metrics with 95% confidence intervals.
IterGP reduces the necessary computations for CG-based GP inference. We compare IterGP to the CG-based GP inference used in the GPyTorch library [26]. For all datasets, we select hyperparameters using the training procedure of Wenger al. [29]. As we show in Theorem S5, the posterior mean of IterGP with (conjugate) residual actions is exactly equivalent to performing CG to compute the representer weights. Therefore, both methods produce the exact same posterior mean estimate and thus achieve the same RMSE as function of CG iterations (Figure 4, bottom). The primary difference between the two methods in the posterior variance. The combined variance estimate of IterGP is essentially "free" in the sense that it reuses terms from the posterior mean calculation. In contrast, computing the posterior variance with CG requires n. additional linear solves (K-ixo.,.... K K-.xon.)) GPyTorch relies on the Lanczos Variance Estimate technique [62] which essentially warm-starts each of these solves by reusing quantities from the linear solve K K-tk( k(x .F,1). While this approach produces reliable variance estimates that converge to the true posterior variance, it requires additional computation: at least one set of additional CG iterations to compute K kkX,, In Figure 4(a) (top). we see that IterGP and GPyTo'ch's 's CGGP achieve nearly identical NLL, suggesting that both methods produce variances that yield similar generalization.
Synthetic
Parkinson's
Bike Sharing
Protein
(11 1,024, di 5)
(n. 5,287, d' 21)
(a 15,641, = 16)
(n 41,157
Any Dataset
16-
1.200
400
20
TeeCIP-O
BerOPPO
1.175
10
woderr prodices
1.150
100
200
1.125
Mutris Muttii-aodorr
800
OGGP Tee:GP-C
RMSE
1.5
0.70
100
0.68
5-
1.0
0.68
10a
30
100
Iteration
Iteration
Iteration
(a) Generalization synthetic and UCT benchmark datasets.
(b) Comp. Cost
Figure 4: Generalization of CGGP and its closest IterGP analog. (a) GP regression using a Matérn kernel on UCI datasets. The plot shows the average generalization error in terms of NLL and RMSE for an increasing number of solver iterations. The posterior mean of IterGP-CG and CGGP is identical, which explains the identical RMSE. However, CGGP performs additional computation for the posterior covariance as (b) illustrates, which is not needed since IterGP-CG has identical NLL.

---

The key difference between the methods is that 1) unlike CGGP, IterGP's variances exactly capture both mathematical and computational uncertainty, and 2) IterGP's variances require no additional solves, resulting in half as much computation GPyTorch's CGGP implementation (see Figure 4(b)).
Quantifying computational uncertainty improves generalization of inducing point methods. To understand the benefits of quantifying computational uncertainty, we compare the linear-time SVGP method (which does not quantify computational uncertainty) with the closest (quadratic-time) inducing point analog from our proposed IterGP framework (see Section 2.1). While the IterGP method is inherently more expensive than SVGP, our goal is simply to demonstrate that inducing points can yield far more accuracy if one has the budget to account for computational uncertainty. To that end, we compare SVGP against IterGP using the same set of randomly -placed inducing points. We identify a set of kernel hyperparameters by optimizing the ELBO of SVGP on the training data, using these for both SVGP and IterGP. As Figure shows, we find that across all datasets that IterGP offers better RMSE and NLL than SVGP, despite the fact that the hyperparameters are chosen to favor SVGP. This suggests that the extra computation needed to quantify computational uncertainty can more "effectively" utilize set of inducing points for predictive models.
Synthetic
Parkinson's
Sharing
KEGGurdir
(1
1,024, d
(a
5,287, d 21)
(n
15,641,
16)
(R 41,157, d -9)
(n 57,247, 26)
3.5
1.25
0.00
SVGP IeerGP-P
1.20
3.0
1,00
26
1.15
0.70
-0.50
1.10
0.65
0.25
SVGP IterGP-PI
RMSE
0.75
0.60
0.70
0.15
2
0.55
102
10
10a
16a
10a
10a
10a
10a
102
10*
1nd. Points / Iseration
Ind. Points Iteration
Ind. Points Iteration
Ind. Points Iterationn
Ind. Points (teration
Figure 5: Generalization of SVGP and its closest IterGP analog. GP regression using Matérn ) kernel on UCI datasets. The plot shows the average generalization error in terms of NLL and RMSE for an increasing number of identical inducing points. After small number of inducing points relative to the size the training data, IterGP has significantly lower generalization error than SVGP.
5 Conclusion
Scalable GP approximations inevitably introduce error, leading to a worse model for the latent function in question. This work demonstrates that possible account for both uncertainty arising from limited data and uncertainty arising from limited computation exactly which as we show improves model performance. IterGP methods return this combined uncertainty which crucially represents a dataset-specifi pointwise worst-case bound on the error to the true latent function. At its core, IterGP performs repeated matrix-vector multiplication resulting in quadratic complexity. Since modern computing architectures (i.e. GPUs) have been specifically designed for this operation at scale, iterative approaches for GP approximation are becoming competitive with theoretically cheaper approximations, like inducing point methods [26. 27]. Finally, in addition to the general utility of IterGP, we expect this class methods to be particularly useful in applications where accurate uncertainty quantification is important or, due to its inherently online nature, where data is acquired sequentially such as in active learning and Bayesian optimization.
10

---

Acknowledgments and Disclosure of Funding
JW, MP and PH gratefully acknowledge financial support by the European Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: funds from the Ministry of Science, Research and Arts of the State of Baden-Württeen authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting JW and MP. JPC and GP are supported by the Simons Foundation, the McKnight Foundation, the Grossman Center, and the Gatsby Charitable Trust. The authors would like thank Hanna Dettki, Frank Schneider, Lukas Tatzel and all anonymous reviewers for helpful feedback on an earlier version of this manuscript.
and The
References
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processex for Machine Learning. The MIT Press, 2006.
[2] Huaiyu Zhu, Christopher KI Williams, Richard Rohwer, and Michal Morciniec. Gaussian regression and optimal finite dimensional linear models. In Neural Networks and Machine Learning. 1997.
[3] Giancarlo Ferrari Trecate, Christopher KI Williams, and Manfred Opper. Finite-limensional approximation of Gaussian processes. In Advances in Neural Information Processing Systems (NeurIPS), pages 218-224, 1999.
[4] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems (NeurIPS), 2007.
[5] Andrew Wilson and Hannes Nickisch. Keree interpolation for scalable structured Gaussian processes KP) In Proceedings of the 32nd International Conference on Machine Learning (ICML), 2015.
[6] Andrew Gordon Wilson, Christoph Dann, and Hannes Nickisch. Thoughts on massively scalable Gaussian processes. arXiv pre-print, 2015. URL http://tseeee //t::///rvvvv//// wtt//a/////aaeeeeeee
[7] Zichao Yang, Andrew Wilson, Alex Smola, and Le Song. carte-lcarning fast kernels. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), pages 1098-1106, 2015.
[8] Pavel Izmailov, Alexander Novikov, and Dmitry Kropotov. Scalable Gaussian processes with billions of inducing inputs via tensor train decomposition. In International Conference on Artificial Inteliigncce and Statistics (AISTATS), 2018
[9] Trefor Evans and Prasanth Nair. Scalable Gaussian processes with grid-structured eigenfunctions (GP- GRIEF). In International Conference Machine Learning (ICML), 2018
[10] Amir Zandich, Navid Nouri, Ameya Velingker, Michael Kapralov, and Ilya Razenshteyn. Scaling up kernel ridge regression via locality sensitive hashing. In International Conference on Artificial Intelligence and Statistics (AISTATS). 2020.
[11] Aldo V Vecchia. Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society: Series (Methodologicall, 50(2):297-312, 1988.
[12] Abhirup Datta, Sudipto Banerjee, Andrew O Finley, and Alan E Gelfand. Hierarchical nearest-neighhor Gaussian process models for large geostatistical datasets. Journal the American Statistical Association, 110(514):00-812,2 2016.
[13] Matthias Katzfuss and Joseph Guinness. General Framework for Vecchia Approximations of Gaussian Processes. Statistical Science, 36(1):124 141, 2021.
[14] Florian Schäfer, Mathias Katzfuss, and Houman Owhadi. Sparse Cholesky factorization by Kullback- Leibler minimization. SIAM Journal on Scientific Computing. 43(3),,2019-A2246, 2021.
[15] Alex J Smola and Peter Bartlett. Sparse greedy Gaussian process regression. In Advances in Neural Information Processing Systems (NeurIPS), pages 598-604, 2000.
[16] Christopher Williams and Matthias Seeger. Using the Nyström method speed up kernel machines. In Advances in Neural Information Processing Systems (NeurIPS), 2001.

---

[17] Petros Drineas and Michael W Mahoney. On the Nyström method for approximating a Gram matrix for improved kernel-based learning. Journal of Machine Learning Research, 6:2153-2175, 2005.
[18] Matthias W. Seeger, C. K. Williams, and N. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS). 2003.
[19] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. Advances in Neural Information Processing Systems (NeurlPS), 18:1257-1264, 2005.
[20] Joaquin Quinooero-Candee and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Research, 6:1939 11999-1959, 2005.
[21] Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), pages 567-574, PMLR, 2009.
[22] James Hensman, Nicolè Fusi, and Neil D Lawrence. Gaussian processes for big data. In Proceedings of the 29th Conference on Uncertainty in Artificial Intelligence (UAI), pages 282-290, 2013.
[23] Mark Gibbs. Bayesian Gaussian processes for classification and regression. University of Cambridge, Cambridge, 1997.
[24] lain Murray, Gaussian processes and fast matrix-vector multiplies. In Numerical Mathematics in Machine Learning Workshop (ICML), 2009.
[25] Kurt Cutajar, Michael Osbore,, John Cunningham, and Maurizio Filippone. Preconditioning kernel matrices. In International Conference on Machine Learning (ICML), 2016.
[26] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinherger, and Andrew Gordon Wilson, GPy- Torch: Blackbox -atrix Gaussian process inference with GPU acceleration. Advances in Neural Information Processing Systems (NeurIPS), 2018:7576-7586, 2018.
[27] Ke Alexander Wang, Geoff Pleiss, Jacob R Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gor- don Wilson, Exact Gaussian processes a million data points. Advances in Neural Information Processing Systems (NeurIPS), 32, 2019.
[28] Artem Artemev, David R Burt, and Mark van der Wilk. Tighter bounds on the log marginal likelihood of Gaussian process regression using conjugate gradients. In International Conference on Machine Learning (ICML), 2021.
[29] Jonathan Wenger, Geoff Pleiss, Philipp Hennig, John P. Cunningham, and Jacob R. Gardner. Precondition- ing for scalable Gaussian process hyperparameter optimization. In International Conference on Machine Learning (ICML), 2022.
[30] Matthias Bauer, Mark van der Wilk. and Carl Edward Rasmussen. Understanding probabilistic sparse Gaussian process approximations. Advances in Neural Information Processing Systems (NeurIPS), 29, 2016.
[31] Jonathan H. Huggins, Trevor Campbell, Mikolaj Kasprzak, and Tamara Broderick. Scalable Gaussian process inference with finite-data mean and variance guarantees. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2019.
[32] Philipp Hennig, Mike A. Osborne, and Mark Girolami. Prohabilistic numerics and uncertainty in computa- tions. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 421(2179), 2015.
[33] Jon Cockayne, Chris Oates, TJ Sullivan, and Mark Girolami. Bayesian probabilistic numerical methods. SIAM Review, 6164):756-789, 2019.
[34] Chris Oates and TJ Sullivan. A modern retrospective on probabilistic numerics. Statistics and Computing, 10 2019.
[35] Philipp Hennig, Michael A. Osborne, and Hans P. Kersting. Probabilistic Numerics: Computation as Machine Learning- Cambridge University Press, June 2022. ISBN 1881116661411 doi: 10.1017/ 9781316681411.
[36] Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaussian processes and kernel methods: A review on connections and equivalences. arXiv preprint, 2018. URL http: // aarxiv org/ a88// 12tp/2888...02.....
12

---

[37] Florian Schäfer, TJ Sullivan, and Houman Owhadi. Compression, inversion, and approximate PCA of dense kernel matrices at near-linear computational complexity. Multiscale Modeling and Simulation, 19 (2):88-770,, 2021.
[38] Magnus Rudolph Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49, 1952.
[39] John P Cunningham, Krishna V Shenoy, and Maneesh Sahani. Fast Gaussian process methods for point process intensity estimation. In International Conference on Machine Learning (ICML), 2008.
[40] Andres Potapczynski, Luhuan Wu, Dan Biderman, Geoff Pleiss, and John P Cunningham. Bias-free scalable Gaussian processes via randomized truncations. In International Conference on Muchine Learning (ICML), 2021.
[41] Bernhard W Silverman. Some aspects of the spline smoothing approach to non-parametrie regression curve fitting. Journal of the Royal Statistical Society: Series B (Methodological), 47(1):1-21, 1985.
[42] Lehel Csató and Manfred Opper. Sparse on-line Gaussian processes. Neural computation, 14(3):641-66) 2002.
[43] Veit Wild, Motonobu Kanagawa, and Dino Sejdinovic. Connections and equivalences between the Nyström method and sparse variational Gaussian processes. arXiv pre-print, 2021. URL http::/arxixi org/ abs/2106. 01121.
[44] David R. Burt, Carl Edward Rasmussen, and Mark Van Der Wilk. Rates of convergence for sparse variational Gaussian process regression. In International Conference on Machine Learning (ICML), 2019.
[45] Andre G Journel and Charles J Huijbregts. Mining geostatistics. Academic Press London, 1976.
[46] James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In International Conference on Machine Learning (ICML). pages 10292-10302. PMLR, 2020.
[47] James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of Gaussian processes. arXiv pre-print, 2020. URL http://arxiv: http://aaxvvv.rr/a/b/// /bs/// 04026.
[48] Benjamin Charlier. Jean Feydy, Joan Alexis Glaunès, François-David Collin. and Ghislain Durif. Kernel operations on the GPU, with autodiff, without memory overflows. Journal of Machine Learning Research (JMLR), 2- 2021. URL attp://jmlz: werg/papers/v22 wv2p////22///2////// 2200/200555 html.
[49] Philipp Hennig. Probabilistic interpretation of linear solvers. SIAM Journal on Optimization, 25(1): 264-260.22 2015.
[50] Simon Bartels, Jon Cockayne, llse CF Ipsen, and Philipp Hennig. Probabilistic linear solvers: A unifying view. Statistics and Computing. 29(6); 1249-1263, 2019.
[51] Jon Cockayne, Chris Oates, Ilse C. Ipsen, and Mark Girolami. A Bayesian conjugate gradient method. Bayesian Analysis, 14(3)937--012. 2019.
[52] Jon Cockayne, Ilse C. F. Ipsen, Chris J. Oates, and Tim W. Reid. Prohabilistic iterative methods for linear systems. arXiv pre-print, 2020. URL http:// wwtg////////eeooee//oo//eee 12615.
[53] Jonathan Wenger and Philipp Hennig. Probabilistic linear solvers for machine learning. In Advances in Neural Information Processing Systems (NeurIPS), 2020.
[54] Tim W. Reid, Isse C. F. Ipsen, Jon Cockayne, and Chris J. Oates. BayesCG as an uncertainty aware version of CG. arXiv pre-print, 2022. URL http://arxiv hwtp:/oog//aaa www////////////e/////oo 2008.. .03225.
[55] Stefan Kaczmarz. Angenäherte Auflösung von Systemen linearer Gleichungen. Bulletin International de I' l'Académie Polonaise des Sciences et des Lettres. Classe des Sciences Mathématiques et Naturelles. Série
A, Sciences Mathématiques, pages 355-357, 1937.
[56] Thomas Strohmer and Roman Vershynin. A randomized Kaczmarz algorithm with exponential convergence. Journal of Fourier Analysis and Applications, 15(2):262-278, 2009.
[57] Robert Gower and Peter Richtárik. Randomized iterative methods for linear systems. SIAM Journal on Matrix Analysis and Applications, 1664),6660-1600, 2015.
13

---

[58] Robert M Gower. Sketch and project: Randomized iterative methods for linear systems and inverting matrices. arXiv preprini, 2016. URL htttt:://aariiior/abbb/////s////
[59] Simon Bartels and Philipp Hennig, Conjugate gradients for kernel machines. Journal of Machine Learning Research, 21(55):1 2020.
[60] Jonathan Wenger, Nicholas Krämer, Nathanael Bosch, Marvin Pförtner, Jonathan Schmidt, Nina Ef- fenberger, Maren Mahsereci, Johannes Zenn, Thomas Gläßle, Alexandra Gessner, Toni Karvonen, Philipp Hennig, and ProbNum Contributors. ProbNum: Probabilistic Numerics in Python, 2021. URL https:/ / /tpp:://ittee thub.
[61] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL BtLLttt///rchhi// iis. uci. edu/ml.
[62] Geoff Pleiss, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson. Constant-ting predictive distributions Gaussian processes. Proceedings of the International Conference on Machine Learning (ICML), pages 4114-4125 4123, 2018.
[63] Yousef Saad, Manshung Yeung, Jocelyne Erhel, Frédéric Guyomare h. deflated version of the conjugate gradient algorithm. SIAM Journal on Scientific Computing, 21(5):1909-199) 1926, 2000.
[64] Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & Business Media, 2006.
[65] Jason Frank and Cornelis Vuik. On the construction of deflation-based preconditioners. SIAM Journal on Scientific Computing, 23(2):4424462, 2001.
[66] Jean-Paul Chiles and Pierre Delfiner. Geostatistics: modeling spatial uncertainty, volume 497. John Wiley & Sons, 2009.
[67] Arnaud Doucet. A note on efficient conditional simulation of Gaussian distributions. Technical report. Departments of Computer Science and Statistics, University of British Columbia, 2010.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] Proofs are included in the supplementary material.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as URL)? [Yes]
(b) Did you specify all the training details (e.g.. data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, or cloud provider)? [Yes]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
14

---

(d) Did you discuss whether and how consent was obtained from people whose data you' re using/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "rXpXrUxHKFm", "reviewer": "Reviewer_pTot", "summary": "The paper proposes an approximate inference technique, called IterGP, for Gaussian progress (GP) regression. Methodologically, the key algorithm (Algorithm 1) provably constructs increasingly more accurate approximations (of the GP posterior) as the runtime increases. Each approximation also comes with a notion of approximation error (in the sense of Theorem 2). Empirically, IterGP can be computationally cheaper, but just as accurate as existing methods (Figure 4), or more expensive but more accurate compared to existing methods (Figure 5).", "questions": "I would like the authors to address my questions regarding the “quality” subsection of the “Strengths and Weaknesses” section. My biggest concern is being confused about the relationship between combined uncertainty and the approximate posterior variance. The answer to this question is important, because if the combined uncertainty is the same as the approximate posterior variance, then the claim in the abstract “[…] due to limited computation, is entirely ignored when using the approximate posterior” is not accurate, since approximate GP automatically accounts for both computational and mathematical uncertainties, through the approximate posterior variance / combined uncertainty. More generally, I don't see the benefit of the so-called combined uncertainty, it it's just equal to the approximate posterior variance. My confusion also spawned my question about whether Theorem 2 applies to other approximate GP techniques. It also makes me question whether the claim in Figure 1 is valid: if it were true that Theorem 2 also applies to SVGP, then the message of modelling computational uncertainty improves GP approximation is incorrect. Perhaps the cause of IterGP-PI performing better than SVGP is due to the later using inducing points, which can be inaccurate.\n\nThe list of typos I’ve found are as follows\n- Lines 83-84: Equation 2 is true if the prior mean is zero.\n\n- Lines 96:97: the contraction only happens for a well chosen set of policies $s_i$, not in general\n\n- Lines 119-120: it’s not true that the algorithm is online. Line 9 of Algorithm 1 requires the whole kernel matrix $\\hat{K}$.\n", "limitations": "I do not see any negative societal implications to this line of work.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I think the paper’s key ideas are original and strong, but I have some concerns about the quality and clarity of the manuscript.\n\n# Originality\nSection 2 is a novel way of using probabilistic linear solvers (PLS) to perform GP regression. Section 2 treats the so called representer weights as unknown/random, whereas existing works treat either the kernel matrix or the inverse of the kernel matrix as random. The randomness here is due to limited computation. Since the posterior mean at a test point is a function of the representer weight (Equation 2), the randomness in the representer weight propagates into the posterior mean and is computable as in Equation 6. The update equations 4 and 5 give the general skeleton of how to reduce the randomness in the representer weights. \n\n# Quality \nIn terms of strengths, the visualizations and plots in the paper are very legible and convey the message well. Figure 1 panel b breaks down uncertainty into two constituent parts. Figure 4 demonstrate how IterGP is cheaper than conjugate gradient (CG) GP but just as accurate. Similarly, Figure 5 is very easy to read and understand the takeaway.\n\nIn terms of room for improvement, I think the paper is missing three kinds of discussions. \n- The first is the impact of the policy. Although Table 1 mentions a list of policies, the rest of the paper does not discuss the pros and cons of different policies, or give recommendations on which to use. I do understand that in a sense the actions s_i are hyper-parameters of Algorithm 1, and tuning hyper-parameters is potentially its own problem.  \n- The second is distinguishing the combined uncertainty (in Equation 6) from the regular approximate posterior variance at test points. What I mean by approximate posterior variance is that, while the true posterior variance at $x_*$ is \n$$k(x_*, x_*) – k(x_*, X) \\hat{K}^{-1} k(X, x_*),$$\nWhen we have an approximation of the inverse of the kernel matrix, say, $L \\approx \\hat{K}^{-1}$, then we announce the approximate posterior variance to be\n$$k(x_*, x_*) – k(x_*, X) L k(X, x_*),$$\nMy current impression is that the combined uncertainty is exactly the approximate posterior variance for the choice $L = C_i$, where $C_i$ is constructed using the IterGP iterative procedure. But more generally, all approximate GP approaches that have approximations to $\\hat{K}^{-1}$ will also have the combined uncertainty, and it’s equal to the approximate posterior variance that they report. \n- The final discussion is whether Theorem 2, in particular Equation 13, applies to other kinds of approximate GP inference techniques. Skimming through Theorem 2’s proof, there is no part of the proof that relies on the actual structure of iterGP, in the sense that the statement will hold if we replace $C_i$ with any other approximation of $\\hat{K}^{-1}$. In that case, since other approximate GP inference techniques have their own approximation of $\\hat{K}^{-1}$, their approximate posterior variance can also be used to give worst-case bounds in the sense of Theorem 2. \n\n# Clarity\nThe paper is well-written. There are some negligible typos that I will highlight in the “Questions” section.\n\n# Significance\nAdvances in scalable and accurate GP inference are important since GP regression is a staple of modern data analysis.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "K3tbPQFcg6", "reviewer": "Reviewer_Cb8D", "summary": "        This paper proposed a method to take into account extra uncertainty arising from using approximate GP solutions to scale GPs to large datasets. The proposed method is based on making inference about some parts of the computation of the posterior predictive distribution of a full GP. The authors indicate that the proposed approach is a generalization of many GP approximate methods. Theoretical results associated to the proposed method and derived and some experiments are carried out on synthetic and real-world data.\n", "questions": "It is not clear that the marginalization of v* below Eq. (2) gives the GP prior. The mean would be zero mean after the marginalization, using Eq. (2). The resulting variance will have to add the variance of the mean. Can the authors clarify this? It seems Eq. (2) is missing the GP mean.\n\nAbove Eq. (3) should not r_i bet r_i-1?\n", "limitations": "The authors have indicated that they have described the limitations of their method. However, they do not indicate where in the manuscript.\n", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Weaknesses:\n\n        - The explanation of the proposed method is poor. The authors have to make a better job at explaining their method. Eq. (3) and the update on p(v*) has to be better explained.\n\n        - Section 2.1 is poorly explained. It is not clear the connections with other methods.\n\n        - The experimental section is a bit weak.\n\nStrengths:\n\n        - Nice theoretical results are derived for the proposed method.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "ElqZ8lzBl8W", "reviewer": "Reviewer_5Mot", "summary": "The authors present a novel method, named IterGP, for doing approximate inference in Gaussian processes. Taking a probabilistic numerics approach, they treat the representer weights, $\\mathbf v_* = \\hat{\\mathbf K}^{-1} (\\mathbf y - \\boldsymbol \\mu)$, as a quantity to iteratively compute. By building on recent advances in probabilistic linear solvers, they derive an algorithm for iteratively updating a distribution over the representer weights, which captures the (computational) uncertainty in them. \nBy reparameterising the GP posterior in terms of the representer weights and then marginalising them out, the authors obtain a new GP posterior expression, which implicitly accounts for the uncertainty in approximation made by the iterative representer weights.\n\nWhile explicitly computing the computational uncertainty is expensive, computing the combined (computational and mathematical) uncertainty has a quadratic cost, meaning that IterGP computationally sits in-between exact GP inference (cubic) and linear-time approximation methods, which, however, do not provide calibrated uncertainty estimates. At the same time, the memory cost of IterGP is only linear.\n\nTo make the iterative updates, IterGP requires a policy for selecting and weighing data. The authors discuss different strategies and show how they correspond to standard approximation methods, such as conjugate gradients and inducing point methods, and extend these to properly account for computational uncertainty.\n\nThe paper concludes with a theoretical analysis, proving general convergence of the posterior mean as well as proving that the combined uncertainty is a worst-case bound on the error, and an empirical evaluation of IterGP, showing that it provides better posterior estimates than commonly used methods.\n", "questions": "NB I have quite limited experience with probabilistic numerics, so the following questions should be read in this context! :)\n\n1. I don't fully understand the definition of the residual in line 90. Specifically, why does it contain $\\hat{\\mathbf K}$ and not just $\\mathbf K$? I would have expected the latter if it's the residual between the observations and the GP posterior mean.\n\n2. I found the discussion of the overconfidence in SVGP extremely interesting, but I don't quite understand the point made in lines 142-146. Since the representer weights of the inducing points are typically calculated exactly, where does the error in the uncertainty estimate (i.e., the overconfidence) come from? (To be clear, this is absolutely not a crucial question; I'm just trying to understand SVGP :) ).\n\n3. I suppose the choice of a stopping criterion will be problem-specific, but do you have some suggestions for this? While the computational uncertainty is too expensive to explicitly compute in general, is it possible to somehow estimate it and stop when we have reached a pre-defined tolerance?\n\n4. Figure S5 (b) in the supplementary hints are instabilities in certain settings. Do you know which parts of the algorithm lead to this, and can you say something more about when this could happen?\n", "limitations": "Potential negative societal impacts have not been discussed, but it is not necessary given the paper's theoretical nature.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\nThe paper overall is of very high quality. The writing is clear, and the figures are fantastic. I especially like the pedagogical use of colours for the different kinds of uncertainties, which are shared between figures, equations, and the text.\n\nFrom a technical perspective, the paper is very thorough, and the authors carefully explain the intuition behind the method where possible. They also provide a good discussion of the different policies for choosing actions and how these link back to commonly used approximation methods, showing a kind of unification (and extension) of these methods.\n\nThe theoretical analysis seems solid, but I did not check it in detail.\n\nWhile the field of probabilistic numerics is not my area of expertise, the work seems original and highly novel to me. Overconfident posteriors are a well-known issue in the approximate GP world, so the contributions in this paper make it a very exciting piece of work. It is definitely a significant contribution to the approximate GP literature, which will be of major interest to the NeurIPS community.\n\n**Weaknesses**\n\nWhile the writing is overall of very high quality, there are a lot of new concepts to digest for someone not in the field of probabilistic numerics (but who does have experience with GPs). A short background on probabilistic numerics (perhaps probabilistic linear solvers in particular) would have been helpful for me, but I understand that I might not be the target audience.\n\nOne thing that the authors did not discuss is how to choose the stopping criterion in algorithm 1. I suppose there are some trivial choices (e.g., if one has a fixed computational budget), but perhaps one can be smarter here. A brief discussion of this would have been nice.\n\nFor the experiments, it would have been interesting to see SVGP and IterGP compare in terms of wall-clock time for different numbers of inducing points/iterations. Since SVGP is faster to compute, I suppose it would achieve better results for very few inducing points. It would be interesting to see when IterGP begins to perform better for a fixed computational budget.\n\nAlso, it would be interesting to see the RMSE and NLL obtained by an exact GP indicated in both figure 4 and 5.\n\n**Various minor things**\n* In equation (2), I think a $+ \\mu$ is missing in the mean.\n* The next-to-last sentence in the caption of figure 2 appears incomplete.\n* A couple of lines in algorithm 1 are coloured grey, but there's no discussion of what this means, as far as I can see.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "5aflQp4nTU", "reviewer": "Reviewer_r5qg", "summary": "Gaussian process regression model is impractical for big data and various approximations have been proposed to alleviate such difficulty. However, people have overlooked the uncertainty due to the numerical approximations for saving computational resources. This paper presented a novel low-complexity algorithm that respects uncertainties arising from the finite number of data observed and the finite amount of computation expended. Besides, theoretical analyses are given to support the effectiveness of the proposed algorithm. Experimental evaluations were also conducted to verify the proposed algorithm with some widely used approximation methods, such as sparse GP. ", "questions": "1. Would the proposed algorithm work for distributed GP as well? If yes, what is the action $s$_i?\n2. The theoretical analysis assumes that the actions $s$_i are linearly independent. Is this always true for the considered approximations?\n3. Have you verified that the \"computational uncertainty\" term will go to zero when the computational resource is infinite?\n4. Do you have more experimental results obtained using other kernels, such as the spectral mixture kernel and some hybrid kernels?\n5. In the conclusion, you mentioned that the proposed algorithm is particularly useful for online data processing, which is not clear to me? \n", "limitations": "Maybe more experiments with different kernel functions and different data lengths would be helpful. ", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. A novel GP approximation that accounts for computational uncertainty.\n2. Well-crafted performance analyses (Th.1 and Th.2)\n3. Experimental evaluations in terms of a few widely used GP approximations. \n\nWeaknesses:\n1. In Eq.(6), I understand that the \"computational uncertainty\" term results in a combined uncertainty that is cheaper to evaluate; however, I doubt this term is solely due to the numerical approximation made to the standard GP. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}]

Author rebuttal:

Decision: Accept