Pointwise uncertainty quantification for sparse variational Gaussian process regression with a Brownian motion prior
Authors: Luke Travis
Kolyan Ray
Venue: NeurIPS 2023 poster (2023)
Keywords: Gaussian process
sparse variational Bayes
uncertainty quantification
theoretical guarantees
TL;DR: We study theoretical guarantees for pointwise uncertainty quantification using a sparse variational Gaussian process method based on a Brownian motion prior.

Abstract: We study pointwise estimation and uncertainty quantification for a sparse variational Gaussian process method with eigenvector inducing variables. For a rescaled Brownian motion prior, we derive theoretical guarantees and limitations for the frequentist size and coverage of pointwise credible sets. For sufficiently many inducing variables, we precisely characterize the asymptotic frequentist coverage, deducing when credible sets from this variational method are conservative and when overconfident/misleading. We numerically illustrate the applicability of our results and discuss connections with other common Gaussian process priors.

Pointwise uncertainty quantification for sparse variational Gaussian process regression with a Brownian motion prior
Luke Travis
Kolyan Ray Department Mathematics Imperial College London rolyyan.ry@iiperr... ac.
Department of Mathematics Imperial College London
luke. trwvisi1eimeerrillll...... acc
Abstract
We study pointwise estimation and uncertainty quantification for a sparse varia- tional Gaussian process method with eigenvector inducing variables. For a rescaled Brownian motion prior, we derive theoretical guarantees and limitations for the frequentist size and coverage of pointwise credible sets. For sufficiently many inducing variables, we precisely characterize the asymptotic frequentist coverage, deducing when credible sets from this variational method are conservative and when overconfident/misleadtie
We numerically illustrate the applicability of our results and discuss connections with other common Gaussian process priors.
Introduction
Consider the standard nonparametric regression model, where observe training samples D, ((11,3),..... (In, (n) arising from the model
f(x;) + +
3d N(O,00),
where o² > 0 and the design points x; E X' Rd are either fixed considered i.i.d. random variables. Our goal is to predict outputs y" based on new input features x", while accounting for the statistical uncertainty arising from the training data. A widely-used Bayesian approach is to endow f with Gaussian process (GP) prior 225, which is especially popular due to its ability to provide uncertainty quantification via posterior credible sets.
While explicit expressions for the posterior distribution are available, a well-known drawback is that these require ( time and Onn²) memory complexity, making computation infeasible for large data sizes n. To avoid this, there has been extensive research on low-rank GP approximations, where chooses < inducing variables to summarize the posterior, thereby reducing computation to Onn²²) time and (mm) memory complexity, see the recent review (00).
We consider here the sparse Gaussian process regression (SGPR) approach introduced by Titsias [1I, which is widely used in practice (see
9 for implementations) and has been studied in the number inducing variables m relatively wedl--nderstood, fewer theoretical results are available understanding how affects the quality of statistical inference. In particular, it is crucial to understand how large m needs to be to achieve good statistical uncertainty quantification, ideally close that of the true computationally expensive posterior. The present work contributes this research direction by establishing precise theoretical guarantees for pointwise inference using a natural choice of SGPR with eigenvector inducing variables.
I
many recent works [13 21 5 38
28
32
22
23 While the computational dependence
main appeal of Bayesian methods their ability to perform uncertainty quantification (UQ) via credible sets, often at the 95% level, see Figure T for visualization involving full GP and two
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Comparison of the Posteriors
0.4
0.2
0.0
02
X
Posterior
SGPR
SGPR (m =006)
Figure 1: Plot of the (sparse) posterior means with the ribbons representing 95% posterior pointwise credible intervals for the i) SGPR with m = 5 inducing variables (red); ii) SGPR with m = 106 inducing variables (green); iii) full GP (blue) and iv) the function (black). Here, n 500 and the covariance kernel is the rescaled Brownian motion with = 0.5.
SGPR approximations. We focus here on the canonical problem of pointwise estimation and UQ for f(x) at single point x E X, which corresponds to the practically relevant problem of quantifying the uncertainty prediction (() based on a new input æ. We study the frequentist behaviour of the resulting SGPR credible sets, with our main conclusion being that for well-calibrated priors and sufficiently many inducing points m, SGPRs can indeed provide reliable, though conservative, uncertainty quantification.
Concretely, consider a rescaled Brownian motion prior that models Hölder smooth functions 33 36 29
and whose pointwise UQ for the corresponding full posterior was studied in [29]. We extend this analysis the SGPR setting using the eigendecomposition of the sample covariance matrix as inducing variables, i.e. the optimal rank-m posterior approximation. We provide rigorous theoretical justification for pointwise estimation and UQ using this SGPR, provided m large enough. Our main theoretical contributions are summarized that in dimension d = 1, for an a-smooth truth inducing variables:
fo. 7-smooth prior, ? (0, 1]. and least >> n (1
1. The SGPR recovers the true value fu(a)) point with convergence rate a 7, this gives the minimax optimal rate n I
2. For truths that are smoother than the prior (a > y). SGPR credible sets are conservative with their frequentist coverage converging and one, which we characterize.
value strictly between their nominal level
3. For truths that are rougher than the prior (a S 7). coverage can be arbitrarily bad.
4. It is not necessary for the KI. divergence between the SGPR and true posterior to vanish as n + oo to get good (though conservative) coverage.
5. We provide a more general small-bias condition under which the coverage SGPR credible sets is at least that of credible sets from the true posterior.
Our proofs exploit specific properties of Brownian motion to obtain explicit expressions, which allow us to obtain the above precise results. Such expressions seem difficult to obtain for other commons GPs (e.g. Matérn, squared exponential). meaning that existing theoretical guarantees for the Matérn SGPR apply enlarged credible sets [32 23. Since our results indicate that pointwise SGPR credible sets are already conservative, such enlargements make the SGPR UQ procedure very

---

conservative, see Section 4 We provide complimentary perspective to these works, studying me credible sets from a related GP, rather than modified credible sets from the Matérn.
The Brownian motion prior we consider here shares many similarities with the Matérn process, meaning our main qualitative conclusions also apply there, as we demonstrate in simulations. We thus verify the applicability of our limit theorems in numerical simulations, confirming they provide good guide to the behaviour of SGPRs for finite sample sizes, including those based on Matérn priors. To the best of our knowledge, our work is the first establishing exact theoretical coverage guarantees for SGPR pointwise credible sets, and can be viewed as guide of what to expect for the other Gaussian processes, see Section
Related work. Rigorous theoretical guarantees for Bayesian UQ using the full posterior have received much recent attention. In parametric models, the classic Bernstein-von Mises theorem [35 ensures that regular credible sets are asymptotic confidence sets, justifying Bayesian UQ from a frequentist perspective. However, such a result fails to hold in high- and infinite-dimensional settings [10] 16] such as GPs, where the performance is sensitive to the choice of prior and model. Theoretical guarantees for UQ using full GPs have been established mainly involving L2-type credible sets (e.g. which, while mathematically more tractable, are less reflective of actual practice than pointwise credible sets. Results for pointwise Bayesian UQ are more difficult to obtain and are correspondingly much rarer [29 40 39], see Section 3 From a technical perspective, we perform
17 7 30)
careful and novel analysis of the gap between the pointwise behaviour of the true posterior and its SGPR approximation, after which we can invoke ideas from the proofs for the full posterior 29
While the frequentist asymptotic properties of VB have increasingly been investigated, this is usually in the context of estimation, convergence rates approximation quality (e.g. (24 2
26]
including for SGPRs [5] 6 (22) rather than UQ. Available results on mean-field VB often show that this factoriazable approach provides overconfident UQ by underestimating the posterior variance 3,4. though there exist more positive results in certain low-dimensional models [37]. Regarding UQ for SGPRs, Nieman al. 2 investigate the frequentist coverage of global Leacredible sets from certain SGPRs, showing these can have coverage tending to I, possibly after enlarging them. Vakili al. [32) similarly show that enlarged pointwise or uniform credible sets can have good coverage in a different bandit-like setting.
Organization. In Section 2 detail the problem setup, including notation and an overview of (sparse variational) Gaussian process regression. Main results on pointwise convergence and UQ for SGPRs based on Brownian motion prior found in Section discussion the connections with other GPs in Section 4 provide the proofs of our results (Section A) and additional simulations (Section B
simulations in Section 5 and discussion in Section In the supplement, we
2 Problem setup and sparse variational Gaussian process regression
Recall that we observe n training examples D,, {I11 B1)...,, }n)} from model II and write y = (y1,-... (U,)"-....)) and f = .,, We denote by Pj the probability distri- bution of from model (1) and by F, the corresponding expectation. For u, v E R we write {u, v) >n u,v; for the inner-product on IR" and ||u|| (u, 2/2 for the usual Euclidean norm. Write C(X) for the space of continuous real-valued functions on X. For a E (0. 1]. we further define the class of Hölder smooth functions on [0, 1] by C C [0, 1] = {f [0, 1] R
supzzy If(x)-f(y) |
cc). Throughout the paper we make the following frequentist assumption:
Assumption 2.1. There is a true fo E (x) generating the data y ~ Pfo according to ①..
We differentiate f coming from the Bayesian model with the 'true' generative fo-
2.1 Sparse variational Gaussian process regression
In model T]. we assign to ʃ GP k) a Gaussian process (GP) prior with mean function Up : X + Th. and covariance kernel k X' X + TR. We henceforth take the prior mean Mo for simplicity since this will not alter the spirit of our results, see Remark 3.3 below. The resulting
3

---

posterior distribution is again GP by conjugacy [25] with mean and covariance functions equal to
(n) = k, (x)² (Knn + I ' y.
(2)
(2. 2') = k( (n)x)
(K))"((nn (Knn + Iz, k,, (z').
where k,, = (k(2, (., x(,, and [mlli = k(r;; Zz). In particular, the posterior variance at a point x E X' is o (2)) km (2t, 2). In Section 3.2 below, we will specifically consider the choice X [0, 1] and rescaled Brownian motion prior kernel k(x, 2') (n+1/2) 1T min(x, x') for y > 0.
We consider here the sparse variational Gaussian process (SGPR) approximation using inducing variables proposed by Titsias [3I]. The idea is to summarize the posterior via relatively small number of inducing variables
{ui..... U. m < n, which linear functionals of the prior process f. By assigning u multivariate normal distribution, one obtains a low-rank variational family. Titsias [3T explicitly computed the minimizer in Kullback Leibler sense between this family of GPs and the full posterior, which yields the SGPR Q" = GP kim) with mean and covariance
Um (x) (om
KmnKnm) ""
km (m, x²) k(2, x') K KmmKmin Kmx + Kem (Kmm
Km Kn)) - K mar's
where Kum = (cov(f(x), u.).....ov((((?? um)) = KT M [Kmmli = cov(u, Uy) and we write [Kamli) cf,(f(x;)) u;). Recalling that cov (f))))), f(x)) k(x, under the prior, for u, of the form új - -im (2xa) (see 4 below) we obtain cov(f(z), my) = Eÿ= (2), The SGPR variance at x E X is o (z) = km (z, z). The computational complexity of obtaining Vm and km is O(nm²), which is much smaller than the O(n²) needed to compute the full posterior for &
2.2 Eigenvector inducing variables
The approximation quality of the SGPR depends on the inducing variables being well chosen. We consider here arguably the conceptually simplest and most canonical sparse approximation based on the cigendecomposition of K.n. Let 2 μ2 2 ... 2 2 0 be the eigenvalues of the positive semi-definite matrix Kan with corresponding orthonormal eigenvectors Vi...., Vn. Writing
v, (u)....,,,) " E R", we take as inducing variables
" ===
Σ sf.)) j j=1,...,m, m,
so that ūj is linear combination of the inducing points placed each data point with weights proportional to the entries of the jar-estt eigenvector of Knn- Setting V.. [Vi... Vm] Raxm, this choice of inducing points implies
Kym (x)) V.n
(
K diag(( Mm), Kmm Vm.diag(m)........ kmm) K mns
see Section C.1 of (5) Writing Knn. = Ei v.v.v? and substituting these expressions into 3 we obtain the SGPR with mean and covariance functions
(x) = (x)
Mr VVVAV
y
(5)
km (x,x²) = k((r,x') kn (x))
nc V k,(,'),
(6)
where 77 J24p We study inferential properties of the SGPR with mean 5 and covariance (6)
Comparing these expressions with those for the full posterior ②p. we see that this SGPR the optimal rank-m approximation in the sense that we have replaced (K + - n/..v?'' in both the mean and covariance (2) with D V corresponding to the largest eigenvalues of K Knn- These inducing variables are an example of interdomain inducing features [19], which can yield
sparse representations in the spectral domain and computational benefits [1]. Computing this SGPR

---

We further define the involves numerically computing the first m eigenvalues and corresponding eigenvectors of K which can be done in O(n" (n) time using for instance Lanezos iteration [8]]
This choice of SGPR using eigenvector inducing variables has been shown to minimize certain upper bounds for the Kullback-Leible (KI.) divergence between order SGPRs and the true posterior [S]. reflecting that it can provide good posterior approximation. Indeed, one can show that for m m. growing fast enough (though sublinearly in n), this SGPR converges the posterior in the KL sense as n + oc [5] Under weaker growth conditions on m, this SGPR still converges to the true generative fo at the minimax optimal rate in the frequentist model [22]. even when the SGPR
does not converge to the true posterior as n + 00. Our results include both scenarios and will show
that one can still obtain good frequentist UQ performance even for m small enough that the SGPR
diverges from the full posterior (see Remark 3.7 unlike in [5].
3 Main results
3.1 General results on pointwise inference
We now present our main results concerning estimation and uncertainty quantification at a point a E X' using the SGPR Q" with eigenvector inducing variables. Since both the posterior and SGPR are GPs, their marginal distributions (2)((z), Tm (x)))) respectively. It thus suffices to study these quantities under the frequentist assumption that there is a 'true' fo generating the data, i.e. Assumption 2.1
satisfy f (N),(z)), (x) and f(x)
frequentist bias and variance of the (sparse) posterior mean:
(x) E - fo(x)),
(x) = Ejo (u (2) fo(z)),
e. (1) = var (un (x));
2 (1) = var (Mm
where the expectation and variance are taken over Pf.. We have the following useful relationships between these quantities for the variational posterior and full posterior, which assumes fixed design, so that all statements are taken conditional on the design points ø .I....,,p take the noise variance to be 0² = 1 in (I) for our theoretical results.
In E X. For simplicity, we
Lemma 3.1. For fo E (2), x E X' and F (x)) = Σ R
Ej=m+i 1 VV,,
kn (x), we have:
bm (x) = (x) - (fm(z), fo)
(7)
e (x) = (x) [x)(())] - (z)
(8)
o (x) = (x) + (x), (x))
(9)
The quantity (x) measures the rank gap due to ignoring the smallest n - mm eigenvalues of K The SGPR approximation increases the posterior variance o (a:) > of (2) while reducing the frequentist variance of the posterior mean (x) S (x).. Thus if one can control the biases in (7). credible sets from the SGPR will have coverage least as large as the full posterior. One can therefore ensure good coverage by taking m = m. sufficiently small, but this can still lead to poor UQ by making the resulting credible intervals extremely wide, and therefore uninformative. We thus first study the effect of m on the quality of estimation, via the posterior convergence rate.
Proposition 3.2 (Pointwise contraction). For fo E C(X), x E X' and m = mn - oC,
F ( - fo(x)] > Mnem - 0,
as n oC, where M, oc is any (arbitrarily slowly growing) sequence and
2 = (x) + (x) + oz + (2)))))))))))))))))()) fo + Kr (2), k,), (x)))
If (Im(()) fo) (lb (r) and (I), (T), (x)) = o(o as n + oc, then the rate matches that of the true posterior b bz (z) (T) of
(z) (z)
The last result says that the SGPR contracts about the true fo at rate Em, L.e. it puts all but vanishingly small amount of probability on functions for which [f(x)= fu(z)] < Ma where fo is the true function generating the data. Such results not only quantify the typical distance between
point estimator fzz (e.g. SGPR mean/median) and the truth ([I)) Theorem 8.7), but also the typical
5

---

spread of Q' about the truth. Ideally, most of the Q" Q"-probability should be concentrated on such sets with radius e proportional to the optimal (minimax) rate, whereupon they have smallest possible size from an nnnrmation-theoreed perspective.
Contraction rates in global I.2-losses have been considered in the variational Bayes literature 24 21 27 26 22].. Such results crucially rely on the true posterior concentrating exponentially fast about the truth (see e.g. Theorem 5 in (26)), ) which happens in losses which geometrize the statistical model [5]. This is known not to be true for pointwise or uniform loss [5], and hence these previous proof approaches cannot be applied to the present pointwise setting. We must thus use an alternative approach, exploiting the explicit Gaussian structure of the SGPR.
The SGPR approach in UQ is to consider smallest pointwise credible set of probability 1 - 8:
Chh = C (D,) [1)(2)) - 21-80m (x)). Mm (2) + 21-50mm (x)];
(10)
with (( 1)J S 21 6) = - 6. We next establish guarantees on the frequentist coverage of Cs. Ó. ic. the probability 1 1-6 8' with n.F Pso (fo(z)) E C. 2 1- 8' over a target function class F. Remark 3.3 (Zero prior mean). Since Jn a-priori unknown, such coverage guarantees must hold uniformly over function class F to be meaningful. For symmetric function classes, where f E F implies- -f E F, taking a good prior mean vo for f will make estimating -f correspondingly more difficult. Thus without loss of generality we may take zero mean = O.
Proposition 3.4 (Pointwise UQ). Suppose that (am( (a)), fo) = o(bn (x)) as n - oo. Then for mn S and any å E (0, 1), the frequentist coverage of the credible sets satisfies
lim inf Pj. (f)()) E C.,) 2 lim Pj, (fo(z) E C.).
and (rm( (x)) olo2(x)), then If in addition (x)] = ol()(=))
lim (fo((( E C lim P%(fo(z) e C.)).
The first conclusion says that if the rank gap smaller than the posterior bias (a 'small-bia condition), SGPR credible sets will have asymptotic coverage at least at the level of the original full posterior. though perhaps larger, including possibly tending to one. The second conclusion gives further conditions on the rank gap under which the asymptotic coverages will be the same. Together with Proposition 3.22 the goal is to obtain a credible set Co of smallest possible diameter subject to having sufficient coverage.
3.2 Fixed design with rescaled Brownian motion prior
We now apply our general results specific GP and find conditions on the number of inducing variables m needed to get good UQ Consider the domain X' = [o, 1] with regularly spaced design points I; = *+172 for i = 1,..., n. For B a standard Brownian motion, consider as prior f =
Vc,, B, where Cn (n + 1/2) F for 7 O. Thus f is mean-zero GP with covariance kernel 6(x,x') = (x A z²), where x A x' := min(x: xx, x²). The scaling factor e, controls the smoothness of the sample paths of the GP and plays the same role as the lengthscale parameter for stationary GPs. The present rescaled Brownian motion is a suitable prior to model Hölder functions of smoothness Y E (0. (e.g.
33
29
). In particular, for a true fu C a E (0. 1], one obtains (full)
BAY posterior contraction rate n sy+1 in the global L-ssoss 333. We sometimes write Q" = Q for the corresponding SGPR to make explicit that the underlying prior is y-smooth.
Theorem 3.5. Let fo (O 1J. a E (0. and E (0, 1). Consider the SGPR Q7, with rescaled Brownian motion prior of regularity (0, and m Mm + do inducing variables. Then
1/3-x + n 1+27 m 140 + n T427 m 3/2)) + 0
170 (((x))](())]]))))) > M,,r 1+27
-
as n - 0o, where Mn. - do is any (arbitrarily slowly growing) sequence. If mn >> n 1+27 (2t) then the contraction rate of the SGPR matches that of the full posterior 1427.
Theorem 3.5 shows that for m >> n 1+ inducing points, the SGPR attains the same pointwise contraction rate as the full posterior. If ạ = a, we then recover the optimal (minimax) rate of

---

convergence for C smoothh function. For instance if a = y= 1, our result indicates that one can do minimax optimal pointwise estimation based on mn > vn inducing variables, a substantial reduction over the full n observations. Note this convergence guarantee can still hold in cases when the SGPR diverges from the posterior as n - oc, see Remark 3.7
The restriction a E (0, 1] comes from rescaling the underlying Brownian motion 33 One can extend this a > 1 by considering smoother baseline process [33, such as integrated Brownian motion, but the more complex form of the resulting eigenvectors and eigenvalues of Knn make our explicit proof approach difficult. Note also that the prior fixes the value f(0) 0. One can avoid this by adding an independent normal random variable to f. but since we consider pointwise inference at a point æ > 0, this will not affect our results and we thus keep f() = 0 for simplicity. Nonetheless, the message here is clear: for sufficiently many inducing variables (but polynomially less than n), an SGPR based on 7-smooth prior can estimate an a-smooth truth at rate 1127, We next turn UQ. Theorem 3.6. Let fo E (a (0, 1]. a E (0, and x E (), Consider the SGPR Q", with rescaled Brownian motion prior of regularity y (0, and m mn >> n n142n (77%) inducing variables. For 5 E (0, 1), let o := Pjo (Jo E C (I)) denote the frequentist coverage of the 1 - credible set C (x) given by 10 Then as nt 0o:
(i) (Undersmoothing case) If a > % then En (O 1/2)| S 21--) ps > 8 for all f E CO 1I. where P([N (0, 1)1 S 21-6) = 1 and (0, 1/2) is the normal distribution with mean zero and variance 1/2.
(ii) (Correct smoothing case) If = 7, then for each p (0,ps), there exists f SECO Co 1] such that en + p.
(iii) (Oversmoothing case) I < % there exists f E C [0, 1] such that gm + 0.
Theorem 3.6 provides exact expressions for the asymptotic coverage when the credible interval Co has diameter O(n -+ン) In the undersmoothing case a > y, the 1 8 SGPR credible sets are conservative from frequentist perspective, i.e. their coverage converges to a value strictly between 1-8 8 and 1. For instance 8 = 0.05 (0.1), ). the 95% (90%) credible set will have asymptotic coverage 99.4% (98.0%). indicating one does not want to enlarge such credible sets if using enough inducing variables. A desired asymptotic coverage can also be achieved by targetting the credibility according to the formula in Theorem 3.6(i). Reducing Y towards zero will generally ensure coverage for Hölder functions, but this will also increase the size of the credible intervals size O(n 1125 ), making them less informative. Note that the convergence in (i) is uniform over -blls of fixed radius.
In the oversmoothing case < 7, the posterior variance underestimates the actual uncertainty and so the credible interval is too narrow, giving overconfident (bad) UQ for many Hölder functions. In the correct smoothing case a = 7, where we obtain the minimax optimal contraction rate n 142a, coverage falls between these regimes- it does not fully tend to zero, but one can find function whose coverage is arbitrarily bad, i.e. close to zero.
The best scenario thus occurs when the prior slightly undersmooths the truth, in which case the SGPR credible interval will have slightly conservative coverage but its width Onn -42x) is not too much larger than the minimax optimal size O(n rizs Our results match the corresponding ones for the full computationally expensive posterior and the main messages are the same: undersmoothing leads
conservative coverage while oversmoothing leads poor coverage 22
Remark 3.7 (Kullback-Leibber)), the number of inducing variables grows like n 1127
< m E
1+2T, then the conditions of Theorems 3.5 and 3.6 are met, but the KL-divergence between the SGPR and the true posterior tends to infinity as n + oo (Lemmal AA66 in the supplement). Thus one does not need the SGPR to be asymptotically exact approximation of the posterior get similarly good frequentist pointwise estimation and UQ. In particular, one can take m polynomially smaller in than is required for the KL-divergence vanish, see [5] for related bounds and discussion.
While (5) consider approximation quality as measured by the KL-divergence between the SGPR and true posterior, we consider here the different question of whether the SGPR behaves well for pointwise inference, even when not a close apppooiooinn The data-generating setup [5) is also not directly comparable to ours. They take expectations over both the data (x, y) and prior
7

---

on f. so their results are 'averaged" over the prior, which downweights the effect of sets of small prior probability. In contrast, the frequentist guarantees provided here hold assuming there is a true function fo which generates the data according to (T) (Assumption] 2.1). meaning that one seeks to understand the performance of the method for any fixed fo. including worst-case scenarios. Such setups can lead to very different conclusions, see e.g. Chapters 6.2-6.3 in [1.
Connections with other Gaussian process priors
The rescaled Brownian motion (rBM) prior we consider here can be thought of as an approximation to the Matérn process of regularity 7 E (0, 1] since they both model 7-smooth functions. The rescaling
factor C, = (n + 1/2) I- in the covariance function k(x, x²) = cn (I A x') plays the same role as the lengthscale parameter for the Matérn and calibrates the GP smoothness. We thus expect similar theoretical results and conclusions to also hold for the Matérn process, which seems to be the case in practice, as we verify numerically in Section 5 Indeed, similar UQ properties to the rBM posterior are conjectured to hold for the full Matérn posterior (39, in particular that the coverage of credible sets satisfies similar qualitative conclusions the three cases in Theorem 3.6
On a more technical level, the success of GPs in nonparametric estimation is known to depend on their sample smoothness, measured through their small ball probability 33
(6]. The posteriors
based on both GPs converge to a - at rate n can 1.2-loss indicating the closeness of their small-ball asymptotics and hence that both GPs distribute prior mass similarly over their supports. This gives a more quantitative notion of the similarities of these GPs.
33 for rBM, K for Matérn),
Theoretical guarantees do exist for the Matérn SGPR. Theorem 3.5 of 32 nicely establishes that in a bandit-like setting and if the true fo is in the prior reproducing kernel Hilbert space H (roughly 7 + 1/2), then inflating credible sets by factor depending on the rank gap and an upper bound for fl ll guarantees coverage. Our parallel results suggest that this inflation is not necessary. since SGPR credible sets in this case will already be conservative. However, since our proofs rely on explicit expressions available for rBM but not for the Matérn, seems difficult to extend our approach to the different setting considered in 32 (and vice-versa); other techniques are required. Our results can thus be viewed as a guide as to what to expect when using a Matérn SGPR and provide different perspective reinforcing the main messages of [32] namely that SGPRs can provide reliable, if conservative, UQ
On the other hand, the squared exponential kernel seems to behave qualitatively differently in numerical simulations, giving different coverage behaviour in the cases considered in Theorem (366 see Section [ This is due to the difference in the smoothness of the sample paths, with squared exponential prior draws being far smoother, in particular analytic. This leads somewhat different UQ behaviour for the true posterior [12). Rescaled Brownian motion is less good approximation for the squared exponential kernel than for the Matérn, and one must thus be cautious about transferring the messages derived here to the squared exponential.
Numerical simulations
We next investigate empirically the applicability of our theoretical results in finite sample sizes and whether the conclusions extend to related settings, such different designs and other GP priors. We consider various nonparametric regression settings of the form II) with C*-smooth truths and 7-smooth priors. We compute both the full posterior (columns marked GP) and SGPR (marked SGPR) with m < n eigenvector inducing variables for (xr) at a fixed point xo E X. and report the root-mean square error (RMSE), negative log predictive density (NLPD), and the length and coverage of the 90%-credibl intervals of the form (10) (i.c. 8 = 0.1), see Section B in the supplement for definitions. The true noise variance in is taken to be o² 1, but is considered unknown and is estimated by maximizing the log- marrinall likelihood as usual. We ran all simulations 500 times and report average values and their standard deviations when relevant. Simulations were Quad-Core Intel Core i5 processor on 2020 Macbook Pro with 16GB of RAM.
2 GHz
Settings of the theoretical results. We consider the setting of our theoretical results in Section with X² [0, 1] and v; = i/(n + 1/2) for = I,.... n. To generate the data, we take fo(x) | - 0.5| x which is exactly a-Hölder zn 0.5, and investigate pointwise inference 20 = 0.5 for
8

---

different choices of a. As Gaussian priors, we consider (i) rescaled Brownian motion of regularity described Section (ii) the Matérn kereel with parameter y and (iii) the square exponential kernel
k(x,x') = expl 2 T x') }² where the lengthscale en - nt 112- is chosen to be appropriate
for estimating 7 HöHöle smooth functions [33]. We take m := n 1 i inducing variables for the SGPR. Results are given in T Fixed Design, n = 1000. We also consider the related case of uniform random design where 24 (0, 1) (Table
Random Design, n = 500).
Prior GP
Coverage
Length
RMSE
NLPD
SGPR
GP
SGPR
GP
SGPR
GP
SGPR
GP
Fixed Design:
n 1000, (a, 7) = (1.0, 0.5)
rBM
0.98
0.98
0.41
0.41
0.09
0.09
-0.90 (0.21)
-0.90 (0.21)
Matérn
0.98
0.98
0.49
0.49
0.10
0.10
-0.68 (0.16)
-0.68 (0.16)
SE
0.91
0.91
0.65
0.65
0.19
0.19
-0.28 (0.30)
-0.28 (0.30)
Fixed Design: n = 1000, (ax 7) = (0.5, 0.5)
rBM
0.74
0.74
0.41
0.41
0.18
0.18
-0.09 (0.35)
-0.09 (0.35)
Matérn
0.84
0.84
0.49
0.49
0.17
0.17
-0.40 (0.32)
-0.40 (0.32)
SE
0.88
0.88
0.65
0.65
0.21
0.21
-0.21 (0.31)
-0.21 (0.31)
Random Design: n
500, (a, 7)
(1.0, 0.5)
rBM
0.98
0.98
0.49 (0.02)
0.49 (0.02)
0.11
0.11
-0.65 (0.15)
-0.65 (0.15)
Matérn
0.96
0.96
0.59 (0.03)
0.59 (0.03)
0.13
0.13
-0.55 (0.14)
-0.55 (0.14)
SE
0.92
0.92
0.76 (0.08)
0.76 (0.08)
0,20
0.20
-0.02 (0.25)
-0.02 (0.25)
Random Design: n = 500, 70, = (0.3, 0.5)
IBM
0.25
0.25
0.49 (0.02)
0.49 (0.02)
0.37
0.37
2.23 (0.86)
2.23 (0.86)
Matérn
0.47
0.47
0.59 (0.02)
0.59 (0.02)
0.34
0.34
0.91 (0.66)
0,91 (0.66)
SE 0.71
0.71
0.77 (0.08)
0.77 (0.08)
0.31
0.31
0.36 (0.52)
0.36 (0.52)
Table 1: Comparison of SGPR and full posterior (marked GP) for 90% pointwise credible intervals for
different values of (a.7). For the SGPR we use m == कमे 24o to 178, 316, 106 and 244 in the four (n, a, 7) regimes above (top-to-bottomm.
inducing variables, corresponding
Theorem 3.6 predicts that for 90% credible sets (8 0.1), m >> m² and Y, the coverage should converge to Po.1 = (O, 1/2)| < 20.9) = 0.98. We see that for n = 1000, the observed coverages for both the full posterior and SGPR based on the Brownian motion prior are very close to this predicted value, and are indeed conservative. In the oversmoothing case < ", coverage is poor and far below the nominal level, while the correct smoothing case (a = 7), coverage is moderately below the nominal level. We see that the asymptotic theory is applicable for reasonable samples sizes.
The Matérn process behaves qualitatively similarly to Brownian motion in these three cases, though the (still conservative) limiting coverage in the undersmoothing case is predicted to be slightly different (39, as reflected in Table On the other hand, the square exponential behaves differently in all cases and generally has wider sets with this lengthscale choice. The wider credible intervals and larger RMSE suggest greater bias of the posterior mean. While rescaled Brownian motion seems a reasonable guide for the Matérn, it does not seem so for the squared exponential, see Section 4
For all three GP priors, the reported metrics are practically indistinguishable between the SGPR and the full GP, confirming that m >> n 1", 1/2 inducing variables seems sufficient to obtain virtually the same pointwise credible intervals as the full posterior. However, since m < 1-22 for any a > O, the KL-divergence between the SGPR and the full posterior grows to infinity as n -> oo in the Brownian motion case (Lemma (0.)maaA.6)..am and hence the SGPR will differ from the posteriors in some ways. This reflects that one does not necessarily need to fully approximate the posterior in order to match its pointwise inference properties. Note that for fixed design, the length of the credible sets is not random since it depends only on the features I1.... .w, and so no standard errors are reported.
Other settings. We turn to simulation settings not covered by our theory. Given the above discussion, we consider only the Matérn process in dimension d = 10 and the corresponding SGPR with v my n2 गाे., As features, consider random design with (i) x; ~ ((([1/2, and (ii) X, 1 Na(0, E,) for [E for j and [E,],: = p otherwise. To better reflect correlation in real-world data, we also consider semi-synthetic dataset with real features but simulated responses.
9

---

As features, we took the first d = 10 columns and randomly sampled 2000 rows from Korean meteorological dataset [ In all cases, we consider estimating the c--öllerr function fo(r) = ||= Toll% at zo 0. and the resulting pointwise 90%-credible sets at Zo coming from the Matérn SGPR. Results are presented in 2 where we see that for m mi the performance of the pointwise credible sets from the SGPR matches that of the full posterior in each design. Due to the higher dimension, credible intervals are less confident (wider) and have larger coverage.
Multidimensional Random Design, n 1000
Design
Smoothness
Coverage
Length
NLPD
Type
P
a
y
SGPR
GP
SGPR
GP
SGPR
GP
Uniform
0.5
0.5
0.96
0.96
1.46 (0.03)
1.46 (0.03)
0.96 (0.26)
0.96 (0.26)
Gaussian
0.0
0.7
0.5
0.99
0.99
2.17 (0.05)
2.17 (0.05)
1.34 (0.11)
1.34 (0.11)
Gaussian
0.2
0,9
1.00
1.00
2.13 (0.06)
2.13 (0.06)
1.90 (0.14)
1.90 (0.14)
Gaussian
0.5
1.1
0.5
1.00
1.00
2.03 (0.06)
2.03 (0.06)
1.36 (0.12)
1.36 (0.12)
Semi-synthetic Data, n = 2000
Smoothness
Coverage
Length
RMSE
NLPD
Cx
Y
SGPR
GP
SGPR
GP
SGPR
GP
SGPR
GP
1.0
0.5
1.00
1.00
2.29
2.29
0.16
0.16
3.22 (0.47) 3.22 (0.47)
Table 2: Comparison of SGPR and full posterior (GP) with Matérn prior for 90% pointwise credible
intervals for different values of 7). The SGPR uses 534 and 1002 inducing variables in the first four rows and fifth row, respectively.
Discussion
In this work, we established the frequentist asymptotic behaviour of pointwise credible intervals coming from sparse variational Gaussian process regression (SGPR) with cigenvector inducing variables based on rescaled Brownian motion prior. We showed that if the prior undersmooths the truth and with enough inducing variables, SGPR credible sets can provide reliable, though conservative, uncertainty quantification, with the coverage converging to a value strictly between the nominal level and one. If the prior oversmooths the truth, UQ can be poor. We further showed that it is not necessary for the SGPR to converge to the true posterior in KL-divergence to have similar behaviour of the credible sets. Our results suggest that properly calibrated SGPRs can perform reliable UQ We verified these conclusions in simulations and discussed connections with other GPs.
Despite the widespread use of SGPR, there are still relatively few theoretical guarantees for their use, particularly for UQ Our work provides some new results in this direction, but most of this research area is still wide open. A key step would be to prove similar results for the most commonly used GP priors, notably the Matérn and squared exponential. Similarly, one would like extend these results to other choices of inducing variables, for instance the eigenfunctions of the kernel operator. While our results give some intuition of what one can expect, they rely on careful analysis of Brownian motion with fixed design, and new ideas and techniques will be needed for these other settings.
It is also unclear what minimal number of inducing variables is needed get good pointwise UQ When a = 7, the threshold for minimax convergence rates for estimation in L is m 22), which is smaller than our bound n 1+2. (?+2), see Section B in the supplement for some related numerical simulations. Our results are also non-adaptive since they assume fixed prior smoothness y, whereas one often selects y in a data-driven way, for instance by maximizing the marginal likelihood or evidence lower bound (ELBO). However, adaptive UQ is a subtle topic which typically requires further assumptions (e.g. self self-similarrit) on the true generative function to even be possible [30]. and extending our results to such settings will require significant technical and conceptual work.
Acknowledgements. We are grateful to four anonymous reviewers for helpful comments that improved the manuscript. The authors would like to thank the Imperial College London-CNRS, PhD Joint Programme for funding which supported Luke Travis in his studies.
10

---

References
[1] M. Abadi, A. Agarwal, Barham, Brevdo, Chen, Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp. G. Irving, M. Isard, Y. Jia, R. Joze- fowicz, L. Kaiser, M. Kudlur, Levenberg. D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Va- sudevan, F. Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. Large-scale machine learning on heterogeneous distributed systems, 2016.
Tensorflow:
[2] P. Alquier and J. Ridgway. Concentration of tempered posteriors and of their variational approximations. Ann. Statist., 480)):1455-4477, 2020.
[3] C. M. Bishop. Pattern recognition and machine learning. Information Science and Statistics. Springer, New York, 2006.
[4] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: a review for statisticians.
J. Amer. Statist. Assoc., 12(518):859-5777 877, 2017.
[5] D. Burt, C. E. Rasmussen, and M. Van Der Wilk. Rates of convergence for sparse variational Gaussian process regression. In Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 862-871. PMLR, 09-15 Jun 2019.
[6] D. R. Burt, C. E. Rasmussen, and M. van der Wilk. Convergence of sparse variational inference Gaussian processes regression. Mach. Learn. Res., 21(131):1 1-63, 1-63, 2020.
[7] Castillo and Nickl. Nonparametric Bernstein-von Mises theorems in Gaussian white noise. Ann. Statist., 1194):199 199992028,,, 2013.
[8] D. Cho, C. Yoo, J. Im, and D.-H. Cha. Comparative assessment of various machine learning- based bias correction methods for numerical weather prediction model forecasts of extreme air temperatures in urban areas. Earth and Space Science, 7, April 2020. doi: 10.1029/ 2019EA000740.
[9] A. G. de G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, Lell-Villagrá, Z. Ghahramani, and J. Hensman. Gpflow: A Gaussian process library using TensorFlow. J. Mach. Learn, Res., 18(40):1-6, 2017.
[10] D. Freedman. On the Bernstein-von Mises theorem with infinite-dimensional parameters. Ann. Statist., 27(4):
1999.
[11] S. Ghosal and A. van der Vaart. Fundamentals of nomparametric Bayesian inference, volume 44 of Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2017.
[12] A. Hadji and B. Szabó. Can we trust Bayesian uncertainty quantification from Gaussian process priors with squared exponential covariance kerne?? SIAM/ASA J. Uncertain. Quantif., 9(1): 185-230, 2021.
[13] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, UAI 2013, Bellevue, WA, USA, August I1-15, 2013. AUAI Press, 2013.
[14] J. Hensman, N. Durrande, and A. Solin. Variational Fourier features for Gaussian processes. J. Mach. Learn. Res., 18(151):1-52, 2018.
[15] M. Hoffmann, J. Rousseau, and J. Schmidt-Hieber. On adaptive posterior concentration rates. Ann. Statist., 43(5):2259-222) 2299, 2015.
[16] I. M. Johnstone. High dimensional Bernstein-von Mises: simple examples. In Borrowing strength: theory powering applications- Festschrift for Lawrence D. Brown, volume 6 of Inst. Math. Stat. (IMS) Collect., pages 87-98. Inst. Math. Statist., Beachwood, OH, 2010.
[17] B. T. Knapik, A. W. van der Vaart, and J. H. van Zanten. Bayesian inverse problems with Gaussian priors. Ann. Statist., 39(5):2626-2655 2011.

---

[18] C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Research Nat. Bur. Standards, 45:255-282, 1950.
[19] M. Lázaro-Gredilla and A. Figueiras-Videe Inter-domain Gaussian processes for sparse infer- ence using inducing features. In Advances in Nearal Information Processing Systems, volume 22.
Curran Associates, Inc., 2009.
[20] H. Liu, Y.-S. Ong, X. Shen, and J. Cai. When Gaussian process meets big data: review of scalable GPs. IEEE Trans. Neural Netw. Learn. Syst. 31(11)=445544423, 2020.
[21] A. G. d. G. Matthews, J. Hensman, R. Turner, and Z. Ghahramani. On sparse variational methods and the Kullback-Leible divergence between stochastic processes. In Proceedings of the 19th International Conference on Artificial Intelligence and Staristics, volume 51 of Proceedings of Machine Learning Research, pages 231-239, Cadiz, Spain, 09-11 May 2016. PMLR.
[22] D. Nieman, B. Szabo, and H. van Zanten. Contraction rates for sparse variational approximations in Gaussian process regression. J. Mach. Learn. Res., 23(205): 2022.
[23] D. Nieman, B. Szabo, and H. van Zanten. Uncertainty quantification for sparse spectral variational approximations in Gaussian process regression. Electron. Stat., 17(2) 2288, 2023.
[24] D. Pati, A. Bhattacharya, and Y. Yang. On statistical optimality of variational Bayes. In Proceedings the Twenty-Fie First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1579-1588. PMLR, 09-11 11 Apr 2018.
[25] C. E. Rasmussen and C. K. I. Williams. Gaussian processes for machine learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2006.
[26] K. Ray and B. Szabó. Variational Bayes for high-dimensional linear regression with sparse priors. J. Amer. Statist, Assoc., 12755))) 2022.
[27] K. Ray, B. Szabo, and G. Clara. Spike and slab variational Bayes for high dimensional logistic regression. In Advances in Neural Information Processing Systems, volume 33, pages 14423-14434, Curran Associates, Inc., 2020.
[28] S. Rossi, M. Heinonen, E. Bonilla, Z. Shen, and M. Filippone. Sparse Gaussian processes revisited: Bayesian approaches to inducing-varian approximations. In Proceedings of The 24th International Conference Artificial Intelligence and Statistics, volume 130 of Proceedings of Machine Learning Research, pages 1845. PMLR, 13-15 Apr 2021.
[29] S. Sniekers and A. van der Vaart. Credible sets in the fixed design model with Brownian motion prior. Statist. Plann. Inference, 166:78 0077788862776 2015.
[30] B. Szabó, A. W. van der Vaart, and J. H. van Zanten. Frequentist coverage of adaptive nonparametric Bayesian credible sets. Ann. Statist., 43(4):199111428, 2015.
[31] M. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Proceed- ings of the Twelfth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 567-574, PMLR, 16-18 Apr 2009.
[32] S. Vakili, J. Scarlett, D.-S. Shiu, and A. Bernacchia. Improved convergence rates for sparse approximation methods in kernel-based learning. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 21960 21983. PMLR, 17-23 Jul 2022.
[33] A. van der Vaart and H. van Zanten. Bayesian inference with rescaled Gaussian process priors. Electron. J. Stat., 1:433 448, 2007.
[34] A. van der Vaart and H. van Zanten. Information rates of nonparametric Gaussian process methods. J. Mach. Learn. Res., 22060))2095-2119 2011.
12

---

[35] A. W. van der Vaart. Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge. 1998.
[36] A. W. van der Vaart and J. H. van Zanten. Rates of contraction of posterior distributions based 2008. ISSN 0090-5364.
on Gaussian process priors. Ann. Statist,,
[37] Y. Wang and D. M. Blei. Frequentist consistency of variational Bayes. J. Amer. Statist. Assoc.,
1 114(527)11147-1161, 2019.
[38] V. Wild, M. Kanagawa, and D. Sejdinovic. Connections and equivalences between the nyström method and sparse variational Gaussian processes. arXiv:2106.0 012.1100.111111 01121, June 2021.
[39] Y. Yang, Bhattacharya, and Pati. Frequentist coverage and sup-norm convergence rate in Gaussian process regression, art. 1708.04753, 2017.
[40] W. W. Yoo and S. Ghosal. Supremum norm posterior contraction and credible sets for nonpara- metric multivariate regression. Ann. Statist., 1469-1106-, 2016.
[41] Zhang and C. Gao. Convergence rates variational posterior distributions. Ann. Statist., 48 (4):2180-2200 2020.
13

---

Supplementary material to "Pointwise uncertainty quantification for sparse variational Gaussian process regression with Brownian motion prior"
A Proofs
Throughout this section, for sequences on and b, we will write On S by if there exists > 0 such that on I Cb, for large enough, and will write an x b,, if both an S bn, and 6n S On. Recall that for simplicity, we take the noise variance o² 1 in the regression model (T).
A.1 General results on pointwise inference
Proof of Lemma 3.7 Using the definitions and (3) of the posterior and SGPR means, have Vm (2) = (x) (Tm (x), from which follows on taking the Es -xxpettaatn Next, recall
that Un (2) = (x))
Σラ= nov,v?]
y. and thus
Σ 1/ Viv
(x)x) = Var(Ux(2)) =
kn
No
kn(x) + +mm(I)
12
m
(x)z)
1m (2)
(z) + I' (=))3
from which (8) follows. Finally,
(x) kn, xx, k( kx)
Σ No
kx
k(x,x) (x)
Σ VJV
kn - (n))
Σ -=++1
Nyvsyy kn(z),
which gives
Proof of Proposition 3.2 We will prove this for the full posterior, with the expressions for the SGPR following the same. Denoting by Eny the expectation with respect to the posterior, Markov's inequality yields,
- -(((()))f(()) MYE, Eny (f(x) fo(x)))
Mas2
Moreover,
Eny (f(x)
- fo(x))) (vn(x)) 5o(x))) + Varnj (f(x)) = -o(x)) + +x)())
Now,
Ejo Emy (f(x) - fo(x))) = Es, - (2)(z) + +x)(x) + (z))).
Thus for []) +f) (x) and any Mn xo, (1f(2) fo(x) MnsnlY)
Proof of Proposition 3.4 We have,
(f(fo(x) C.h
- fo(x)] (1-50((z)) P1o([mml S 21-6),
14

---

for Vm N
bm(
t (x:) on (tt)
Now. if (Tm (2), fo) (bb we have 6m()) (1+
(m(
we have f (x) om (z)
1 (a) and thus,
and further, by expressions
and
o (w))
inf Pa. (|vml S 61 2 lim pso ( ml $ 21-6) n-000
lim inf Pfo (fo(z) E C.)) n-000
n-000
in addition, we have ||rm (x)] 0 (x), and (x); kn (x)) = a( [ (2)) then by expressions
g and p we have en (z) (=) o(1) and (I) (z) + a()) and thus,
lim Pio ([Vmll S =1 lim Ps ([V,) =1- lim Pso (Jo(2) E C.))
A.2 Rescaled Brownian motion prior
Recall that we take regularly spaced design points 2; n+11/2 .,..,, n. and that consider the rescaled Brownian motion prior, which mean-zero Gaussian process with covariance kernel k(r, x') = (x/Ay) C. min(2, y) with 6x (m+1/2) FT for > 0. For notational convenience, we write N (n + 1/2)/c,, (n 1/2) 1.27 The kernel matrix K.n evaluated at the sample points then equals
1/ 1/(n+ + 1/2)
+ 1/2)
//n+1/22) 2/(n 1/2) 2/(n + 1/2)
1/(n + 1/2)
1
1
2
2
K
N
1/(n++ 1/2) 2/(n 1/2)
n/(n + 1/2)/
2
For this specific choice of prior and design, we can compute fairly explicit expressions for the eigenvalues and eigenvectors of K which will tum allow to obtain the precise pointwise asymptotics of the VB method needed for our results.
Lemma A.1 (Eigenvalues and eigenvectors). The eigenvalues (11) of the covariance matrix K given by
j 3-11 -1/2 bh n 1/2 n with
2N(1 cos(();)) for j 1.....7. n. Moreover, their size satisfies
Mg x
No3
x (n 1/2)cn (j 11222²
The corresponding orthonormal eigenvectors Vi...., Vn where v, (u),...,,?)) given by
2 (in((y))
o/
v 2n +1
Proof. Let be an eigenvector of K with corresponding eigenvalue μe, so that K pec
Applying elementary row operations to this equation, we arrive the following equivalent equations:
v
2v - 2
20 - v - 3
1
nn
20 2 2
n-1
From this, one can see that to obtain a non-trivial eigenvector v of O, we require v o O. We thus consider the unnormalized eigenvector w (w ...., w" E R", for which we w 1 for simplicity. Further define w" = 0 and w = -1, so that the first (n - 1) equations in the last display are equivalent (after rearranging) to the recurrence relation
u = (2
1
w i -
j j=
(11)
uN
15

---

while the last equation gives the boundary condition
JUN
W
(12)
The linear recurrence relation (TT) has solution
Φ+ - -- w
004 -b).
(13)
for
4
++ 1
J.N =: x bt,
N
2
N2N2
with x
2pN and y
V JN N (computation of the eigenvalues below reveals that
4 0 for all eigenvalues J.) Note that 13) provides the solution 11 for any arbitrary
N= pN
value of μ., but that (12) will later be needed ensure that μ is an eigenvalue of Knn- Since
x² + y² = 1, we can conveniently rewrite this in complex exponential form as
Φ +ip
where
COS va Re((1) = x =
(14)
2μtN'
and db is an alternative parametrization of u.. This form makes it easy to compute expressions for oi and hence evaluate (13). Indeed, we have
1
④ (0) x²)= sin(
w
(15)
②+
sin t0
giving the form of the unnormalized eigenvector w = (1, w² ..... w²) corresponding to an eigen- value μ via (T... To normalize the eigenvector, using that cos (2je) = 1 - 2 sin (j+))
sin sin (j+b
1 =
(w))2 1=1
2 (1 cos(()t))
sin
sin V0 j=1
2n + sin((2n 1)
4sin²
eb
sin 1b
where in the last line we have used the Lagrange trigonometric identity that Σj= cos(22 (2jz))
sln(2n+1)w) for any vi + +(mmod т)) (else equals n + 1). On computing the values & corresponding 2sin
to the eigenvalues in (77 below, one has sin + 1)x²) = 2 sin((n + 1/2)))) cos((n+112) 1/)
+
= 0.
and hence ||wl| Tain
Thus the normalized eigenvector with w/ //|lw||2 has the required
coordinates.
2sin(j
(16)
22n +
Turning to the computation of the eigenvalues, the boundary condition
12 defines whether is an
eigenvalue. Substituting the solution (15) into [12,, any that satisfies
sin(nd)
WN sin((n 1)t)
sin 2/ WN - sin
equivalently
(aN 1) sin(nop) = nN sin((n 1) and sin + 0,
is an eigenvalue of K Indeed, if one approaches the problem in the traditional way, by computing the determinant of (Kn uII,)) (see Lemma A.2 below), then one has
(-1)a
1-1
K nm μlpp
110
sin(nu)) n1" sin((n 1)
sin
N
sin(nep) JNN sin((n 1)4b) and sin # 0
N
(μN 1) sin(no))
WN sin((n 1)6b 0 and
sin zb * 0,
16

---

exactly as above. Recalling that db and are related by cos = -1- 1 equivalent to
this last expression is 2aN
1 HN sin(n - sin((n 1)15)
(2 1b 1) sin sin((n 1)xb)
sin((n 1))) sin(n
22in(4//2) cos((n + 1/2)x));
where we have used the identity cos Asin B = sin(A - B) sin(A B) twice. By (4). the eigenvalues j are only defined by cos v and hence we may restrict & to [- σ) and look solutions within this range. For & in this range, the condition sin ei * 0 then implies sin( (2/) # 0. and hence we need only look for solutions to the equation cos((n + 1/2)r/ = O. For & E [-,, TT, these are given by
j- n + 1/2
- -
(17)
By symmetry, cos 1 cos 11-; for 1, ,...., n and hence we may discard the roots
-n + 1...., 0. Moreover, sin th-n 0 violating the restriction sin & o O. This leaves us with the eigenvalues
Mj = 2N(1 COs v)))
= i=1,...,7
as desired.
It remains only to quantify the size of the eigenvalues py. Expanding cos 1 x²/21 + 24/41-
o(**), ), we have x² x²/ /12 < 2(1 cos cos I x² for all x, and thus
x²/12-
1
1
Nx
2N cos
Nx x²/12)
Since 1/y
x+1/2" n we have 1 < 5/6 (on bounding n² by 10 and the fraction by 1) and so (1 01-01)
j-1/2
No. Hence we have Noz S Mj S Noy and My x Nv.
Let |A| denote the determinant of a matrix Lemma A.2. The solutions to the characteristic equation of K Knn p.ll O, given by the solutions to
(-1)a
(11- (100 N(-11 ) sin(mup) 1²" "inssin((n - 1)(b) = 0,
sin 1b
where cos 1b
2NK
Proof. Using row operations, the desired determinant equals
1
N
N
AY
N
p
N
N
K nnn
pln
N
N
0
:
11
0
11
N
which can be further be reduced to
0
0
2
NN
K
paaal=
=: das
O
0
N
11
17

---

the determinant of a symmetric tridiagonal matrix. Using the Laplace expansion of the determinant terms of the determinants of its (n 1) x (n submatrices, we obtain the recurrence relation dn, ( pi) d - - - The solution this recurrence is given by
1
dp
(o7 1
Φ+ Φ_
where #+ satisfy x² - (3 1)) Φ 2²² 0. That is,
(2 Φ
2
Similar to when we solved the recurrence in the proof of Lemma A.1 we end up with Φ = μ ip where cos(el") 2Nz 1. We the notation v/' to emphasise that these angles are different to the etss found the proof of Lemma A.1 indeed, is easy to see that e/ 1b + n, so that sin(ny) (nx)) = (-1)" sin(re)). This yields
1
dn = (-1)
sin((n + 1)v),
sin(
and the result is shown by seeing
1/4
1
Knn - 11|
N 14
dn-1 - 22
1 sin(me)) (111)=-2 1=m²
= (-1)m-1
N 11)
sin((n 1)b)
sin(
sin(v)
1-1)
n* N 1111 sin(ne)) - 11" sin((n 100--))0. 1)1²)
sin x
Given the explicit expressions derived for the eigenvalues and eigenvectors of the covariance matrix K we derive bounds for the size of the remainder terms presented in Lemma 3.1
Lemma A.3. Suppose that fo E CO 1], E (0, 1I, and consider the rescaled Brownian motion prior of regularity 7 E (0. The corresponding SGPR Q" based on the first m mn + co eigenvector inducing variables satisfies, as + xc,
(m) = bx(x)) + (ne,m -1+aa) (x) = (x))) + O(n,, -3) o (x) = (x)) + O 3))
(18)
(19)
(20)
Proof. First, suppose that x x; a design point. Starting with the bias, we have,
bm(x)= = bn(x)) + (I)TT
Σ
fo = bn (x)
Σ
(x)
j=m+1
j=++1
In this case kn (I,)) v, = p.y of and the next step is to bound vř fo for each j. Using the explicit form for v, given in Lemma A.I
''v,
2
2 sin ((6 (1111/2)
1/2)T n 1/2
T v fo
sin (in(1)))ffoz))
fo(z))
/2n +1
/2n +
1=1
1=1
2
sin (z) - 1/2) fo(x1))
2n
18

---

Now, writing Ty = (j 1/2)JT, we have sin() = (e -e- and thus
sin (zi(j 1/2)m) fo(ai)
x
sin(t(j 1/2) + +n
e fo(x)d S'" e fo(a))d + inlaa
2i
nem n 1- x nī;
for fo E (0, where the first inequality follows from Lemma A.4 below. We have thus shown
that [v3 ful Vn(j 1/2) -a Using this, that M (n + 1/2)cn 1/2)22 by Lemma A.1 and that
My S I,
(1) j=m+1
E
1/2)) nen
Σ
j-2-a m (1+a)
j=m+1
required. Next, we have,
(x) - (1(()) - I|'m (z)) (x()) kn (x))
Σ (vyjzjk(())
jmm+1
Using again that Hj and so Kn (v;)) vv; = H.D have,
n
n²c2
(x)²
Σ
nivovz kn(z)= Σ T
Σ
(j- 1/21/ S ne,
j=m+1
-=
-=++1
as required. Finally,
n Σ
o (x) 0x + (x)²)
ny kn(z),
(21)
j=m+1
and again at design point =
n
kn (x)) x
Σ jmm+1
kn(x) = Σ nnz(t()) S negm-
jum+1
This establishes the three relations at the design points Mi.
Turning the general case E [0, 1]. let í be such that vi < 31+1 Recalling that the covariance
function equals k(z, x²) (x A x'), we obtain
x1
0
(x) = Cn
(x;) + Cn 32g
where 0 < x - 24 <
1 Consider first the case of the SGPR variance (x). The last display
n+1/2
implics that (k,(x) - (x:)) = cn(x xx) Eisi1 u/. and hence the quadratic term in (2) can be rewritten
k, (
Σ f=m+1
kn(x) = kT()))
Σ no
no
kx()))
f=m+1
+ 20n
(x - xz)
v. kn
(2x)
(22)
j=m+1
l=++1
+
T
[
xx
j=
i=++1
19

---

The first term is exactly the quadratic term evaluated at the design point which we showed above
has size On m 3). Using Lemmas A.1 and A.4
2
sin( - 1/2)=)
2n
1mit1
1mi+1
2n
(L/ sin(t(j 1/2) +0(n-1)
22n +1
2n
cos(1;;j - 1/2)m²
+ +0(n-112)=0
vn 1-122 + +O(n +(nn11/2).
V2n +1
(j - 1-/2
By Lemma A.1 of k, (24) 1122} -((av/nj-2)) Since also N]j S 1, the second term in bounded by a multiple of
22
j kz)) &
Σ j=m+1
S
Onncm
-=
The third term in (22) can similarly be shown be OC O We have thus shown that (22) is of size m -3) exactly as in the case where a was design point, which implies the desired bound for the SGPR variance o (z)).
The case e, (z) follows in a similar fashion, simply replacing ni by nī. while the bias 6m (z) uses additionally the bound vv]"fol fol K vn(j 1/2) derived above. In conclusion. the three quantities bm (x), e (z) and of (I) satisfy the same bounds when is a design point, which completes the proof.
The next standard bound controls the error when approximating C function by its Riemann sum. Lemma A.4. If E C 1] with E (0, 1]. then
S f(t)dt *E(i)) S (n 1/2)
with 2, 2+1/2 for i I....,
The following lemma shows that for large enough, the bias, frequentist variance of the SGPR mean and SGPR variance can each be related to those of the full posterior derived in (29
Lemma A.5. If m >> n n=myy (==))
then the SGPR satisfies
(x) by ((())
(x()) = (1(z)((++((1))
o. (x) = +=,(=)(+++([)
Proof. From 29 we have the following expressions for the full posterior:
[xx()) S 1627
(x()) x n Li2.
o (z)
X n 21
Thus, if we can show that the O terms given in equations 18
20 are smaller than their respective
terms above, we have the result. By Lemma A.3 we have bm (z)
bn + ((+o)))
and mC, (1+a) - n m m-(1mm))) (1+a) S r if and only if n 1 (Hte) m. Thus, if
1427 (1ta) then ncC, (1+aa) = o(bx(z)) and (x) = (z) + o(1)). Similarly,
nc m -3 < (1) (x) x x - m >> पमे.
Thus, (x) (1) + o(1) and (x) = (z)( +01) m > n 212) (in fact these are implied even at the smaller bound n n.)5)
Proof of Theorem 3.5 The proof follows from Proposition 3.2 combined with the two groups of expressions 7 and [8 (00).
20

---

Proof of Theorem 3.6 The coverage equals
(f) E ) = Pso (1)m(r) - fo(2)) S 21-50m(z)) = P/o | I 21-6),
for Vm N Um (w) L. (m)
Gin (a) Now, since m > n 11424 (ttt) we have (m(
o(1))
om (2))
en (z) e (1) + ((1) and o (I) of (1) + +L)) and thus,
lim Pso ((vml S =1-6)= lim Pso (|vnl S =1-)) lim Pso (Jo(z) C.). n-000
n=000
n-000
Lemma A.6. The Kullback-Leibeee divergence between the SGPR Q' and the full posterior satisfies
j=m+1 Mj + V (
2KLQQ'I1(X))) - yy Σ Mj
Σ log o2 p3
Pj.
j=m+1
j=m+1
Suppose further that fo C a E (0, 11, and consider the fixed design setting with the rescaled Brown motion prior of regularity E (0, If m > ni then
Eg LQQ X n m-
in 2:2: m 1-20
In particular, nt <m m < n niTz then E KLQQ" [I(()))) + 00.
Proof. We have from 31
KI(QIIT(-IT)) (y" 1Q -0=-) + log
+ σ² tr((K, - Qn)):
Kn,
where K,, Knn + 121, Qn Qnn + 0221, Knn ΣE= VJV and Qun E
We thus have
Q_ K
1 0² 0² 02+ ) Vjv
Σ
P4 v,vY,
0² +
j=m+1
j=m+1
Q
K. j= o² : and tr(K,, - Qn) = Ej=m++ Hý, giving the results. We now consider the
three terms given in the first part of the lemma. Defining M == Ej=m+1 uyto V.y we have
E yT
My + 02 Viv3
Hj
y
1 σ²
Eso (fo + +e) + E)
02 j=m+1
1
E -
fo Mfo +
Me]
o²
Recalling from the proof of Lemma A.3 that VV S
vnj
we can control
1 Oo Mr
Σ
Mj
1o)
2
x
nuyj
3-22 -1-2a n
1- m
j=m+1
j=m+1
Next, we have
E= Me] = Es tr("" Me) = E (M(E")) = tr (,ME = (r(M) x m
For the second term in the first part of the lemma.
log 02 +2+13 14 =
Σ
log
1 +
H4g
Σ
+
0
02
o²
j=
j=m+1
j=m+1
0
j=m+1
j=m+1
/mm
O m 7)
Hence, adding these expressions the third term in the first part of the lemma yields
E ([Q x n + n n m 10 On m 3)
3427
The last term is of strictly smaller order than the first two if m > ामन.
21

---

B Additional numerical simulations
We provide here additional simulations to those in Section 5 investigating further the minimal number of inducing variables needed for good pointwise UQ and the sensitivity of the results to misspecification of the noise distribution. To illustrate the pointwise approach to UQ using full posterior inference and SGPRs, one refer to Figure T
in the paper.
Threshold for the number of inducing variables. We showed in Sections and 5 that for m >>
1+2y (145) inducing variables, the SGPR behaves similarly to the true posterior for pointwise inference. When a = 7; the threshold for minimax convergence rates for estimation in L is the
smaller bound m >> T+20 [2 We next investigate how the SGPR behaves near this threshold.
Return the fixed (X, i/(n 1/2)) and uniform random (X; rid U0, 1)) design settings on [0, 1] and consider the same Gaussian processes as before: rescaled Brownian motion, Matérn and rescaled square exponential. We consider three different values of m given by m = n ािशन / log n,
m4 11+22 log n, so that these values of m are just below and above the threshold 1+27 respectively, and the full posterior case m n2. We consider the 1 Hölderr function fo(x) = |22 1/2
and 2 Hölder function fo( = sign(r x 1/2) 1/2
Coverage
Length
RMSE
NLPD
GP
m
m+
n
m=
m
m=
mf
m=
m4
Fixed Design:
n
1000,
(ax, 7)
= (1, 0.5)
rBM
1.00
0.98
0.98
0.52
0.42
0.42
0.06
0.09
0.09
-0.85
-0.89
-0.89
Matérn
1.00
0.98
0.98
0.69
0.50
0.50
0.07
0.11
0.11
-0.59
-0.70
-0.70
SE
1.00
0.93
0.93
2.73
0.66
0.66
0.08
0.18
0.18
0.74
-0.28
-0.28
Fixed Design:
n
1000,
(a,7)
= =2,
1.5)
rBM
0.98
0.98
0.98
0.18
0.17
0.17
0.04
0.04
0.04
-1.72
-1.74
-1.74
Matérn
1.00
0.95
0.95
0.56
0.22
0.22
0.04
0.05
0.05
-0.82
-1.51
-1.51
SE
1.00
0.91
0.91
2.29
0.31
0.31
0.05
0.09
0.09
0.56
-0.99
-0.99
Random Design: n
500, (a, 7)
=
0.,)
rBM
1.00
0.97
0.97
0.59
0.50
0.50
0.08
0.10
0.10
-0.70
-0.75
-0.75
Matérn
1.00
0.96
0.96
0.87
0.59
0.59
0.07
0.12
0.12
-0.38
-0.58
-0.58
SE
1.00
0.94
0.94
2.70
0.77
0.77
0.10
0.21
0.21
0.73
-0.13
-0.13
Random Design: n
500, (a, 7)
(2,
1.5)
rBM
0.98
0.95
0.95
0.23
0.23
0.23
0.05
0.05
0.05
-1.49
-1.49
-1.49
Matérn
1.00
0.96
0.96
0.57
0.29
0.29
0.05
0.07
0.07
-0.79
-1.19
-1.19
SE
1.00
0.95
0.95
2.10
0.39
0.39
0.05
0.10
0.10
0.47
-0.86
-0.86
Table 3:
Comparison of SGPR and full posterior for 90% pointwise credible intervals for different values of (n, a; 7) in dimension d 1. The column headers describe how many inducing variables
were used, with m := n n / log(n)) and m == n n2-22 log(n) and the full posterior represented by n. The hypothesised bound n ामे. equals 32, 6, 22 and 5 in the four different regimes (top-to- bottom).
The results of our experiments for different values of (n, 0, 7) are given in Table we note here that we exclude estimates of the standard error in this Table due to space constraints, but they are comparable to those presented in Table the main text. In all cases, one can see that the coverage, length and RMSE are practically indistinguishable between the cases m m and the full posterior
n. suggesting that the choice of m >> 1422 inducing variables is sufficient to achieve the same performance as the full posterior for pointwise inference. Furthermore, for fixed design, we again see the Brownian motion with both choices = m and m achieves the expected coverage of 0.98 predicted by our theory.
Conversely, for m = m inducing variables, the coverage and credible interval length was larger in almost all cases, sometimes significantly so. This suggests that below this threshold, the SGPR is no longer good approximation for the true posterior, matching the findings in 22 23 for 1,2-type inference. On the other hand, the RMSE is usually reduced, indicating that additional smoothing helps the sparse posterior mean for estimation, which is not surprising since in our simulation choices
22

---

the prior undersmooths the truth. The main message is that taking fewer inducing points than this threshold leads a significant increase in the credible interval width, making the UQ less informative and overly conservative. It is interesting that this threshold shows so clearly, since m and m² differ only by a subpolynomial factor.
The same story holds for random design, with some additional fluctuations due to this additional source of randomness. We recall that for fixed design, the credible intervals have deterministic length, whereas for random design these are also now random.
It would be interesting to extend our theoretical results to the threshold n निक as suggested by these simulations. However, the mismatch between pointwise loss and the I.2-distance induced by the likelihood makes using standard tools from Bayesian nonparametrics difficult, see [15] for further discussion.
Multidimennional setting. Consider again d = 10 dimensional covariates with the Matérn kernel
of smoothness 1/2, and the corresponding threshold मन needed for the minimax contraction rate in L [22]. We take uniform random design X; لان 4[111/2, 1/2]*), true function fo (z) = ||x]||
for a 0.9 and investigate the pointwise 90% intervals the point zo (0, ................................ , 0). We
now take three values of m given by m net2./ log n. m = n nd42y log n and m = n. The results are presented in Table Indeed, one sees that the variational posterior performs identically to the full posterior once at least m2+ = n2((d+))) log n >> n क*2द inducing variables have been used, while m inducing variables do not seem be enough to get the same performance.
Coverage
Length
RMSE
NLPD
m m2+
m
m
n
my
m
12
m
m
π
0.93
1.00 1.00
2.59
2.40
2.40
0.91
0.73
0.73
1.35
1.10
1.10
Table 4: Comparison of SGPR and full posterior with Matérn prior for 90% pointwise credible
intervals for (d, a, 7) = (10, 0.9, 0.5). The column headers describe how many inducing variables
were used with m²- न,न / log(n) and m n
loo(n) and the full posterior represented
by n. Here, n n442y 284.
Misspecified noise. We now take the same fixed design setting as Section of the main paper, but sample Y fX,) + Ex, where E are independent standard Laplace random variables. Based on Table 5 the results are much the same as in the well-specified case, so that the credible sets demonstrate some robustness noise misspecification. Note that while the coverage remains broadly similar, the credible intervals are slightly wider and have larger RMSE reflecting the larger noise.
Smoothness
Coverage
Length
RMSE
NLPD
GP
7
SGPR
GP
SGPR
GP
SGPR
GP
SGPR
GP
rBM
0.5
0.98
0.98
0.49
0.49
0.11
0.11
-0.71 (0.15)
-0.71 (0.15)
Matérn
1
0.5
0.98
0.98
0.58
0.58
0.12
0.12
-0.58 (0.17)
-0.58 (0.17)
rBM
0.5
0.5
0.69
0.69
0.49
0.49
0.22
0.22
0.11 (0.32)
0.11 (0.32)
Matérn
0.5
0.5
0.88
0.88
0.58
0.58
0.20
0.20
-0.17 (0.20)
-0.17 (0.20)
Table 5: Misspecified noise. SGPR summary statistics for 90% pointwise credible intervals for
different values of 7a,7) with regular fixed design on [0, 1] with = 1000 and Laplace noise. The
SGPR uses m = n 1-27
242
inducing variables, which is 126 and 316 in the two regimes.
B.1 Definitions of simulation statisties
For completeness, we define here the statistics reported in the simulation sections based on M Monte-Carhe samples. For Y-O) (X,) + +,,") i I,..., and j 1...., M, let ((ه) (zo) denote the credible interval for f*o), Fl (mr) and (نا (zg)² the posterior mean and variance from the SGPR with m inducing variables computed from the joh realisation of the data. We then consider the Monte-Carlo estimates of the coverage, credible interval length, root-mean square error (RMSE)
23

---

and negative log-pre rrgepeddittive density (NLPD) given by
Cov ([f(zz)) (zo)),
Σ
RMSE-
2l (xo) fo(2x))22
M
M
j=1
Length Ch (zo)) NLPD Σ
log N (fo(zo)lf (zo), oi (z0)²)
Len
M
M j=1
respectively, where A (fu(zo)lsm)) (zo), om (xr)²) with mean ! (mr)) and variance an (zro)² evaluated at the point fo(x))
denotes the density of the Gaussian distribution
24


Reviews:
[{"review_id": "eHHFHhZFDm", "reviewer": "Reviewer_tkgm", "summary": "This paper considers the problem of uncertainty quantification for SVGP models with a Brownian motion prior for a non-parametric estimation problem. Theoretical results are presented that show ", "questions": "What precisely is Q*. Can you please define it formally in section 2. Perhaps Q* should include a gamma in the subscript. Is Q* the eigenvector SVGP or any SVGP as defined above (3).\n\nSimilarly please also provide the definition of rescaled BM in the preliminaries. \n\nA formula for the Mattern kernel in Section 5 would be useful.\n\nTable 1 columns currently labeled m* and n need to be renamed to something more informative and which doesnt require reading the paragraph above the table before reading the table. The caption can also be more informative (and better spaced from the table). The explicit value of m* would also be useful to know instead of the formula in the caption. A brief description of coverage, length and RMSE is necessary to read the table clearly. Alas, these are relegated to the appendix. (Similarly for Table2)\n\nn,2 subscript in Proposition 3.4 is unclear.\n\nIs Lemma 3.1 a noteworthy result or can it be part of section 2? Similarly the paragraph above Prop 3.2\n\nIn Section 2.1, is it GP(nu0, k) instead of m? please dont use m for mean function since m is the number of inducing points.\n\nThe space in which u_i belong is quite unclear. It appears they are actually part of the target space (ie same as space of y_i). Then the choice of the naming them Eigenvector inducing features is perhaps misleading, and Eigenvector inducing \"targets\" is perhaps more appropriate.\n\nThe term cov(f(x), u_i) is confusing. How is this covariance function defined? since the kernel requires two points in the domain as input. Please clarify this immediately following the definitions of Kxm etc., rather than point the reader to [38].\n\nThe bolding can be more consistent. I think r_m is a vector which is not bolded unlike other vector quantities.\n\nIn Section 3.1 end of line 150, shouldn't this also be f(x)|Dn since SVGP also is conditional on n samples?\n\nThe choice of naming tn as \"frequentist variance\" is unclear. Maybe a comment on that (at least in OpenReview) would be useful, but also to the reader of this paper. As suggested earlier, the discussion leading upto Prop 3.2 can be part of Section 2 since everything follows naturally from the preliminaries and does not require any special mention in the section on Main results.\n\nThe role of Mn in all results is a little unclear. An example Mn and its implications would be useful for the reader.", "limitations": "No limitations envisioned.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The results are fundamental in nature and of broad interest, and clearly explained except in places where quantities are not clearly defined.", "weaknesses": "The main contributions of the paper are a little unclear. It is a bit unclear which results are known and which results are new.\n\nIt is also unclear what is special about the rescaled Brownian motion prior. Does the paper recommend using the rBM prior over other priors because of some special properties. Otherwise the choice of rBM in the title is not well motivated.\n\nDoes the proof technique also extend to more general prior models beyond rBM.\n\nThe choice of rBM prior as the main focus is also not well motivated. If this is a standard choice, this must be explained clearly.\n\nIt would perhaps be useful to provide a brief sketch of the proof. For eg the discussion below Proposition 3.2 is somewhat of a proof sketch but isn't labeled as such. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive suggestions and helpful comments. Our response:\n\n- The main contributions of the paper are a little unclear/What is special about the rescaled Brownian motion prior?/Does the proof technique also extend to more general prior models beyond rBM?\n\nPlease see our general reply.\n\n- It would perhaps be useful to provide a brief sketch of the proof.\n\nWe would like to include this, but unfortunately do not have space in the main text.\n\n- What precisely is Q*?\n\nThis is the SGPR based on $m$ inducing variables. This is defined in Section 2 and we have added reminders at various points in the text to help the reader, e.g. \"the SGPR $Q^*$. We have also written $Q_\\gamma^*$ for our main results for rescaled Brownian motion to more clearly differentiate when our results apply to general GPs, and when to rescaled Brownian motion.\n\n- Similarly please also provide the definition of rescaled BM in the preliminaries. \n\nThis has been added to the start of Section 2.\n\n- Fixing Tables.\n\nWe have followed your suggestions, and hopefully the table is now easier to read with a more informative caption. We have kept the definitions of the summary statistics (RMSE, etc.) in the supplement due to space constraints.\n\n- Is Lemma 3.1 a noteworthy result or can it be part of section 2? \n\nWe expect that most of this is known to experts, but we have not seen it clearly written down. Given it is central to our proof, we have kept it in Section 3. In particular, the part studying the frequentist variance of the SGPR mean may not be known given there is relatively little overlap between those working on SGPR and statistical theory. While not a difficult result, we find it informative.\n\n- The space in which u_i belong is quite unclear. It appears they are actually part of the target space (ie same as space of y_i). Then the choice of the naming them Eigenvector inducing features is perhaps misleading, and Eigenvector inducing \"targets\" is perhaps more appropriate.\n\nThank you for point this out. We changed this to \"eigenvector inducing variables\" to avoid confusion.\n\n- The term cov(f(x), u_i) is confusing. How is this covariance function defined? since the kernel requires two points in the domain as input. Please clarify this immediately following the definitions of Kxm etc., rather than point the reader to [38].\n\nWe have added the definition. In our case, $u_j = \\sum_{i=1}^n v_j^i f(x_i)$ and hence $cov(f(x),u_i) = \\sum_{i=1}^n v_j^i k(x,x_i)$.\n\n- The bolding can be more consistent.\n\nWe have tried to make this more consistent throughtout following your suggestion.\n\n- In Section 3.1 end of line 150, shouldn't this also be f(x)|Dn since SVGP also is conditional on n samples?\n\nWhile true, we prefer not to write the conditioning to clearly differentiate this from the full posterior, which is more traditionally written this way.\n\n- The choice of naming tn as \"frequentist variance\" is unclear. Maybe a comment on that (at least in OpenReview) would be useful, but also to the reader of this paper. \n\nOne can consider the (sparse) posterior mean as a frequentist estimator, i.e. the randonmess comes through the data generating distribution $y \\sim P_f$. Then you can study the variance of this estimator with respect to the data generating distribution. Together with the bias of the posterior mean, this measures how far this posterior centering will be on average from the truth, and hence is a key quantity to understand whether the credible intervals are correctly centered.\n\n- The role of Mn in all results is a little unclear. An example Mn and its implications would be useful for the reader.\n\n$M_n$ can be any sequence tending to infinity, arbitrarily slowly. One should think of this as 'a large enough constant' in front of the rate \n$n^{-\\frac{\\min(\\alpha,\\gamma)}{2\\gamma+1}}$\n, but for technical reasons when proving posterior convergence rates, it is often easier to take this `large constant' as something tending to infinity very slowly.\n\n- All other minor points.\n\nWe have corrected these subject to space."}, {"review_id": "BdeV92hyhE", "reviewer": "Reviewer_BMyX", "summary": "\nThe paper gives algorithms for point wise estimation and uncertainty quantification (error bars) for \nsparse variational Gaussian processes with the Brownian motion as prior. The SGVP model though widely used in practice, lacks theoretical guarantees and this paper is the first to establish error guarantees for the model with Brownian priors. Experiments are presented to validate the results. The priors for SVGP models used in practice are not Brownian, extending the results in this paper to these priors remains an open question. \n", "questions": "\nPg 2, the wedge notation for the minimum, the definitions of smoothness and rate should be introduced before stating main results. \n\nIt may be helpful to look at the stochastic process literature where the various properties of BM were originally developed, there may be other stochastic processes like the Ornstein-Uhlenbeck process considered in the stochastic process literature whose spectra are similar to the BM, however these do not seem to have not been used for SGVP. \n\n", "limitations": "Limitations are discussed clearly by the authors. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper is among the first to address rigorously the statistical inference properties of the low rank approximations used in Gaussian processes. \n\nThe main result establishes rigorous bounds on the size of the confidence interval for the SGVP prediction given bounds on the smoothness of the prior and the ground truth and the number of eigenvalues retained in the approximation of the kernel matrix. This is a new and useful result for SGVPs. \n\n- The analysis is particular to the case of Brownian motion priors. BM priors have eigenvectors that can be expressed in closed form (they correspond to the trigonometric functions) and have strong self-similarity properties. However, in the experiments it is observed that the qualitative behavior that is observed for the BM priors is also observed for the Matern and squared exponential kernels used in practice. \n", "weaknesses": "\n- The result does not extend in a straightforward way to the Matern and squared exponential kernels used in practice with SGVPs. Brownian priors have many special properties that are not shared with other Gaussian processes used in ML. \n\n- The result although a useful addition to the literature on SGVPs do not seem not be of enough general interest for a clear accept \nand need further development for analyzing the instances of SGVPs used in practice. \n\n- The results are non adaptive as they require a pre-specified prior smoothness which is not a realistic assumption in practice. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive suggestions and helpful comments. Our response:\n\n- The result does not extend in a straightforward way to the Matern and squared exponential kernels used in practice with SGVPs. Brownian priors have many special properties that are not shared with other Gaussian processes used in ML.\n\nWe fully agree it would be better to derive results for the Matern or squared exponential, but we do not currently have the tools to do so. We compromised by studying a more tractable prior sharing many similarities with the Matern. This allowed us to derive precise theoretical results, whose main messages seem to also apply to the Matern, as we verify numerically. Please see our general reponse for more details.\n\n- The results are non adaptive as they require a pre-specified prior smoothness which is not a realistic assumption in practice.\n\nWe fully agree and plan to extend our results to adaptation in future. However, it is typical to first consider the non-adaptive case both for simplicity and also to understand the effect of the prior smoothness choice on the results. Indeed all the SVGPR theory papers we have cited (Burt et al ICML 2019, Nieman et al. JMLR 2021, Nieman et al. arXiv 2022, Vakili et al. ICML 2022) are non-adaptive.\n\nThe natural approach would be to select the parameter $\\gamma$ by maximizing the marginal likelihood or evidence lower bound (e.g. how the lengthscale is usually picked for the Matern). We have run some simulations in this way, which seem to match the messages we derive in this paper. So our results should extend to adaptive procedures for 'nice' truths.\n\n However, UQ theory for adaptive procedures is a subtle topic, which typically requires further assumptions (e.g. self-similarity) on the true generative function to even be possible. There exist functions that can trick an adaptive procedure into thinking they are smoother than they really are, and consequently make the resulting UQ procedure overconfident (= bad/misleading), see e.g. the discussion paper of Szabo, van der Vaart and van Zanten (AOS 2015). Given such functions will also cause the adaptive SGPR to fail here (like any adaptive procedure), we prefer not to include partial simulations in this manuscript that may not give the full picture of what happens.\n\nExtending our results to the adaptive setting will require significant technical and conceptual work, which is beyond the scope of our manuscript. In particular, one might hope that taking fewer inducing variables might counteract the overconfidence in such settings when faced with bad truths.\n\n- Pg 2, the wedge notation for the minimum, the definitions of smoothness and rate should be introduced before stating\nmain results.\n\nWe changed $\\wedge$ to $\\min$ in the introduction. We prefer not to introduce the rather lengthy definition of Holder smoothness and rate here since we feel this would break the reader's flow in the introduction.\n\n- It may be helpful to look at the stochastic process literature where the various properties of BM were originally developed, there may be other stochastic processes like the Ornstein-Uhlenbeck process considered in the stochastic process literature whose spectra are similar to the BM, however these do not seem to have not been used for SGVP.\n\nThank you for this is an interesting point. We will have a look to see how far these results can be extended, but have not had a chance to do so yet. "}, {"review_id": "i2vVQmLYds", "reviewer": "Reviewer_sbVc", "summary": "The paper gives some theoretical results for frequentist coverage of the credible sets for sparce variational Gaussian processes using a class of Brownian motion based priors. It provides some conclusions relating the smoothness of the prior vs the smoothness of the true underlying function and the resulting ability of the Bayesian inference to produse reliable frequentist uncertainty estimates. ", "questions": "It would be useful for my own understanding if the authors commented more on the statement that $L_2$-type credible sets are less reflective of actual practice than pointwise credible sets. \n\nI would also be curious why the authors chose this venue for publishing this work, given the theoretical nature of the results. ", "limitations": "It is possible that the paper will not be very accessible to the broader audience that it is addressing due to the theoretical nature of the main results. That said, I think the authors did a good job highlighting the exact contribution in the introduction and it should be accessible for practitioners in informing them of the possible consequences in the estimates of credible intervals based on the choice of priors. As the choice of prior is typically difficult in practice, these types of results are valuable in practice.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper is very well written and easy to follow at a high level. The theoretical results are summarised well and give a good indication of the overall contribution of the paper. \n\nThe contribution seems to be an important starting point for providing guarantees for pointwise uncertainty estimates in SVGPs. It gives enough details to indicate the exact extent to which the results are valid. \n\nThe theoretical results are provided with a sufficient level of detail to be accessible by a more general (less statistically inclined) audience.\n\nThe discussion provides a number of directions for future work that can be taken up by the community.", "weaknesses": "The main motivation for the work is the wide use of SVGPs in practice. However, the results in the paper consider only very rough (BM) priors which are rarely used in practice, and only touch on smoother, more commonly used options such as Matern kernels. Admittedly, it is still an important contribution that further work can build on. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive suggestions and helpful comments. Our response:\n\n- The main motivation for the work is the wide use of SVGPs in practice. However, the results in the paper consider only very rough (BM) priors which are rarely used in practice, and only touch on smoother, more commonly used options such as Matern kernels. Admittedly, it is still an important contribution that further work can build on.\n\nWe fully agree it would be better to derive results for the Matern or squared exponential, but we do not currently have the tools to do so. We compromised by studying a more tractable prior sharing many similarities with the Matern. This allowed us to derive precise theoretical results, whose main messages seem to also apply to the Matern, as we verify numerically. Please see our general reponse for more details.\n\n- It would be useful for my own understanding if the authors commented more on the statement that $L_2$-type credible sets are less reflective of actual practice than pointwise credible sets. \n\nMany theoretical results consider $L_2$-credible sets for the entire function $f$ of the form $C_n = \\\\{ f: ||f-\\hat{f}||_{L_2} \\leq q_n \\\\}$ with $\\Pi(C_n|D_n) = 0.95$. The structure of theses $L_2$-sets interacts nicely with Gaussian processes, since a GP can be naturally expanded in terms of an $L_2$ basis via its Karhunen-Loeve expansion. This allows an easier mathematical analysis. However, such credible sets are not so intuitive in practice, e.g. they contain unbounded functions and one cannot actually plot an $L_2$-ball. Indeed, two functions can be extremely close in $L_2$, but arbitrarily far apart at any given point $x_0$, which does not match the more practical intuition of functions being 'close'. The more usual approach in practice is to plot credible bands about the posterior mean, which is a pointwise concept (see the figure in the pdf file attached to the general author response). Pointwise UQ also has the interpretation of measuring the uncertainty of a prediction $f(x_0)$ for a given input $x_0$, which the global $L_2$ set does not.\n\n- I would also be curious why the authors chose this venue for publishing this work, given the theoretical nature of the results.\n\nSparse variational GPs are widely used in the ML community, but much less so in the traditional statistics community. Hence we chose to submit to an ML venue to try and reach the community actually working with such methods. Given the length of our work, we felt an ML conference would be an appropriate venue for our submission.\n\n- It is possible that the paper will not be very accessible to the broader audience that it is addressing due to the theoretical nature of the main results. That said, I think the authors did a good job highlighting the exact contribution in the introduction and it should be accessible for practitioners in informing them of the possible consequences in the estimates of credible intervals based on the choice of priors. As the choice of prior is typically difficult in practice, these types of results are valuable in practice.\n\nThank you for positive comments. We have indeed tried to clearly summarize our main messages, which is always a challenge for theoretical results."}, {"review_id": "eXXLESiG1h", "reviewer": "Reviewer_oD8s", "summary": "Broadly, the paper addresses the issue of the quality of the approximate posterior of a sparse GP model in terms of uncertainty qualification (UQ). The paper chooses the setup with eigenvector-inducing variables and rescaled Brownian motion prior. However, as shown, the results also extend to other kernels (squared exponential (RBF) and Matérn kernels). Precisely, the paper's main conclusion is that with a well-calibrated prior and a sufficiently large number of inducing variables, the uncertainty quantification (UQ) obtained is reliable, though it can be conservative. The paper presents theoretical results on the rate and quality of SVGP posterior, nature of credible sets, and contraction rate of approximate posterior. The theoretical results are shown to match the empirical results based on experiments on simulation data and semi-simulation data.", "questions": "* Is the metric for the quality of the approximate posterior different from the metric of UQ? Is having a “good” quality approximate posterior but not a well-calibrated UQ possible? As the authors also point out in the paper, there has been a lot of recent research evaluating the quality of approximate posterior and convergence rates. Is it much different?\n* What are the connections with Burt et al. (2019), and where do the authors agree and disagree with them? In general, I believe there should be more connections drawn with them.\n* A minor comment, references to the appendix need to be included. I do not see any reference to the appendix in the main paper, so connecting derivations and their corresponding sections becomes tricky.\n* More practical question: how to select the value of $\\alpha, \\gamma$ for a real-world dataset, let us say a UCI data set. I believe more evaluations should be done by extending section 4. It would make the work more applicable, in my opinion.\n* More questions are written in the weakness section.\n", "limitations": "The authors do not discuss them, and I do not see a limitation directly of this work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The topic of the paper is relevant as Gaussian processes are go-to models for many tasks, including Bayesian optimization, active learning, and various sequential tasks. Uncertainty quantification is one of the crucial properties of the Gaussian process models that make them popular for these tasks. Recently, a lot of research has been done on the quality of the approximate posterior of sparse GP models, and this paper rightly fits in there in terms of evaluating the quality of approximate posterior in terms of uncertainty quantification. The claims of the paper are sensible, and the detailed derivations in the appendix also look logical (I would like to see what other reviewers think about the proofs as that is not my core area). Apart from a few places, the paper's notations and flow are good and convenient to follow. The experiments are also chosen wisely and demonstrate the theoretical results.", "weaknesses": "A couple of weaknesses that I find in this paper are:\n\n* The paper talks about sparse variational Gaussian process models (SVGP) throughout, but actually, the model considered in the paper is sparse Gaussian process regression (SGPR). It can be clearly seen from Eq (3). It is better to make this clear in the paper. With a non-Gaussian likelihood, the equations for the posterior are different from those presented here.\n* Section 2.1 has issues with the notations. \n  - As pointed out earlier, the equations are of SGPR and not SVGP. \n  - L104 mean of the GP `m` is a typo?\n  - L117-L118 starts writing cov. Should it be `k`?\n  - $K$ is the gram matrix and should be bold as I believe the paper tries to follow that notation by making $k_n(x)$ bold.\n* A key takeaway message of the paper is missing. There are a bunch of theorems, derivations, and contributions, but I am still trying to figure out the main takeaway message from the paper that I should remember next time I plan to use a sparse GP model.\n* As we discuss uncertainty quantification, does it make sense to report negative log predictive density (NLPD) in the experiments along with the RMSE?\n* A plot of the true posterior, approximate posterior, and credible set would be good to demonstrate the setup and the output. \n* Theorem 3.5 derives a bound for the number of inducing variables `m`. How does it relate to Burt et al. (2019), as they also derive a bound for the number of inducing variables?\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive suggestions and helpful comments. Our response:\n\n- The paper talks about sparse variational Gaussian process models (SVGP) throughout, but actually, the model considered in the paper is sparse Gaussian process regression (SGPR).\n\nWe have changed our terminology to SGPR as per the reviewer's suggestion.\n\n- Notational issues.\n\nWe have fixed all of these, thank you.\n\n- A key takeaway message of the paper is missing. There are a bunch of theorems, derivations, and contributions, but I am still trying to figure out the main takeaway message from the paper that I should remember next time I plan to use a sparse GP model.\n\nWe have tried to clearly summarize these via bullet points in the introduction. In summary, we give some first ideas of the various regimes one can expect for UQ using SGPR (reliable but conservative UQ in some cases, bad UQ in others). The practical messages also seem to apply in practice to the widely used Matern kernel. Please see our general response for more details.\n\n- As we discuss uncertainty quantification, does it make sense to report negative log predictive density (NLPD) in the experiments along with the RMSE?\n\nWe have added this to the simulations, thank you.\n\n- A plot of the true posterior, approximate posterior, and credible set would be good to demonstrate the setup and the output. \n\nThis is a very good idea. We would have liked to add this to the main text, but were unable due to space constraints. We have hence added this image in Section 8 of the supplement (see the figure attached with our general author response).\n\n- Theorem 3.5 derives a bound for the number of inducing variables m. How does it relate to Burt et al. (2019), as they also derive a bound for the number of inducing variables?\n\nPlease see our general response.\n\n- Is the metric for the quality of the approximate posterior different from the metric of UQ? Is having a “good” quality approximate posterior but not a well-calibrated UQ possible? As the authors also point out in the paper, there has been a lot of recent research evaluating the quality of approximate posterior and convergence rates. Is it much different?\n\nYes, these are fundamentally different questions, please see our general reponse (\"Connection with Burt et al. 2019\"). If the true posterior is bad for UQ, a good quality approximation will also be bad. One can also find examples (outside the GP setting) where a variational posterior performs better than the true posterior, e.g. Section 5.2 of Zhang & Gao (Annals of Statistics 2020).\n\n- What are the connections with Burt et al. (2019), and where do the authors agree and disagree with them? In general, I believe there should be more connections drawn with them.\n\nPlease see our general response. We have added some brief comments in the text, but these were sadly limited due to space constraints.\n\n- A minor comment, references to the appendix need to be included. I do not see any reference to the appendix in the main paper, so connecting derivations and their corresponding sections becomes tricky.\n\nWe placed all proofs in Section 7, as indicated in the 'Organization' paragraph on p3. Since all proofs are in the supplement and the section is not too long, we prefer not to add specific references after each result.\n\n- More practical question: how to select the value of $\\alpha,\\gamma$ for a real-world dataset, let us say a UCI data set. I believe more evaluations should be done by extending section 4. It would make the work more applicable, in my opinion.\n\nThe natural approach would be to select $\\gamma$ by maximizing the marginal likelihood or evidence lower bound (e.g. how the lengthscale is usually picked for the Matern). We have run some simulations in this way, which seem to match the messages we derive in this paper. So our results should extend to adaptive procedures for 'nice' truths.\n\nHowever, UQ theory for adaptive procedures is a subtle topic, which typically requires further assumptions (e.g. self-similarity) on the true generative function to even be possible. There exist functions that can trick an adaptive procedure into thinking they are smoother than they really are, and consequently make the resulting UQ procedure overconfident (= bad/misleading), see e.g. the discussion paper of Szabo, van der Vaart and van Zanten (AOS 2015). Given such functions will also cause the adaptive SGPR to fail here (like any adaptive procedure), we prefer not to include partial simulations in this manuscript that may not give the full picture of what happens.\n\nRegarding theory, we plan to extend our results to adaptation (i.e. picking the smoothness parameter in a data-driven way) in future. However, it is typical to first consider the non-adaptive case both for simplicity and also to understand the effect of the prior smoothness choice on the results. Indeed all the SVGPR theory papers we have cited (Burt et al ICML 2019, Nieman et al. JMLR 2021, Nieman et al. arXiv 2022, Vakili et al. ICML 2022) are non-adaptive.\n\nExtending our results to the adaptive setting will require significant technical and conceptual work, which is beyond the scope of our manuscript. In particular, one might hope that taking fewer inducing variables might counteract the overconfidence in such settings when faced with bad truths.\n\nRegarding $\\alpha$, this is the smoothness of the unknown true function but is not needed for the method ($\\gamma$ plays this role and can be picked as described above). While our bound for the number of inducing variables required does depend on $\\alpha$, one can use a weaker but still meaningful bound. Since $\\alpha \\mapsto \\frac{2+\\alpha}{1+\\alpha}$ is decreasing on $(0,\\infty)$, our bound satisfies $n^{\\frac{1}{1+2\\gamma}\\frac{2+\\alpha}{1+\\alpha}} \\ll n^\\frac{2}{1+2\\gamma}$, so that it is enough to take at least this many inducing variables, which only depends on $\\gamma$. Note that for $\\gamma>1/2$, this will be less than linear in $n$."}]

Author rebuttal:
We thank all reviewers for their constructive suggestions and helpful comments. We answer here the main issues raised by reviewers. The main changes we made are to better motivative the use of the rescaled Brownian motion prior and further discuss the connections with the Matern.

- Motivation behind using the rescaled Brownian motion (rBM) prior.

Our main motivation for using this prior is that it provides a compromise of tractable formulas allowing us to obtain precise theoretical results, while simultaneously sharing many properties with the widely used Matern GP, meaning our main qualitative conclusions also apply there. We do not necessarily advocate using rBM in practice, since the prior draws will be very rough, but our results do provide a template for what to expect when using the Matern GP in practice, which we verify is true numerically in simulations.

The tractability of rBM allows us to derive precise results, e.g. that pointwise credible sets from properly calibrated SVGPR's will be conservative (larger than optimal), and thus require no inflation as has been proposed elsewhere (e.g. Vakili et al. ICML 2022). Simulations indicate similar results holds for the Matern as predicted, but we do not know how to prove this rigorously. Since this work is, to the best of our knowledge, the first establishing exact theoretical coverage guarantees for pointwise credible sets for any SVGPR and seems to carry the right practical message for the Matern, we hope this is already an interesting message.

We expanded this motivation in our introduction and Section 4, further discussing connections with the Matern, e.g. the rescaling we consider plays the analogous role to the lengthscale parameter for the Matern. 

- Extending these results to other Gaussian process priors beyond rescaled Brownian motion.

Our analysis is particular to rescaled Brownian motion since our proofs rely on explicit expressions for the eigenvectors/values for this GP. Since these are not typically available for other GPs our method of proof will not extend to these, and new ideas are required. We note that even for the true posterior, few pointwise UQ results are available, e.g. we are not aware of any exact rigorous coverage guarantees for the Matern process in nonparametric regression, let alone for sparse approximations. We first need to develop tools to deal with the true posterior before tackling the harder case of SVGPR approximations. While our proofs are indeed specific, these are the first exact pointwise inference results for SVGPRs and seem to give the correct message for the Matern prior, as we verify numerically.

- Clarify the main contributions of the paper. What results are known and which are new?

We tried to summarize our main conclusions in bullet points in the introduction. From a technical perspective, our main contributions are to derive exact asymptotic UQ guarantees and convergence rates for pointwise inference using SVGPR based on the rBM prior for different regimes depending on the smoothness of the prior and truth. As far as we know, none of our main results (theorems, propositions, etc.) are known elsewhere.

As pointed out by a reviewer, we expect Lemma 3.1 (general formulas for bias and variance of SVGPR) may be known to experts, but we haven't seen it written down exactly in this form. Given it is central to our proof, we have kept it.

- Connections with Burt et al. (ICML 2019).

Burt et al. consider approximation quality as measured by the KL-divergence between the SVGPR and true posterior. We consider the different question of whether the SVGPR behaves well for pointwise inference, even when it is not a close KL-approximation. In particular, the approximation quality may be good but the SVGPR performance will be bad if the true posterior also performs badly. Similarly, the approximation quality may not be good in KL-sense, but the SVGPR posterior can still perform well for pointwise inference.

Our results provide an example of the last situation (see Remark 3.7), where we take sufficiently few inducing points such that the KL-divergence with the true posterior increases with the sample size $n$, but the SVGPR still provides as good pointwise inference. Of course, if the true posterior performs well and the approximation is good (as in Burt et al.), then the SVGPR will also perform well. But our results show this is not necessary.

Their data-generating setup is also not directly comparable to ours. They take expectations over both the data $(x,y)$ and prior on $f$, so their results are 'averaged' over the prior. This is a more Bayesian way of thinking and such results do not translate into frequentist guarantees about convergence rates or UQ, which hold assuming the data comes from a fixed function $f_0$ (see our Assumption 2.1). This Bayesian approach allows one to ignore the method's performance for certain functions as long as they occur with small prior probability, whereas the frequentist (minimax) perspective seeks to understand how it behaves for a given (often worst-case) function. Such setups can lead to different conclusions, see e.g. Chapters 6.2-6.3 in the Bayesian nonparametric textbook of Ghosal and van der Vaart (2016). 

In conclusion, we answer a fundamentally different question and work in a different setting from Burt et al, so our results are not directly comparable with theirs. We have only added some brief discussion in our manuscript due to space constraints.

Regarding frequentist results, Nieman et al. (JMLR 2021, arXiv 2022) study convergence rates and UQ for SVGPR with $L_2$ credible sets. These results crucially exploit the $L_2$ structure, which allows easier mathematical analysis, but is less reflective of practice. For pointwise inference, Vakili et al. (ICML 2022) obtain nice guarantees in a different bandit-like setting by enlarging credible intervals. In contrast, our results suggest the uninflated sets are already larger than needed.

Meta-review: This paper addresses a challenging but important topic: obtaining uncertainty quantification guarantees for commonly used heuristic GP approximations. Although the rescaled Brownian motion prior is non-standard, the theoretical results nevertheless provide valuable qualitative insights.
Decision: Accept (poster)