Sparse Gaussian Process Hyperparameters: Optimize or Integrate?
Authors: Vidhi Lalchand
Wessel Bruinsma
David R. Burt
Carl Edward Rasmussen
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian processes
Sparse Gaussian processes
regression
Hamiltonian Monte Carlo
probabilistic inference

Abstract: The kernel function and its hyperparameters are the central model selection choice in a Gaussian process (Rasmussen and Williams, 2006).
Typically, the hyperparameters of the kernel are chosen by maximising the marginal likelihood, an approach known as Type-II maximum likelihood (ML-II). However, ML-II does not account for hyperparameter uncertainty, and it is well-known that this can lead to severely biased estimates and an underestimation of predictive uncertainty. While there are several works which employ fully Bayesian characterisation of GPs, relatively few propose such approaches for the sparse GPs paradigm. In this work we propose an algorithm for sparse Gaussian process regression which leverages MCMC to sample from the hyperparameter posterior within the variational inducing point framework of (Titsias, 2009). This work is closely related to (Hensman et al, 2015b) but side-steps the need to sample the inducing points, thereby significantly improving sampling efficiency in the Gaussian likelihood case. We compare this scheme against natural baselines in literature along with stochastic variational GPs (SVGPs) along with an extensive computational analysis.  




Sparse Gaussian Process Hyperparameters: Optimize or Integrate?
Vidhi Lalchand
Wessel P. Bruinsma Microsoft Research Al4Science
Department of Physics University of Cambridge
Coruissma@mccrooofft.. com
vr3080cam.adt uc
David R. Burt LIDS Massachusetts Institute of Technology dburt edu
Carl E. Rasmussen Department of Engineering University of Cambridge cer540cam ac uk
Abstract
The kernel function and its hyperparameters are the central model selection choice in Gaussian process [Rasmussen and Williams, 2006]. Typically, the hyper- parameters of the kernel are chosen by maximising the marginal likelihood, an approach known as Type-Il maximum likelihood (ML-II). However, ML-II does not account for hyperparameter uncertainty, and well-known that this can lead to severely biased estimates and underestimation of predictive uncertainty. While there are several works which employ a fully Bayesian characterisation of GPs, relatively few propose such approaches for the sparse GPs paradigm. In this work we propose an algorithm for sparse Gaussian process regression which leverages MCMC sample from the hyperparameter posterior within the varia- tional inducing point framework of [Titsias, 2009]. This work is closely related Hensman al. [2015b], but side-steps the need sample the inducing points, thereby significantly improving sampling efficiency in the Gaussian likelihood case. We compare this scheme against natural baselines in literature along with stochastic variational GPs (SVGPs) along with an extensive computational analysis.
Introduction
Gaussian processes (GPs) are prominent class of models for supervised learning which can quantify uncertainty and incorporate inductive biases function space via the kernel function. Hand-crafting kernel function is powerful way incorporate prior knowledge. In many instances not all properties of kernel function can be specified from prior knowledge alone, and parameters are chosen via ML- II. However, defining complex kernel function with large number of hyperparameters the marginal likelihood prone to multiple local optima and overfitting. Further, several local optima may correspond to priors that do not sensibly model the data. Weakly identified hyperparameters can manifest in flat ridges in the marginal likelihood surface making gradient based optimisation extremely sensitive to starting values [Warnes and Ripley, 1987]. Overall, the ML-II point estimates for the hyperparameters are subject to high variability and underestimate prediction uncertainty.
make
The problem of ridges in the marginal likelihood surface also does not necessarily go away more observations are collected. For example, fi and f2 are Brownian motions, a fi (x/E) is equal in distribution to aaa f2(x//f), which means observations do not provide any information about the product oé. More generally, for a greater class of kernels, including the Matémn1/2 kernel, ff (x/E)
where different combinations of hyperparameters give very similar marginal likelihood values
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Train marginal likhhooo Lengthscalel
size
szee
Lenghoods
Lengheaae
Lenghhcae
30
10
10a
19a
109
10a
10a
102
10F
10
102
5kg
5kg
5ig
Figure I: Negative log marginal likelihood surface as a function of two hyperparameters: of and for squared exponential kernel and Id function. The red cross indicates the truc hyperparameters. The hyperparameters selected via gradien- based optimisation are sensitive the initialisation due to the long ridge of almost identical height at values the hyperparameters not concordant with the ground truth.
equivalent to yāo (aa//()) which means that not possible to consistently estimate the product ol from data, no matter how many observations are collected in a fixed domain [Chapter 6, Stein, 1999]. This implies that one cannot estimate the individual hyperparameters (o, E) consistently. It also motivates why there can be benefits to estimating the hyperparameter posterior even in large data regimes, and ML-II may be insufficient. A more satisfactory treatment of hyperparameters involves placing a prior over the hyperparameters and performing Bayesian inference to compute thyper)posterior For large datasets, this motivates using scalable GP inference (e.g. sparse methods) in conjunction with Markov chain Monte Carlo (MCMC) for the hyperparametrs.
The Bayesian treatment of weakly identified hyperparameters may also be fraught with difficulties. Gradient-based samplers like Hamiltonian Monte Carlo [Neal et al., 2011] and its variants have difficulty navigating regions of high curvature and flat ridges where the gradient offers no information for transition [Betancourt, 2017]. This leads to over-concentration of samples from the flat region. (Usually, this can at least partially be rectified with informative priors.) These pathologies are also typical of other hicrarchical models [Betancourt and Girolami, 2015]. Figure/ 1 shows that the GP marginal likelihood surface can manifest these pathologies. The evidence lower bound (ELBO) used in hyperparameter selection for variational sparse GPs relying inducing points inherits similar, or even less favourable, characteristics to the exact marginal likelihood. As a result, for weakly informative priors, gradient based samplers are susceptible to getting stuck the boundary of these pathological regions hence biasing the sample estimates. The effective sample size metric, used for diagnosing mixing in MCMC, indicative of this behaviour when directly observing the phase space of the target distribution, but is infeasible in high-dimensions.
Historical justification for ML-II (also called the evidence framework) comes from [MacKay, 1994] which highlighted several conditions for ML-II to yield reasonable estimates for the hyperparameters. Crucially, the evidence is unlikely to manifest multiple local optima for model well-matched the data and with a high signal-to-noise ratio. Transferring this insight to the Gaussian process paradigm show how the evidence can have a significant tail or no well-defined maximum in settings with a low signal-to-noise ratio which arises with high aleatoric uncertainty or sparse data, frequently both. In these settings, a point estimate may not adequately summarise the hyperparameter posterior and the benefits of marginalisation stand out. The situation is only exacerbated in high-dimensions (i.e. when there are many hyperparameters) where increasingly more volume of the posterior is captured in thin shell making the density peak extremely unrepresentative of the posterior. Unlike the likelihood of parameters, the marginal likelihood inherently contains trade-off between the fit and complexity penalty term. This is one of the main properties that makes the marginal likelihood objective a viable choice for model selection [Rasmussen and Williams, 2006]. For example, for the Gaussian process regression model,
Yn = (xx) + Ens 6n ~N(0, 02), f GP ko)
the marginal likelihood takes the form,
data fit term
complexity penalty
log 0) log plulf)((110dd -- 1/y (Kg 12 k + 021
(2)
2

---

where c is constant, p(y]f) denotes the data likelihood and 0 denotes kernel hyperparameters. This trade-off is a well-established idea which embodies the automatic Occam's razor effect [Rasmussen and Ghahramani, 2001] where models well-suited to the data are automatically selected just by using the marginal likelihood objective.
This may seem to contradict earlier claims regarding overfitting, but by shying away from dealing with the hyperparameter posterior we risk overfitting even with the marginal likelihood approach. In other words, the evidence framework subdues the overfitting effect induced by the canonical maximum likelihood approach which does not have any complexity penalty term. The parameters in the canonical approach are free to fit the data as well as possible making them prone to overfitting and poor generalisation. Overparameterized kernels based on neural networks like deep kernel models [Wilson et al., 2016] are well known to exacerbate overfitting [Ober et al., 2021].
The main motivation for this work is to highlight that fully Bayesian schemes in sparse Gaussian process models are practically beneficial. While several works in the literature employ fully Bayesian scheme of integrating out the hyperparameters (see table I), in this work we attempt to analyse them in an orthogonal direction, focusing on comparison with the evidence framework and other benchmarks by extending the main sparse variational formulation in the literature Titsias [2009]. We present generalised inference scheme for fully Bayesian GP regression and counteract some of the computational cost of sampling both inducing variables and hyperparameters by deriving a doubly collapsed bound which selects the optimal distribution over the inducing points analytically and targets the kernel hyperparameters with HMC.
2
Related Work
Table 1: Existing literature on fully Bayesian inference in GPs, sparse GPs and generic likelihoods.
Index
Reference
Sparse
Posterior (f/u)
Posterior (0)
Methods
1.
Murra and Adams [2010s]
x
sampling /
sampling
Slice Sampling
2.
Filippone et al. [2013]
X
sampling / NA
sampling
MH + HMC + MA-LA
Filippone Girolami [2014]
Gaussian
sampling
Deterministic: Deeuio-Marginal
4.
Hensman et al. [2015b]
Goussian / sampling
sampling
HMC
Bui al. [2018]
Gaussion / (sampling & VFE)
sampling VFE
MCMC
6.
Lalchand and Rasmussen [2020]
X
Gaussian NA
sampling / VFF
NUTS /VI
7.
Rossi al. [2021]
Gaussian / sampling
sampling
SG-IIMC
8.
Simpson et al. 12021]
x Gaussian / NA
sampling
NUTS / Nested Sampling
9.
This work
Gaussian / (sampling & VFE)
sampling
HMC/ NUTS
Fully Bayesian Gaussian processes have been used by several authors spawning several variants. In early accounts, Neal [1998]. Williams and Rasmussen [1996] explore the integration over covariance hyperparameters using HMC in the regression setting. Barber and Williams [ [1997] extend this the classification setting using HMC for sampling in the hyperparameter space and Laplace approximation for the integrand over function values. Murray and Adams [2010a] and Filippone et al. [2013] focused on MCMC schemes to sample covariance hyperparameters in conjunction with latent function values, mainly mitigating the coupling effect through reparameterisation. Hensman et al. [2015b] considered joint sampling of inducing variables and hyperparameters from the optimal variational posterior distribution while [Bui et al., 2018] consider inference schemes for fully Bayesian sparse GPs in streaming setting. More recently, Rossi et al. [2021] studied fully Bayesian sparse GPs using SG-HMC. Rossi et al. [2021] modify the generative model by adding a prior over the inducing inputs, and peform inference using SG-HMC over the joint (Z, u, 8) space. We list the most recent works in Table 1.
Background
Let f ~ GP (0, ko) be a Gaussian process prior with kernel function Ko depending on hyperpa- rameters 0. We are given noisy observations y (v,)) R of f (f(x))) at input data = ( C R.D. We consider a Gaussian likelihood which factorises over the data, p(ylf) = N (N,,, We wish to compute the posterior p( ly, 0). In this section, we recapitulate the canonical inducing variable approximation of p(Jly, , 0) by Titsias [2009] and its extension to Bayesian treatment of the hyperparameters.
3

---

3.1 Sparse variational inference in Gaussian processes
Following Titsias [2009], we consider a set of inducing variables u [mm)) S R inducing
inputs Z
[2m) m= mm 2m IRd. The complete generative model can then be factored as,
1M
ply f. () p(ylf, (),(f ()p(u(0)
(3)
We approximate the posterior p(,, uly, 0) with a variational distribution:
p(f uly. 0) ~ a(s, 0) = p(f|u, 0))(u)
(4)
where q(u) is chosen minimise the Kullback Leibler divergence (L( p(fly, 0)). Min- imising this KL divergence corresponds to maximising the evidence lower bound [Matthews et al., 2016], henceforth called the ELBO:
ELBO(g(u), 0) [og f. 0)] (L(g(u(0) 11 p(u()))
Because the ELBO still depends on ( this bound is called uncollapsed. Hensman et al. [2013, 2015a] let q(v) be Gaussian, which is optimal the likelihood Gaussian [Titsias, 2009], approx- imate the expectation using Monte Carlo, and maximise the ELBB using stochastic optimisation. On the other hand, the likelihood is Gaussian, Titsias [2009] computes the optimal form for ((u) directly:
q'' () = argmax,yyaa (ax) ELBO(q(u), x p(u(0) exp
p(ylf, 0)]
(6)
Plugging q" (u) back into equation (5), the resulting bound is called collapsed, because it now only depends on and Z. The collapsed bound, denoted La is the objective that Titsias [2009] proposes:
(og((y))) log Nyy: 0, Knm Km1 h K mnn +021)
Tr( Tr(Knn Kym K=1 20
mmm K mmx) =: Le, (7)
where Km is the prior covariance matrix of f, Kmm is the prior covariance matrix and K crosss-covariace matrix formed by f and u. Using the collapsed bound, approximate ML-II
consists of finding,
0" E argmax Lo.z.
Predictions at new functions values can be made in O(M2) after an initial cost of O(N Under certain assumptions on the data generating process, even when M < N, the approximate posterior closely resembles the posterior, and equation (7) is a provably accurate approximation to equation (2) [Burt al., 2020].
3.2 Bayesian treatment of hyperparameters and sparse methods
The extension of the sparse variational framework to a Bayesian treatment of the hyperparameters has been previously considered by Hensman et al. [2015b]. Extend the generative model with prior p(0) over the hyperparameters 0, and let the variational approximation of the posterior p(f. u, Oly) be glf, u. 0) p(f (),(u, 0). The analogue of equation (5)
EL u. 0)) (. [log p(ylf. 0)] Eo (L(()()) [( (L(g(0) ll(())
(9)
and the optimal form for q(u, 0) can again be determined:
q" (u, 0) x p(u, 0) exp
(10)
The distribution q" (u, 0) does not have closed form, and for general likelihoods, Hensman et al. [2015b] propose approximate the expectation in (10) with quadrature and to sample from q (u, 0) using HMC. While this approach quite general, in the case of Gaussian regression, it vastly increases the dimensionality of the state space over which HMC must be relative to HMC in GPR, since the u are sampled addition to the 0. This increases the cost of the procedure, and impacts the success of the sampler.
alternative approach to approximately inferring hyperparameters is to assume a parametric form for q(u, 0) and maximise equation (9) with respect to the variational parameters. Bui et al. [2018]

---

Table 2:
Comparison of variety of approaches to approximating the posterior over hyperparameters in Gaussian process regression. Compares the quality of the posterior (QUALITY); the time complexity per iteration (TIME/IT.); the memory complexity per iteration (MEM./IT.): the number of parameters and/or variables (PARS/VARS); and whether the approach supports non- non-Gaussian likelihoods (LIK.).
APPROACI
QUALITY TIME/IT.
MEM.IT.
PARS/VARS
LIK.
Maximum posteriori [MacKay, 1994]
X
Inducing points; non collapsdd [Titsias and Lázzar-Grediille Greditla, 2014] Inducing points; collapsed [Bui et al., 2018]
+-t
num
m2
m
+
m
n
x
SAMPLING
Exact with Gaussian lik. (Simpson et al., 2021]
+
n2
Exact with non- non-Gaussian [ik[[Muraa und Adams, 2010b] Inducing points; non-collapsed [Hensman et al., 2015b]
+
12
+
m
my?
1he + 172
Inducing points; collapsed (ours)
m
Vhg
took such an approach, assuming that q(u, 0) (()))(()) with both distributions Gaussian. Similar approaches have been applied to variational inference in state-space modelling, sometimes leveraging the optimal form of g( O)) discussed earlier.
The QUALITY column in Table 2 indicates the ability of the method to faithfully represent the hyperparameter posterior. If VI is run to convergence, a potentially significant amount of error will be incurred by the Gaussian approximation to the non-Gaussian posterior over the hyperparameters (red). On the opposite extreme, no sparsity assumption is made, MCMC over the hyperparameters without sparse approximations is asymptotically consistent (green). The inducing point approximations combined with MCMC lie somewhere in-between these methods (yellow).
3.3 Making predictions
The predictive posterior distribution for unknown test inputs X integrates over the joint posterior,
pffl)) ~ f /fff's u. ()p(fl(( ())((())))
(11)
where we have suppressed the conditioning over inputs Xx for brevity. The inner integral simplifies to f If u, ()f(f 0)df = =(ff /u, 0). We discuss the predictive posterior in such models in section 4.3.
4 Fully Bayesian SGPR with HMC: Doubly collapsed formulation
In the previous section, we observed that major drawback of the approach taken in Hensman et al. [2015b] is the need to sample u, which for nigh-dimensional inputs or cases where many inducing points are needed could introduce thousands of additional variables to sample. In this section, we leverage the optimal form of q( (u((0) derived in Titsias [2009] to alleviate this sampling problem.
4.1 Collapsing the evidence lower bound (again)
We first derive the lower bound for this formulation and provide pseudo-code for the algorithm in Algorithm 1. Following the usual derivation of the ELBO,
logp(y) Z / ((0) logr(y 0) (L( (())
(12)
/ (0) zd0 (1) f ( Mo.zp(O) (0) d0 Lz((()), (13)
where log,lyy (0) 2 La.z with La,z defined in equation (7), and where we assign Mo.z eo..2.
4.1.1 Deriving q (0)
We can interpret (Q)) as negative KL divergence as long as we account for a normalisation constant Cz Mo. dd for the un-nomaalised numerator Mo Hence, we can re-write
L (Q)) as,
C2 (a()) log Cz - KL(g(0) (0))
(14)
5

---

where q'(0) -e.zz(())/OCC By inspecting equation (14), we observe that the optimal variational distribution over 0 is given by qº (0) Crucially, by sampling directly from q" (0) using MCMC we climinate the need sample the variables u. By evaluating L3 q" (0), we find the doubly collapsed ELBO 2 : L2 (g) = log Cz. Although the value of L% is computationally intractable, given samples (0,) from 9º (0), gradients of L% with respect to can be estimated using the stochastic estimate of the canonical ELBO equation (7): using the chain rule,
( + (5g hq C2 (0) oz
dz Lz = 02 (0)
(2
~ J J=1
Lo,
(15)
O7
where 8 (2) is the functional derivative of (q) with respect to q, which is zero at q" (0). 5q because q" optimises Ly (it is critical point, so the derivative is zero). Further, the partial derivative of C2 with respect to concerns just the first term of the LHS of equation (13) as the KL term (L(g(() ((O)) is independent of Z.
4.2 Performing approximate inference
We deploy HMC to (approximately) sample from the optimal variational posterior q" (0) along with optimising the inducing inputs Z in hybrid scheme. We alternate between the two steps allocating longer intervals for optimising Z for every HMC sampling run for the hyperparameters. We note that this hybrid scheme is much more computationally efficient than sampling u and 0 jointly where one has to tackle the coupling between inducing variables and hyperparameters joint space. Further, joint sampling is only feasible for moderate number of inducing variables while this scheme can scale to much larger datasets as the efficiency of sampling in the hyperparameter space is only dependent on the dimensionality of the hyperparameter space rather than the number of inducing variables. The entire inference scheme is summarized Algorithm 1, The warm-start strategy (of optimizing both (Z, 0) jointly for few gradient steps) is used to find good region for the sampler to initialise O.
36PR HMC
XOCX0300X
SGPR HMC
Langfwodde CATS CIICE
06-41
MW
M
4
2
040
D12 0150
Figure 2: Top- 1d regression with Left: SGPR, Middle: SGPR + HMC. Right: Joint HMC [Hensman et al., 2015b), Below Left: Samples from the mixture posterior, Middle: Length-scale distribution under SGPR+HMC and MI-II. Note that the data is generated through parameteric function and hence there is ground truth lengthscale. Right: Noise std. deviation distribution from SGPR+HMC and MI-II.
KKL divergence reaches its minimal value of zero when the two input probability distributions are equal, and seek maximise L' (g( which entails minimizing the KL.

---

Algorithm Fully Bayesian Sparse GPR with HMC
1: Input: ELBO objective La.2 L(O, 2) (equation (7)), gradient based optimiser optim() ()
2:
3: procedure WARM-START
4: for fixed number of steps do
5: Gradient step: Z, 0 t- optim(C(0,
6:
return initial values Zine and Anst
7:
8: procedure TRAIN
9: ## Initialisatinn protocol
10: Initialise Co. warm-start values C(0mits Zanie), lets call this C
Freeze kernel hyperparameters
ELBO objective by setting requires grad-False.
12:
13:
while not converged do for 1.... T do start of training loop
14:
16:
Gradient step: Zopt - optim(L) (## equation (15) shows the validity of taking the derivative of the stochastic ELBO)
17:
if t mod == 0 then (## For every gradient stepa)
18:
20:
Draw samples from the optimal hyperparameter variational distribution log g" (0) L(O, Zapi) + p(g) keeping Zopt fixed.
21:
(0g.Ps) P.)
HMC H(0. p); (where is the Hamil monintn
22:
where p denotes the zero-mean auxilliary momentum variable in phase space with the same dimensionality O.
23:
24:
Compute stochastic ELBO 2 j Ein [(0,, Zopt). where 0, ~ (0)
25:
return Zopt, (0)
4.3
Predictive posterior
is ultimately the posterior predictive (PP) distribution that is of interest rather than point estimates or the hyperparameter posterior. The Bayesian sparse GP predictive posterior entails integrating out the posterior over inducing variables and hyperparameters Once we have performed inference, we can approximate this directly,
p(5"lo) /
p( p(f O),(u, Oly)dud0 /
p(f ())((0))))))))(0))))
(16)
mm 1 p( 0)g 004" 0 v)p(0 ddddd NNF K., + A(S* Kmm) ) )A (0
where q- (ul0, 3) = (m,", ( is the optimal Gaussian variational distribution and A Kam Km under the SGPR scheme (Section 4) is available in closed form with m "Kmm 02 kmm) and S" Kmm (Kmm + oZKKn Knm) Kymn- In either case, the internal integral with respect to the inducing variables is analytic and the outer integral be estimated using the samples collected from HMC to perform Monte Carlo estimation. This yields mixture of Gaussians.
0, ~ (0), 1 E N s I0:), p(f =
(17)
n = A (8,) m (8))
yd =
(0,) A -K(8)) - ((
(18)
where J samples are (approximately) drawn from (0) via HMC. The distribution inside the summation the Gaussian posterior predictive distribution for fixed hyperparameters with identical mixing proportions. The compute cost for the predictive posterior scales the sparse GPR cost linearly the number of samples. The Monte Carlo approximation costs O(JN for M inducing points and / hyperparameter samples.
7

---

Deigy
Wodel
lound)
101
5037
1077
1007
2.5
2.5
3.2
0.5
08
230
058
0.2
8.50
0.1
Figure 3: Negative test log likelihoods with standard error of mean across 10 splits with 80% of the data reserved for training. Our method is SGPR + HMC.
5 Experiments
In the previous section, we discussed a hybrid scheme which leverages MCMC within the variational sparse inducing variable formulation leading to fully Bayesian sparse Gaussian processes. In the experiments we demonstrate the feasibility of this scheme relative to several benchmarks and assess regression performance on a 1-dimensional illustrative example and a range of other datasets. We also compare with exact GPs where the cost of using a gradient based sampler is prohibitive for even moderately large datasets requiring several inversions of the full covariance matrix. We show that using the doubly collapsed scheme proposed in this work is a much more attractive alternative for large datasets as compared to direct HMC in GPR, and sampling is more efficient than in the uncollapsed bound used in Hensman et al. [2015b].
We henceforth refer benchmark methods as follows: SGPR HMC refers Bayesian GPs with doubly collapsed variational inference with NUTS, as described in section 3 (compatible with Gaussian likelihoods); we benchmark this model against sparse GPs (SGPR) [Tittsis,, 2009] and Stochastic Variational GPs (SVGP) both using approximate ML-II [Hensman et al., 2015a]. Additionally, we consider the sparse, joint sampling inference scheme proposed [Hensman et al., 2015b) 5b] which gives natural benchmark. We call this model JointHMC. The FBGP method extends the Bayesian treatment to the inducing locations similar to Rossi et al. [2021] 1]; we use NUTS to sample from the joint posterior over (2, 0). It is not direct comparison to Rossi et al. [2021] as the latter explores free-form sampling of u along with (Z, 0) while we work with the collapsed bound incorporating the optimal Gaussian variational distribution q" (u).
We also present analysis where we fix inducing point locations a random subset of the training data (as opposed to interleaving as per Algorithm 1) and only learn hyperparameters using NUTS. We provide several details about the experimental set-up in the supplementary.
5.1 One dimensional synthetic data
We sample noisy observations from f(x) sin(3z) + 0.3 cos(mm)) (TT) with the constraint (r < -2) and (a > 2). Figure 2 shows the results for SGPR along with the fully Bayesian schemes. We keep data split and noise identical across the three models to facilitate a comparison. While there is significant data to identify the hyper- parameters we notice that the models mainly differ their extrapolation abilities away from the training data. SGPR with ML-II overfits to the training data and recovers Jow lengthscale, low noise solution while the SGPR + HMC scheme recovers a more moderate fit and performs
Table 3: Prediction performance in ID synthetic re- gression across SGPR, SGPR HMC and Joi methods with identical number of inducing points and train/test split.
Method
SGPR
SGPR HMC
JointHMC
RMSE NLPD
0.580
0.537 0.065
0.682
0.214
0.74
significantly better in terms of RMSE and NLPD on unseen data Table We note that the MMi
scheme which samples both (u, 0) overfits in the central missing data region. We use M 25 inducing locations across all methods which are optimised according the protocol of each method and recover a very similar spatial distribution.

---

Table 4: A comparison of Sparse GP approaches for UCI benchmarks. RMSE (+ standard error of mean) evaluated on average of 10 splits with 80% of the data used for training. 8 indicates that the posterior over hyperparameters approximated by a point estimate under the respective scheme.
Dataset
GPR HMC
SGPR
SGPR HMC
SVGP
JointHMC
FBGP
[M]
100
100
100
((0)
free-form
free-form (ours)
free-form
free-form (2,0)
Boston
506
13
3.049 (0.14)
3.291 (0.11)
3.286 (0.09)
3.619 (0.11)
3.28 (0.11)
3.845 (0.103)
Concrete
1030
8
4.864 (0.12)
5.459 (0.09)
5.402 (0.05)
5.617 (0.09)
5.612 (0.09)
6.084 (0.11)
Energy
0.441 (0.01)
0.477 (0.008)
0.469 (0.0019)
0.500 (0.01)
0.755 (0.02)
0.490 (0.011)
WineRed
1599
11
0.640 (0.01)
0.636 (0.008)
0.635 (0.008)
0.641 (0.007)
0.641 (0.007)
0.642 (0.007)
Yacht
308
6
0,353 (0.03)
0,412 (0.03)
0.387 (0.03)
0.606 (0.04)
0.794 (0.07)
0.569 (0.037)
Wall dock soconds for total MCMC camples
Johna
SOPR HMC
RMSE
NLPD
0.115
-0,850
SGPR + HMC
SGPR HMC
0.110
Jand--MC
-0.875
JohnHHMM
1DP
-0.900
www.gggaooooooooooooo pert
0.105
-0.925
0.100
-0,950
0.095
-0.975
1IP
0.080
-1.000
M-100
M--200
M=300
M-:00
M-500
M-100
M-200
Mc:30
Mc:400
M-550
Doston
Concrete
Cnengy
Winelee
Figure 4: Left: Test RMSE and NLPD for range of inducing points for the Elevator dataset. Right: Sampling performance measured in terms of the time it took to draw the combined set of samples during the training phase (excluding tuning) as defined by the python standard library cime.terf_counte (wall time). We use the pymc3 pm. NUTS sampler for GPR HMC and SGPR + HMC, and tfp. mcme. Hamil i toni anMonteCarlo for JointHMC [Matthews et al., 2017]. All experiments were conducted on Intel Core 17-10700 CPU @ 2.90GHz x 16.
5.2 UCI regression benchmarks
We compare our approach across methods on 5 standard small to medium-sized UCI benchmark datasets. Following common practice, we use a 20% randomly selected held out test-set [Rossi et al., 2021, Havasi et al., 2018] and scale the inputs and outputs to zero mean and unit standard deviation within the training set (we restore the output scaling for evaluation) [Salimbeni and Deisenroth, 2017]. While could use any kernel, we choose the RBF-ARD kernel with a lengthscale for each dimension. For consistency we initialise all the inducing locations (Z) identically across the methods, i.c. by using the same random subset of training data split. We note that adapting the inducing locations brings serious gains in prediction performance versus keeping them fixed (Figure 3). Further, the JointHMC scheme underperforms SGPR (with ML-III) and SGPR + HMC. This is not surprising given that the JointHMC bound equation (10) does not incorporate the optimal setting for quu) and was originally motivated by the need for fully Bayesian scheme for generalised likelihoods. The method SGPR + HMC significantly improves upon JointHMC in the specific Gaussian likelihood case.
Sensitivity to M: We benchmark SGPR+HMC and JointHMC on the Elevator dataset (N = 16599, D 18) which demands a larger M. SGPR+HMC outperforms JointHMC for this dataset across different M but the advantage is more pronounced at smaller M. Further, our method took 1248 vs. 2109 wall clock sec. for the joint scheme for the same number of hyperparameter samples and 500 inducing points.
5.3 Ablation study
In order to understand the efficacy of Algorithm we conduct an ablation study where we perform inference in the same manner, but keeping inducing locations fixed. Algorithmically, this implies that we don't need to compute the stochastic ELBO 2 and just conduct single sampling run for the hyperparameters. The results across 10 splits are summarised in Table 5.

---

Table 5: An ablation study for the doubly collapsed Sparse GPR scheme comparing performance with and without adapting the inducing locations during training. We report test NLPDs and RMSEs over 10 splits.
Dataset
Metric
Boston
Concrete
Energy
WineRed
Yacht
Fixed Z
3.624 (0.110)
6.021 (0.12)
0.499 (0.014)
0.640 (0.007)
0.533 (0.036)
RMSE
Adapt Z (ours)
3.286 (0.090)
5.405 (0.07)
0.469 (0.009)
0.635 (0.008)
0.387 (0.030)
Fixed Z
2.640 (0.020)
3.200 (0.06)
0.696 (0.014)
0.969 (0.012)
0.791 (0. 130)
NLPD
Adapt Z (ours)
2.524 (0.022)
3,065 (0.01)
0.644 (0.013)
0.961 (0.011)
0.391 (0. (0.146)
5.4 Runtimes
While it is possible to train exact GPR with HMC for datasets of this size (in terms of N) it is important to look at the trade-off in terms of compute cost. In Figure 4 we record the average number of wall clock seconds to draw 500 samples under each method. The cost of sampling is (NN) for SGPR + HMC but (NN)) for Exact GPR + HMC. JointHMC deals with higher dimensional phase space (u. 0) hence requires more tuning. We don't include tuning time for fair comparison. For further context we report the total training run-time for our scheme alongside ML-II, GPR + HMC and FBGP in Table 6. The hybrid scheme we propose is significantly cheaper to canonical alternatives with virtually no degradation in predictive performance.
Table 6: Wall clock seconds (this counts all the CPU time, including worker processes in BLAS and OpenMP as defined by the python standard library time for single training split).
Dataset
Boston
Concrete
Energy
WineRed
Yacht
SGPR (ML-II)
22.17 (0.21)
33.06 (0.07)
30.36 (0.09)
39.96 (0.87)
20.41 (0.22)
SGPR HMC (ours)
29.47 (0.34)
53.85 (2.36)
61.60 (1.47)
60.65 (0.63)
24.50 (0.40)
GPR + HMC
78.05 (2.36)
977.40 (13.82)
326.18 (15.87)
1426.25 (39.49)
31.71 (0.59)
FBGP
72.63 (8.29)
156.31 (4.30)
259.81 (11.58)
175.45 (13.14)
101.92 (2.27)
6 Discussion
The evidence framework continues to be the pre-dominant method for training Gaussian processes since their inception into modern machine learning [Rasmussen and Williams, 2006]. While the marginal likelihood is a compelling model selection objective as it offers an inherent trade-off between data-fit and complexity, it is susceptible overfitting and other pathologies leading to biased inference. This work builds on existing methods that combine sparse Gaussian process regression based on inducing variables with Bayesian hyperparameter inference.
Bayesian hyperparameter inference in GPs is however intractable and one has to consider balancing between objectives of computational cost, prediction accuracy and robustness of uncertainty intervals. While in straightforward conditions the fully Bayesian approach might be counter-productive, most real-world applications of GPs rely on engineering sophisticated hand-crafted kernels involving many hyperparameters where there risk of overfitting is pronounced and further, harder to detect. A more robust solution is to incorporate prediction intervals that reflect these uncertainties in the model choice. Studying full Bayesian inference in more sophisticated GP models like deep [Damianou and Lawrence, 2013], warped [Snelson et al., 2004] and convolutional GPs [Van der Wilk et al., 2017] will offer greater insight to this question and an imminent direction of future work.
Acknowledgements
This research was conducted while WPB and DRB were students at the University of Cambridge. During that time, WPB was supported by the Engineering and Physical Research Council (studentship number 10436152). VL acknowledges funding from the Qualcomm Innovation Fellowship (Europe).
10

---

References
David Barber and Christopher K.I. Williams. Gaussian processes for Bayesian classification via hybrid Monte Carlo. In Advances in neural information processing systems, pages 340-346, 1997.
Michael Betancourt. A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint TrX1:::701.02434, 2017.
Michael Betancourt and Mark Girolami. Hamiltonian Monte Carlo for hierarchical models. Current trends in Bayesian methodology with applications, 2015.
Thang D. Bui, Cuong V. Nguyen, Siddharth Swaroop, and Richard Turner. Partitioned variational inference: A unified framework encompassing federated and continual learning. arXiv preprint arXiv: I. 1226, 2018
David R. Burt, Carl Edward Rasmussen, and Mark van der Wilk. Convergence of sparse variational inference in Gaussian processes regression, Journal of Machine Learning Research, 2020.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In Artificial intelligence and statistics, pages 207-215. PMLR, 2013.
Maurizio Filippone and Mark Girolami. Pseudo-marginal Bayesian inference for Gaussian processes. IEEE Transactions Pattern Analysis and Machine Intelligence, 36(1 1):2214-22222 2014.
Maurizio Filippone, Mingjun Zhong, and Mark Girolami. A comparative evaluation of stochastic- based inference methods for Gaussian process models. Machine Learning, 93(1):93-114 14, 2013.
Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matri. Gaussian process inference with GPU acceleration. arXiv preprint wwXiv:1809.1/166 2018.
Marton Havasi, José Miguel Hernández-Lobato, and Juan José Murillo-Fuentes. Inference in deep Gaussian processes using stochastic gradient Hamiltonian Monte Carlo. Advances neural information processing systems, 2018.
James Hensman, Nicolo Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, 2013.
James Hensman. Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. 18th International Conference on Artificial Intelligence and Statistics (AISTATS), 2015a.
James Hensman, Alexander G Matthews, Maurizio Filippone, and Zoubin Ghahramani. Mome for variationally sparse Gaussian processes. In Advances Neural Information Processing Systems, pages 1648-1656, 2015b.
José Miguel Hermández-Lohan and Ryan Adams. Probabilistic backpropagation for scalable learning of Bayesian neural networks. In International conference on machine learning, pages 1861-1869. PMLR, 2015.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint 14X2..1412.6911 2014.
Vidhi Lalchand and Carl Edward Rasmussen. Approximate inference for fully Bayesian Gaussian process regression. In Symposium on Advances in Approximate Bayesian Inference, 2020.
David JC MacKay. Hyperparameters: optimise or integrate out. Maximun entropy and Bayesian methods, 1994.
A. G. G. Matthews, J. Hensman, R. E. Turner, and Z. Ghahramani. On sparse variational methods and the Kullback-Leibler divergence between stochastic processes. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, volume 54 of Proceedings of Machine Learning Research. PMLR, 2016.

---

Alexander G de G Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujui, Alexis Boukouvalas, Pablo León-Vinn Villagrá. Zoubin Ghahramani, and James Hensman. GPflow: Gaussian process library using TensorFlow. Jorunal of Machine Learning Research, 2017.
lain Murray and Ryan P Adams. Slice sampling covariance hyperparameters of latent Gaussian models. In Advances in neural information processing systems, pages 1732-1740, 2010a.
lain Murray and Ryan Prescott Adams. Slice sampling covariance hyperparameters of latent Gaussian models. In Advances in Neural Information Processing Systems 23. Curran Associates, Inc., 2010b.
Radford Neal. Regression and classification using Gaussian process priors. Bayesian statistics, 6: 475, 1998.
Radford M. Neal et al. MCMC using Hamiltonian dynamics. Handbook of Markov chain Monte Carlo, 2(11):2, 2011.
Sebastian W Ober, Carl Edward Rasmussen, and Mark van der Wilk. The promises and pitfalls of deep kernel learning. arXiv preprint arXiv:2102.12108, 2021.
Carl Edward Rasmussen and Zoubin Ghahramani. Occam's 's razor. Advances in neural information processing systems, pages 294-300, 2001.
Carl Edward Rasmussen and Christopher K.I. Williams. Gaussian processes for machine learning (adaptive computation and machine learning), 2006.
Simone Rossi, Markus Heinonen, Edwin Bonilla, Zheyang Shen, and Maurizio Filippone. Sparse Gaussian processes revisited: Bayesian approaches to inducing-variant approximations. In International Conference on Artificial Intelligence and Statistics, pages 1837-1845. PMLR, 2021.
Hugh Salimbeni and Mare Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. Advances in neural information processing systems, 30, 2017.
John Salvatier, Thomas Wiecki, and Christopher Fonnesbeck. Probabilistic programming in python using pymc3. Peer: Computer Science, 2:e55, 2016.
Fergus Simpson, Vidhi Lalchand, and Carl Edward Rasmussen. Marginalised Gaussian processes with nested sampling. Advances in neural information processing systems, 2021.
Edward Snelson, Carl Edward Rasmussen, and Zoubin Ghahramani. Warped Gaussian processes. Advances in neural information processing systems, 16:337-344, 2004.
Michael Stein. Interpolation of Spatial Data. Springer- Verlag New York, edition, 1999.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
Michalis Titsias and Miguel Lázaro-Gredilla. Doubly stochastic variational Bayes for non-conjugate inference. In International conference machine learning, 2014.
Mark Van der Wilk, Carl Edward Rasmussen, and James Hensman. Convolutional Gaussian processes.
arXiv preprint www.:v:100.011100 779.01894, 2017.
Mark Van der Wilk, Vincent Dutordoir, ST John, Artem Artemev, Vincent Adam, and James Hensman. A framework for interdomain and multioutput Gaussian processes. arXiv:2003.01155 2020. URL /etps://arxiv..oo wtw.:orr///bbb/oooeoeeeee 01115.
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Burkner. Rank- normalization, folding, and localization: An improved R for assessing convergence of MCMC (with discussion). Bayesian analysis, 2021.
JJ Warnes and BD Ripley. Problems with likelihood estimation of covariance functions of spatial Gaussian processes. Biometrika, 74(3):440-64,, 1987.
Christopher K.I. Williams and Carl Edward Rasmussen. Gaussian processes for regression. In Advances in neural information processing systems, pages 514-520, 1996.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep kernel learning. In Artificial intelligence and statistics, 2016.
12


Reviews:
[{"review_id": "e4O39mx9lNn", "reviewer": "Reviewer_vWTc", "summary": "The paper proposes a novel framework for tackling Gaussian Process (GP) inference by combining sparse GP (based on inducing variables) and Bayesian hyperparameter inference. The proposed method offers great potential to avoid overfitting compared to previous methods that rely on maximizing the marginal likelihood. Experiments include both synthetic data (one-dimensional) and five small/medium-sized UCI benchmarks. Empirical results show that the proposed methods improve existing sparse GP models.", "questions": "N/A", "limitations": "Yes, the authors have adequately addressed the limitations and potential negative societal impact of their work", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- Considering hyperparameter uncertainty in GP inference is a meaningful and valid direction. The paper thoroughly examined existing methods and proposed a practical method to advance this direction, which brings value to GP's applications.\n- All contributions are well motivated and related works are well cited and documented. The entire paper is self-contained and the presentation is logically ideal.\n- The proposed method is technically correct, with Algorithm 1 helping understand how the model works.\n- Empirical results indicate that the proposed model obtains strong performance on a diverse set of tasks, in particular with reasonable running time.\n\nLimitation: In my opinion, a limitation of the method is that the compared methods are relatively old (SGPR in 2009 [1] and SVGP in 2015 [2]). So it would be interesting to see how this work compares to recent GP inference works (or at least discuss the opportunities of adopting the proposed Bayesian hyperparameter inference to recent works).\n\nReference\n1. Titsias, Michalis. \"Variational learning of inducing variables in sparse Gaussian processes.\" Artificial intelligence and statistics. PMLR, 2009.\n2. Hensman, James, Alexander Matthews, and Zoubin Ghahramani. \"Scalable variational Gaussian process classification.\" Artificial Intelligence and Statistics. PMLR, 2015.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "XJIkksT_HE", "reviewer": "Reviewer_gGgw", "summary": "The paper contributes a specific inference method for sparse variational Gaussian process regression in a Bayesian setup allowing inference of  hyperparameters which doesn't require sampling the inducing variables u. The new algorithm is derived by exploiting the exact bound derived in Titsias, 2009 to \"collapse\" the full ELBO (with hyperarameters) leading to a more effective algorithm in the Bayesain setting.\n\nThe algorithm is evaluated on a relevant set of synthetic and benchmark regression problems which demonstrates that the proposed method outperforms existing sparse inducing points methods in terms of RMSE, predictive log-likelihood, and wall clock time. It is not clear that the approach consistently outperforms a basic SGPR+MLII approach.", "questions": "\nComments/suggestions/questions:\n- l105 mentions the classification case, although only regression with the Normal likelihood is considered going forward. I'd suggest removing or clarifying.\n\n- Table 2: I am not sure I follow what's going on in the table; e.g. how is the \"quality\" determined? Why does MAP not support a non-Gaussian likelihood?... I'd suggest adding a reference per entry where these facts are extracted unless clearly explained e.g. in the appendix. Also, perhaps explicitly indicate which approach is yours.\n\n- Algorithm 1 (page 6): I am slightly confused by the comments and notes in the algorithm box (line 2-3 and 11); are the kernel hyperparameters part of the optimization or not? Additionally, I think it would be beneficial to include reference to the relevant equations (specifically for the \\mathcak{L} elements )\n\n- Figure 2: Please clarify the caption; I (initially) found it difficult to work out exactly what was happening in the various panels.\n\n- Figure 3: There are several SGPR+HMC variants; I assume the authors prefer the reader to focus on \"Adapt 2\" but it would be helpful with some guidance.\n\n- Table 4: $\\delta$ is not explained/defined\n\n- Sec 5.2: I believe the discussion fails to comment on the performance of SGPR (ML-II) vs. the proposed SGPR+HMC (Adapt 2) as observed in table 4 and figure 4? Would this not suggest that the benefit of a fully Bayesian approach is not evident in most cases?\n\n- l250 (figure on the right): It would be informative with the GPR+HMC and SGPR+HMC performance shown as horizontal lines.\n\nVery minor:\n- Conditioning on X seems a bit arbitrary and briefly introduced in 101-108, then disappears in sec 3.1, but back in l130... I'd suggest more consistency/clarity.  \n\n- l101 vs introduction: the text switches from \"kernel\" to \"covariance function\" and it is used interchangeably throughout; perhaps consider if a consistent use would be beneficial\n\n- Eq (1) missing comma before f\n\n- l151: \"...sample the u...\" -> \"...sample u....\"\n\n- l190: check missing punctuation prior to  \"...the outer...\"\n\n- l234: I am not sure where this common practice is used; could you provide a reference?\n\n- l248, l266: \"don't\" -> \"do not\" I would suggest avoiding contractions.\n\n- l261: not sure what \"this\" refers to\n\n\n\n", "limitations": "See points made elsewhere regarding the discussion of results etc.\n\n\n\n------------------------\nFollowing the authors' response, I have updated my score to 6.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\n\n- Novelty and signiificance: The approach appears to be leverage a new insights about the ELBO and contributes a sensible addition to the sparse GP regression toolbox. The paper provides a better understanding of GP inference in the sparse regime which is valuable in its own right, although it is arguably not very far reaching and significant from a practical point of view in part because many people simply resort to MLII which also seems to perform on par with the proposed method on most of the UCI benchmarks.\n- Quality: The basic idea appears sensible and the core algorithmic work well-executed. I find that the experiments could have been presented and discussed better (see specific suggestions/questions below).\n- Clarity and presentation: I find that there are several questions (see below) that I would like to see addressed in order for it to be easily accessible and comphrensive.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "8dMV3WHBYwT", "reviewer": "Reviewer_TWYS", "summary": "The paper addresses the problem of hyperparameter optimization for sparse GP kernels. The proposed solution utilises MCMC to keep the hyperparameter posterior, accounting for uncertainty in the hyperparameter space. The paper introduces a series of decisions, which improve computational efficiency (separating sampling of Zs and \\thetas) and help with handling high-dimensional inputs (by using a doubly collapsed ELBO). The experiments show improvement over alternative techniques.", "questions": "Minor issues:\nEquation 2 is a little bit weird in a sense that \\theta is missing after the first equal sign. I think you should include the dependency on \\theta there for p(y|f) and p(f) too.", "limitations": "Yes.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: The proposed fully Bayesian sparse GP is suitable for larger datasets and has more efficient sampling scheme that the alternatives. The proposed approach is theoretically sound.\n\nThe problem is well motivated and presented. The authors nicely demonstrate the problem of kernel hyperparameter optimization with Figure 1 and further compare relevant works and alternative approaches with Tables 1 and 2.\n\nWeaknesses: It would be interesting to see more comparison methods in the experiments section. The benchmark datasets could be better motivated, e.g. by explaining their main differences, adding their descriptions.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "0Ml-hzYlko", "reviewer": "Reviewer_uTXg", "summary": "This paper proposes a new extension of the variational sparse GP framework (Hensman et al., 2015b) to bypass the need of sampling inducing points. This is achieved by a clever analysis of the optimal variational distribution of hyper-parameter $\\theta$, which then reveals a way to collapse the inducing points $u$ from the ELBO. This doubly collapsed ELBO is intractable, but its gradient can be estimated via sampling $\\theta$ from the optimal variational distribution $q^{\\ast}(\\theta)$ using HMC. This is a significant contribution for inference with high-dimensional data, since it improves the sampling efficiency. The proposed method achieves promising empirical results on several benchmarks.", "questions": "- The warm start strategy was introduced in algorithm 1, but was never described anywhere else. Can the authors elaborate?\n- Why do SGPR + HMC and SVGP use different sampler implementations? Could it be that the runtime/RMSE differences are partly due to this choice?\n- The latest benchmarks included in the empirical study were SVGP (Hensman et al., 2015a) and JointHMC (Hensman et al., 2015b). Since SGPR+HMC extends these methods, I believe it is appropriate to include them as benchmarks. However, since the introduction of SVGP there have been plenty of other sparse GPs. Table 1 even listed a few, such as Rossi et al., 2021 and Bui et al., 2018. I think this paper should additionally compare to these benchmarks to fully showcase the relevance of the proposed method in light of more recent developments.\n- Since I'm not entirely clear about Eq. (14), I'm tentatively putting a score of 5. I will upgrade my score once the authors address my confusion.", "limitations": "This is a purely theoretical work. I do not foresee any potential negative societal impacts.", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is clearly written on a high level, but some details are not quite clear to me. For instance, I do not quite follow what happens at Eq. (14), which was written with extreme brevity. Where does the addition come from? If I understand correctly, it should be an equality there as we substitute $\\mathcal{L}^{\\ast\\ast}$ with $\\mathcal{L}^{\\ast}$ evaluated at $q=q^{\\ast}$. Then, we apply chain rule to get the inner product. Still, I'm quite confused by the final approximation. How do we get to $\\frac{1}{J}\\sum_{j=1}^{J} \\mathcal{L}_{\\Phi,\\theta_j}$ ? The rest of the math seems to check out, so maybe the authors can help to clarify this part. \nMinor issue: In line 120, it should be plugging back to the Eq. above Eq. (5), which is unlabeled.\n\nAs far as I am aware, this work is novel. I really like the clever idea of doubly collapsing the ELBO to avoid sampling inducing inputs. I can see that it would be a problem for inference with high-dimensional data, as we have to draw samples from a very large state space. This was explicitly pointed out by the authors: \"While this approach is quite general, in the case of Gaussian regression, it vastly increases the dimensionality of the state space over which HMC must be run relative to HMC in GPR since the u are sampled in addition to the θ\". However, I find it quite disappointing that the authors only demonstrate on fairly low-dimensional datasets (e.g., the largest was 13-dimensional). I do not see how these benchmarks justify the use of this method.\n\n----------\nI have read the author's rebuttal, which clarified some of my concerns. My overall score will be adjusted to 6. However, I would encourage the authors to continue to revise the empirical study by expanding the range of benchmarks and adding more ablation study to confirm the expected advantage over other methods.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept