Scalable Thompson Sampling using Sparse Gaussian Process Models
Authors: Sattar Vakili
Henry Moss
Artem Artemev
Vincent Dutordoir
Victor Picheny
Venue: NeurIPS 2021 Poster (2021)
Keywords: Thomsson sampling
Gaussian processes
posterior sampling
regret

Abstract: Thompson Sampling (TS) from Gaussian Process (GP) models is a powerful tool for the optimization of black-box functions. Although TS enjoys strong theoretical guarantees and convincing empirical performance, it incurs a large computational overhead that scales polynomially with the optimization budget. Recently, scalable TS methods based on sparse GP models have been proposed to increase the scope of TS, enabling its application to problems that are sufficiently multi-modal, noisy or combinatorial to require more than a few hundred evaluations to be solved. However, the approximation error introduced by sparse GPs invalidates all existing regret bounds. In this work, we perform a theoretical and empirical analysis of scalable TS. We provide theoretical guarantees and show that the drastic reduction in computational complexity of scalable TS can be enjoyed without loss in the regret performance over the standard TS. These conceptual claims are validated for practical implementations of scalable TS on synthetic benchmarks and as part of a real-world high-throughput molecular design task.

Scalable Thompson Sampling using Sparse Gaussian Process Models
Sattar Vakili* Henry Moss Artem Artemev2 Vincent Dutordoir Victor Picheny?
Abstract
Thompson Sampling (TS) from Gaussian Process (GP) models a powerful tool for the optimization of black-box functions. Although TS enjoys strong theoretical guarantees and convincing empirical performance, it incurs a large computational overhead that scales polynomially with the optimization budget. Recently, scalable TS methods based on sparse GP models have been proposed to increase the scope of TS, enabling its application to problems that are sufficiently multi-modal, noisy or combinatorial to require more than a few hundred evaluations to be solved. However, the approximation error introduced by sparse GPs invalidates all existing regret bounds. In this work, we perform a theoretical and empirical analysis scalable TS. We provide theoretical guarantees and show that the drastic reduction in computational complexity of scalable TS can be enjoyed without loss the regret performance over the standard TS. These conceptual claims validated for practical implementations of scalable TS on synthetic benchmarks and as part of a real-world high-throughhut molecular design task.
Introduction
Thompson sampling [TS, 1] is popular algorithm for Bayesian optimization [BO, 2]- sequential model-based approach for the optimization of expensive-to-cvaluae black-box functions. typically characterised by limited prior knowledge and access to only a limited number of (possibly noisy) evaluations. By sequentially evaluating the maxima of random samples from a model of the objective function, TS provides a conceptually simple method for balancing exploration and exploitation.
TS is often paired with Gaussian Processes (GPs), which offers spectrum of powerful and flexible modeling tools that provide probabilistic predictions of the objective function. The resulting GP-TS algorithms [3] have been found to provide highly efficient optimization under heavily restricted optimization budgets, with numerous successful applications including aerodynamic design [4], route planning [5] and web-streaming [6]. While most popular BO algorithms cannot query more than a handful of points at a time [7-10] without employing replicating designs [see 11, 12]. TS has a natural ability query large batches of points. Therefore, TS a popular solution for optimization pipelines enjoying a large degree of parallelisation. for example in high-thie sightthooughout chemical space exploration [13] and for the distributed tuning of machine learning models across cloud compute resources [14].
As BO incurs substantial computational overhead between successive iterations, while updating models and choosing the next set of query points, standard BO methods are limited to optimization problems with small evaluation budgets [2). However, with large batches, the computational overhead incurred by BO per individual function evaluation is considerably reduced. Therefore, considering large batches is promising lactic expand BO larger optimization budgets, which are required to
Eqqull contribution, MediaTek Research, Secondmind, 3 Imperial College London, University of Cambridge. Correspondence to Sattar Vakili <sattar. vakil 1 entkreeaacch com.> Henry Moss <henry www.ss@seconnmmmnddccom ai>>
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

---

optimize highly noisy problems with rougher optimization landscapes [11, 12] or high dimensional and combinatorial search spaces [15, 13, 16]. Consequently, the promising candidate for BO under large optimization budgets.
TS is a
Unfortunately, practical implementations of GP-TS suffer from two key computational bottlenecks that prevent the method from scaling in terms of total optimization budget. Not only does each update of the GP posterior distribution require matrix inversion that incurs a cubic cost w.r.t. the number of obscrvations [17], but even sampling from this postcrior can be a daunting task - the standard approach of drawing a joint sample across N point discretization of the search space has an ( complexity [due to a Cholesky decomposition step, 18]. Alternative existing approaches for BO under large optimization budgets include using Neural Networks in lieu of GPs [15, 13] or to use local models [19] and ensembles [16].
A natural answer to the scalability issues of GP-TS is rely on the recent advances in Sparse Variational GP models (SVGP, 20J. SVGPs provide a low rank O(m²t) approximation of the GP posterior, where m is the number of the so-called inducing variables that grows at rate much slower than f. Successful applications of SVGPs for BO under large optimization budgets include optimizing free-electron laser [21], molecules under yythesis-abiilit constraints [22]. and the composition of alloys [23]. Furthermore, [24] introduced an efficient sampling rule (referred to as decoupled sampling) which can be used to efficiently perform TS with SVGPs. particular, [24] decomposes samples from the SVGP posterior into the sum of an approximate prior based on M features (see
Sec. 3.3) and an SVGP model update, thus reducing the computational cost of drawing a Thompson sample to O ((m + N). Leveraging this sampling rule results in scalable GP-TS algorithm (henceforth S S-GP-TS) that can handle orders of magnitude greater optimization budgets.
While [3] proposed a comprehensive theoretical analysis of exact GP-TS, it does not apply to S-GP- TS. Indeed, using sparse models and decoupled sampling introduce two layers of approximation. that must be handled with care, as even a small constant сιтor in the posterior can lead to poor performance by encouraging under-exploration in the vicinity of the optimum point [25]. Our primary contributions can be summarised as follows. First, we provide a theoretical analysis showing that batch TS from any approximate GP can achieve the same regret order as an exact GP-TS algorithm as long the quality of the posterior approximations satisfies certain conditions (A (Assumptions and 4). Second, for the specific case of S-GP-TS (batch decoupled TS using a SVGP), we leverage the results of [26] to provide bounds in terms of GP's kernel spectrum for the number of prior features and inducing variables required to guarantee low regret. Finally, we investigate empirically the performance of multiple practical implementations of S-GP-TS, considering synthetic benchmarks and a thgoghhhuu molccular design task.
2 Problem Formulation
We consider the sequential optimization of an unknown function f over a compact set X' C Rd. sequential learning policy selects a batch of observation points ( at each time step t 1,2,..... .T and receives the corresponding real-valued and noisy rewards {yr,b = f(x.,)) + E4,bbocn where εg denotes the observation noise. Throughout the paper, we use the notation [r] 21,2,... ,,,, for n E N. As is common in both the bandits and GP literature, our analysis uses the following Suu-Gaussyy assumption, a direct consequence of which that [[ 0, for all t, E N.
Assumption 1. E, are i..dd, over both t and b, R sub-Gaussian random variables, where 0 is fixed constant. Specifically, Eleheca]] S exp( (2R²2 Vh E R, Vt, b E N.
Let x² E ff (z) be optimal point. We can then measure the performance of a sequential optimizer by its strict regret, defined as the cumulative loss compared to fa* over a time horizon T
T "E((*))
R(T, B: E
-
((rap) ((ra))
where the expectation is with respect to the randomness in noise and the possible stochasticity in the sequence of the selected batch observation points b Note that our regret measure (1) is defined for the true unknown f. In contrast, the alternative Bayesian regret [see e.g. 27, 14] averages over a prior distribution for f. As upper bounds on strict regret directly apply to the Bayesian

---

regret (but not necessarily the reverse), our results are stronger than those that can be achieved when analysing just Bayesian regret, for example when applying the technique of [28] that equates TS's Bayesian regret with that of the well-studied upper confidence bound policies.
Following [3, 29, 30J, our analysis assumes regularity condition the objective function motivated by kernelized learning models and their associated reproducing kernel Hilbert spaces [RKHS, 31];
Assumption 2. Given RKHS Hg, the norm of the objective function is bounded: || H3 S B. for some B and k(r, x') ) I L, for all x, x' E X.
In the case of practically relevant kemelss, Assumption 2 implies certain smoothness properties for the objective functions.
3 Gaussian Processes and Sparse Models
GPs are powerful aan-aaametric Bayesian models over the space of functions [17] with distribution specified by a mean function puz) (henceforth assumed to be zero for simplicity) and positive definite kernel (or covariance function) (xx, z'). We provide here a brief description of the classical GP model and two sparse variational formulations.
3.1 Exact Gaussian Process models
Suppose that we have collected a set of location-observation tuples Hr (Xx, Yt), where X; is the tB d matrix of locations with rows [X.](s-1)B+b = &a,b, and y: is the t-B-dimensional column vector of observations with elements [yel(s-1)
Ys.b, for all E [t]. b [B]. Then, assuming a Gaussian observation noise, the posterior of the GP model f given the set of past observations He, is
also a GP with mean (t(+), variance of() and kernel function lie:;)) specified as
(x) = -a, (Kx,,, + -I)="yy. (')x, = k(z, x²) k x (Kx,xx + ++1) +-1)-1ix,, -ix,,
(2)
and (z) Ki (x) with kx... the tB dimensional column vector with entries (
k(t,,,,. r), and Kx,X, the tB x tB positive definite covariance matrix with entries x, (=-1) B))++ "1 We directly see from (2) that accessing the pos-
terior expressions require an O((B))) matrix inversion, which is computational bottleneck for large values of tB.
Note that in problem formulation f is fixed and observation noise has unknown sub- Gaossaassin distribution. Using a GP prior and assuming a Gaussian noise is merely for ease modelling and does not affect our assumptions on f and 6,,6- The notation f is thus used to distinguish the GP model from the fixed
3.2 Sparse Variational Gaussian Process Models with Inducing Points
To overcome the cubic cost exact GPs, SVGPs [20, 32] instead approximate the GP posterior through set of inducing points Z, (z1,..., 2m, (z, E X, with my << tB). Conditioning on the inducing variables
f (Z.) (rather than the tB observations in y:) and specifying prior a S.), yields an approximate posterior distribution that, crucially, still GP but with the significantly reduced computational complexity of O(mît). The posterior mean and covariance of the SVGP are given in closed form as
Gaussian density qe (u,) = - (m,,
M4 (+) (z) lzz,, K Z,,, m. k(( (I, = k(z,x') + ZZ., Z.ZZ (S, Kz,z,, (22.z,)'' ) kiz mz
The variational parameters m, and S, are set as the maximizers of the evidence lower bound (ELBO, see Appendix for the details) and can be optimized numerically with mini-batching [32]. There are various standard ways in practice to select the locations of the inducing points Z. e.g. by using an experimental design, sampling from k-DPP (that stands for determinantal point process), or by optimizing them along with the inducing variables.
3.3 Sparse Variational Gaussian Process Models with Inducing Features
An alternative approximation strategy is using inducing feature approximations [33, 26. 34]. Here, define inducing variables as the linear integral transform of f with respect to some inducing

---

features [35] 1/1 (x),.., Vm, (2), i. we set our 1th inducing variable as W4, fx (b) (r) Courtesy of Mercer's theorem, we can usually decompose our chosen kernel ls as the inner product of possibly infinite dimensional feature maps (see Theorem 1.1 in [36]) to provide the expansion k(x,x') - = Dy= == (x)) (x') for eigenvalues {>, E 1 and eigenfunctions {Ø) E Hx }== we set our inducing features be the my eigenfunctions with largest eigenvalues, can be shown that cov(us,s; vt,5) = Xو6ر,ر and cov(ur,j, f()) 2y0 (a)) yielding an approximate Gaussian Process model with posterior mean and covariance given by
( (x) T (r)m, 14
Ry (x, (z, = -(r,)"))((rx))))) + @ mi
(z) Amm Pm. (x').
Here, m, and S are inducing parameters (as above), Φm. (I) (z)) feature vector and Am, is the m, x m, diagonal matrix of eigenvalues. li.j A.ōt.5.
..., Per (x)] is the truncated
Inducing feature approximations have strong advantages, in particular a reduced computational cost and the fact that inducing points need be specified. However, accessing these eigenfeatures require the Mercer decomposition of the used kernel, which is available for certain kernels on manifolds [37, 34], but limited to low dimensions for others [38, 39].
4 Scalable Thompson Sampling using Gaussian Process Models (S-GP-TS)
At each BO step t, GP-TS proceeds by drawing B i.i.d. samples (frbhe)) from the posterior distribution of f and finding their maximizers, i.e. we select samples Fe,6 satisfying
{az,5 f.
(3)
However, since fe,o is an infinite dimensional object, generating such samples is computationally challenging. Consequently, it is common resort approximate strategies, the most simple of which is to sample across an N₂ point discretization D, of X [14] which can be obtained with an ((N7)) cost (due required Cholesky decomposition).
To improve the computational efficiency of TS, a classical strategy [40, 41] is to rely on kernel decompositions, For instance, sample f from GP can be expressed as a randomly weighted sum of the kernel's 's cigenfunctions fxx) = Σ/=1 ()yu;d)(2)) or, in the case of shift-invariant kemels, the kernel's Fourier features (x) (see [42]) as f(x)= = Eÿ= myyy (x)) By truncating these infinite expansions to contain only the M cigenfunctions with largest cigenvalues or M random Fourier features, we have access to approximate but analytically tractable samples. For both expansions, the weights wy are sampled independently from a standard normal distribution. Conditioned on current tB observations, the posterior distribution of wy are Gaussian with mean and covariance functions that can be calculated with (M computations, resulting in O(M3 + BNM) cost to draw B Thompson samples.
Fast approximation strategies described above avoid costly matrix operations and work best only when sampling from GP priors. Posterior GP distributions are often too complex be apll-approximated by finite feature representation [16, 43, 30]. The recent work of [24] tackled this issue by using truncated feature representations only to approximate the prior GP and a separate model update term to approximate posterior samples. For SVGP models, this has been shown to yield more accurate Thompson samples whilst incurring only an O((m, + (M)mmM+M)BB))) BN), on top the (t.Bm)) SVGP model fit, per optimization step t.
For our theoretical analysis, we consider two distinct decoupled sampling rules inspired by [24] one for cach of the two SVGP formulations presented above [see 24, for derivations and similar expressions for Fourier decompositions). The first rule is referred to as Decoupled Sampling with Inducing Points and is defined as
M
j((x) = Σ a oay (oVAy100s(()) + Σ u 2k 2z),
(4)
where we have coefficients U1 = [KZz,,, (an(u, mi) + m a, mL m.. A war)l, for 2m... ( (z1), ) ... &M (2m,)) TTT and wM [wi,... wNs The weights w, are drawn i.dd from (/ 1).

---

(4) is modification of the sampling rule of [24] where we have added scaling parameter E IR (with or = 1, the sampling rule [24] is recovered). When set be greater than one, 04 serves to increases the variability of the approximate function samples (without changing their mean) and is used our analysis to ensure sufficient exploration.
To efficiently sample from our second class of SVGP models, we also consider Decoupled Sampling with Inducing Features:
M f((x) = Σ as LavAyuujfee + ==1
Σ Ue
(5)
where Vt. = [Amy (a
me) + m - at a Wmilly for Am, defined in Section 3.3.
5 Regret Analysis of S-GP-TS
Here, we first establish an upper bound on the regret of any approximate GP model (Theorem 1) based the quality of their approximale posterior, parameterized in Assumptions and 4. We then discuss the consequences of Theorem for the regret bounds and the computational complexity of S-GP-TS methods based on SVGPs and the decoupled sampling rules (4) and (5).
5.1 Regret Bounds Based the Quality of Approximations
Consider a TS algorithm using an approximate GP model. In particular, assume an approximate model is provided where k,, ãg and it are approximations of kt, a, and Pt, respectively. At each time i, a batch of B samples ( is drawn from a GP with mean Ao- and the scaled covariance kf The batch of observation points (5t,,)o=1 selected as the maximizers of { over a discretization D, of the search space.
We start our analysis by making two assumptions on the quality of approximations jú, ã of the posterior mean and the standard deviation. This parameterization agnostic the particular sampling rule (governing f4 and ō,) and provides valuable intuition that can be applied to any approximate method. When it comes to S-GP-TS (as the model governing fr.. ot), we show, Sec. 5.2, that these assumptions are satisfied under some conditions on the value of the parameters of the sampling rules.
Assumption 3 (quality of the approximate standard deviation). For the approximate o, the exact ot and for all x X.
1 - [x()) - Et S 5x(x) S (g)((x) + Et.
where 1 S S Sa 1a.lsa de S for all t 2 1 and some constants a, ā E R, and 0 S 6 S for all t and some small constant
Assumption 4 (quality of the approximate prediction). For the approximate it, the exact μr and ot. and for all x E X,
[ī)(x)) - (x(xa) S cor(a))
where 0 S at S cfor all 1 and some constant e E R.
The following Lemma establishes concentration inequality for the approximate statistics using the one for exact statistics [3, Theorem 2].
Lemma 1. Under Assumptions 2, 3 and 4, with probability at least If(x) - [x(x) S
(x( + E)) where (g( = 2. B R2/2(1 + log(1/6)) + +a).
Proof is provided in Appendix B. Here, Ys is the maximal information gain: ?'
MaxxcxX Al= I([v()))zea) =aXAxx//A==<<V(t>>)eeeee. (f()) where [[[v(=))me (f(laaa;a((())))))) [f(x))=<A) denotes the mutual infor-
mation [44, Chapter 2) between observations and the underlying GP model. The maximal information gain can itself be bounded for specific kernel (see Sec. 5.3).
Following [29] and [3], we consider discretization D, of the search space satisfying the following assumption.
5

---

Assumption 5. The discretization D, is designed in a way that If - f(x()) S 1/t2 for all & X, where x( argmin Ir æ'll is the closest point (in Euclidean norm) to in De- The size of this discretization satisfies [ = N, < Cdd, B) where C(d, B) is independent of ([3, 29]).
We are now in a position to present regret bounds based on the quality of GP approximations:
Theorem 1. Consider S-GP-TS with 04 28c(1/(12)) Under Assumptions 1, 2, 3 and 5, the regret defined in (1), satisfies
YT log(1 +)
30 B, R B: &
+ (313r or) (TT + 15BB +
= ( LaBR der(TrB (Tg( (ogg))) + acTBR d(rro + log(T)) log(T)) log(T))oog(T) (6)
where 5.
(bb + } with bt (o
See the proof in Appendix B. This regret bound scales with the product of the ratios a and a, with an additive term depending on the additive approximation error in the standard deviation.
5.2 Approximation Quality of the Decoupled Sampling Rule
For S-GP-TS with inducing points, assume, in [26], that the inducing points sampled according to discrete k-DPP. While this might be costly in practice, [26] showed that Ze can be elficiently sampled from Eo close sampling methods without compromising the predictive quality of SVGP. For both sampling rules, we also assume in our analysis that the Mercer decomposition of the kernel is used.
The quality of the approximation can be characterized using the spectral properties the GP kernel. Let us define the tail mass of eigenvalues وAs = Σ -i=M+1 à where o; maX (z).. With decaying eigenvalues, including sufficient eigenfunctions in the feature representation results in a small d In addition, [26] showed that, for SVGP, sufficient number of inducing variables ensures that the Kullback-Leibler (KL) divergence between the approximate and the true posterior distributions diminishes. Consequently, the approximate posterior mean and the approximate posterior variance converge to the true ones. Building on this result, we are able to prove Proposition I on the quality of approximations.
Proposition 1. For S-GP-TS based on sampling rule (4) with ar = 1 and an SVGP using an to close k-DPP for selecting Ze, with probability at least 1 6, Assumptions 3 and 4 hold with parameters
Q = /R, Qt VI
Ot -1 V3kr and 6 VCim,dss where C₁ is a constant
Vike
specified in the appendix and K₁ 2366m + A18en.
78
For S-GP-TS based sampling rule (5) with ar 1, Assumptions and hold with parameters ā. /1 + BBke and e = VCmm,ōM, where C₁ the same constant
G. = /op a, V1- /ane 21.6mn
as above and Ki.
Note that our proposition requires extending the results of [26] in two non-trivial ways. First, the decoupled sampling rules introduce an additional error. Secondly, [26] built their convergence results on the assumption that the observation points Fr,5 are drawn from prefixed distribution, which not the case S-GP-TS, where T4,6 are selected according to an experimental design method. detailed proof of Proposition is provided in Appendix B.
5.3 Application of Regret Bounds to Matérn and SE Kernels
We investigate the application of Theorem kernels, widely used in practice [see, e.g., 17, 45). (.. In the case of a Matérn kernel with smoothness parameter
the Squared Exponential (SE) and Matérn
it is known that A O اويية [46]. For the SE kernel, have Xラ = O(exp(-ja) )))) [47, 48]. With these bounds on the spectrum of the kernels and the specific bounds on the maximal information gain [c.g., 7s < (log(ss) d+)) for SE and 7s S ((sd/(2rdd) log(s) for Matém,, 49], Theorem and Proposition I result in the following theorem.

---

Theorem 2. Under Assumptions I and 2, with the algorithmic parameters, kernels and sampling rules specified in Table I, S-GP-TS offers R(T, B; f) = YTBBBTTTTTTT log(T)).
With a batch size = Theorem 2 recovers the same regret bounds as the exact GP-TS [3]. We also note that for fair comparison in terms of both the batch size and number of samples we should consider
TB as the number of samples. In that case, our regret bound becomes O (BBBrrr n log(T" (B)(T//)), B)), which scales at most with v vB. That is v B tighter than the trivial scaling with B.
In order to prove Theorem 2, the algorithmic parameters M and my must be selected large enough such that approximation parameters 2, 0, c. in Assumptions 3 and 4 are sufficiently small. Using the relation between the algorithmic parameters, the approximation parameters and m, provided by Proposition 1, the regret bound follows from Theorem 1, See Appendix B for detailed proof.
The values of M and m, required for Theorem 2 are summarized in Table 1. We also show the resulting computational cost of cach sampling rule (as given by O (B(M + mr)NTT + ))), explicitly demonstrating the improvement of S-GP-TS over the O + B) computational cost of the vanilla GP-TS. Note that, for the Matérn kernel under sampling rule (4), required to be sufficiently larger than //2 in order for m. grow slower than
Inducing points (4)
Inducing features (5)
Matérn
Condition
my
T Tu M ~ T (2w-dd)d
24
d2 djy
m2 T M ~ (avid)d
12 BNT 2h op
BNTT 10)2
Cost
+ BT min (T a T
Condition
me M ~ (log(T))²
my, M (log(T))²
SE
Cost
(BNrT (T) + BT (T))
log + BT log (T))
Table 1: Conditions on the number of features m, and inducing variables My required for Theorem alongisde the resulting cost of each decoupled sampling method.
6 Experiments
We now provide an empirical evaluation of S-GP-TS. As [24] have already comprehensively demon- strated the practical advantage of decoupled sampling for problems with small optimization budgets, we focus here on scalability of S-GP-TS, and in particular a) its efficiency with large batch size, b) its ability to handle large data volumes. We first investigate collection of classical synthetic problems for BO, before demonstrating S-GP-TS in a challenging real-world high-throughput molecular design considered by [13]. Our synthetic experiments focus on multi-modal problems with substantial observation noise, as these cannot be solved accurately with a s small budget yet are still unsuitable for local, exhaustive, or deterministic optimization routines. Our implementation is provided as part of the open-source toolbox trieste [50] and relies also on goflow [5] and gpflux [52].
As is often the case, our regret- based analysis applies to version of the algorithm that is slightly different to a practically viable BO method. Rather than focusing on recreating our algorithm exactly calculate eigen-features exactly)
in the very limited settings (e..g 1-d RBF kernel for which that are of little interest to the BO community, we instead choose to demonstrate the practical strength and unprecedented scalability of S-GP-TS by investigating an implementation that could be used by practitioners. The resulting algorithms demonstrated in this section are still well-aligned with our work through their use of sparse GP surrogate models and decoupled Thompson sampling.
6.1 Synthetic Benchmarks
We first consider two toy problems: Hartmann (6 dim, moderately multi-modal) with a large additive noise and Shekel (4 dim, highly multi-modal) with moderate noise, see Appendix C for the full description. Our SVGP models use inducing points and a Matérn kernel with smoothness parameter
Phttps://githin www/secooddindd h labs ttte/ttriteeeeee
7

---

== 2.5. As eigenfunctions for this kernel are limited to small dimensions [39], we implement decoupled TS using the easily accessible random Fourier Features (RFF). Note that [24] have shown decoupled sampling signilicantly alleviate the variance starvation phenomenon (underestimating the variance of points far from the observations [16, 43]) that lypically hampers the efficacy of RFFs. We use M 1000 features and maximise each sample as in (3) using L-BFGS-B, B, starting from the best point among a large sample.
As sampling inducing points from a k-DPP is prohibitively costly for the repeated model fitting required by BO loops, we use the the greedy variance selection method of [26] which is ep close to k-DPP and has been shown to outperform optimisation of inducing points in practice. We also consider the practical alternative of choosing inducing points chosen by k-means clustering of the observations. As the optimisation progresses, observations are likely to be concentrated in the optimal regions, so clustering would result in somehow "targeted" inducing points for BO, In order to control the computational overhead of S-GP-TS and to allow efficient computational implementation (i.e. avoiding Tensorflow recompilation issues), we use a fixed number m, of points, set to either 250 or 500. Similarly, we set the covariance scaling parameter az to avoid having dynamic tunable parameters, like those that plague UCB-based approaches.
For each experiment, we rum t 50 steps of S-GP-TS with B = 100 (i.e. 5,000 total observations). For baselines, we compare against t 750 steps of standard sequential non-batch BO routines with an exact GP model: Expected Improvement [EI, 53], Augmented Expected Improvement [AE], 54), and an extension of Max-value Entropy search suitable for noisy observations [GIBBON, 10]. Due to the large number of steps, we only consider low-cost but high-performance acquisitions, following the cost-benefit analysis of [10], and exclude the popular knowledge gradient [9] or classical entropy search [55, 40J. Popular existing batch acquisition functions do not scale to batches as large as B = 100, however, we present their performance on smaller batches across additional experiments in Appendix C. We report simple regret of the current believed best solution (maximizer of the current model mcan) across the previously queried data points. All results are averaged over 30 runs and reported function of either the number of function evaluations Bt for S- S-GP-TS and for the baselines), or the number of BO iterations, in Figure 1.
SGPTS
6x 10
10a
Greedy
M=250
10
GP
3x 10
Greedy
2 10
M=500
GPG
1000
2000
3000
4000
10
1000
2000
Kmeans
Objective Function Evaluations
# Objective Function Evaluations
M=250
SGGTS
Kmmeas
6 10
10a
M=500
108
3x 10a
AEI
GIBBON
2x 10
10
300
700
300
# BO Iterations
90 Rerations
Figure 1: Simple regret on Shekel (4D, left) and Hartmann (6D, right). When considering regret with respect the total number of objective function evaluations B (top panels, purely sequential setting), all S-GP-TS methods are initially less efficient (Shekel) or match the performance (Hartmann) of the best baselines, however the best S-GP-TS approach is able efficiently allocate its additional budget to achieve lower final regret. When considering regret with respect to the BO iteration (bottom panels, idealised parallel setting), S-GP-TS achieves low regret in a fraction of the iterations required by the standard BO routines.
The fact that S-GP-TS is able to find solutions on both benchmarks with substantially improved regret than found by standard BO, provides strong evidence that S-GP-TS effectively leveraging parallel resources. Moreover, as these higher-quality solutions were only found after large number of total evaluations, Figure also highlights the necessity for BO routines, like S-GP-TS, that can handle these larger (heavily parallelized) optimization budgets. We reiterate that the existing BO baselines cannot handle as many evaluations as S-GP-TS, becoming prohibitively slow once we surpass 750 data-points), When considering the regret achieved per individual function evaluation, expect batch routines to be less efficient than purely sequential BO routines. However, in the case of the Hartmann function (the benchmark with the largest observation noise), see that our best
typically
8

---

S-GP-TS exactly matches (before going to exceed) the performance of the sequential routines, suggesting that S-GP-TS a particularly effective optimizer for functions with significant levels of observation noise.
Note that the performance of S-GP-TS is sensitive to its chosen inducing points, with k-means providing the most effective routines. On Hartmann, 250 inducing points is sufficient to deliver good performances, while on Shekel, which much more multimodal, using a larger number is critical.
6.2 High-throughput Molecular Search
Finally, we investigate the performance of S-GP-TS with respect to an established baseline for high-throughput molecular screening. Although molecular search has been tackled many times with BO [56, 22, 57], only the approach of [13] standard (non-decoupled) TS over Bayesian neural network (BNN-TS) is truly scalable. We now recreate the largest experiment considered by [13], where the objective is uncover the top 10% of molecules in terms of power conversion efficiency among a library of 2.3 million candidate from the Harvard Clean Energy Project [58]. Molecules are encoded as Morgan circular fingerprints of Bond radius 3 (i.e. 512-dimmnni bit vectors, see [59]).
As the standard GP kernels considered above are not suitable for sparse and high-dimensional molecule inputs [60], we instead build our SVGP with zeroth order ArcCosine kernel [61], chosen due to its strong empirical performance under sparsity and as it permits a random decomposition that can be exploited to perform decoupled TS. In particular, we use the M-feature decomposition investigated by [62] of
Luf xx)) = 2 / dw (2 O(w ())("))""x
2 M
Σ [(w]x)0(w))x)
where (()) is the Heaviside step function and w ~ N( 1). In our experiments,
M = 1 000 random features and, avoid memory issues, compute our GP samples over a random subset of 100 molecules (renewed at each sample). We run S-GP-TS twice, once with me 500 and once with 2000 inducing points. We chose inducing points as uniform samples from the already evaluated molecules (for each model step). as preliminary experiments showed that neither the k-means nor greedy selection routines discussed above were effective when applied to sparse and high-dimensional molecular fingerprint inputs.
Following [13], report the recall (fraction of the top 10% of molecules far chosen by the BO loop) for S-GP-TS, along with the performance of BNN-TS, a greedy BNN (that queries the B maximizers of the BNN's posterior mean), and random search baseline (all taken from [13)). All routines (including our S-GP-TS) are ran for = 250 successive batches of B = 500 molecules. Figure 2 shows that S-GP-TS is able to perform effective batch optimization over very large optimization budgets (120,000 total evaluations) and, when using m 2000 or even just m 500 inducing points, S-GP-TS matches the performance of [13]'s BNN-based TS and greedy sampling approaches, respectively. Note that due to the high computational demands of this experiment, we report just single replication of S-GP-TS (a limitation also of [13]'s 's results). However, we stress that an additional realization of the m = 500 experiment returned indistinguishable results.
7 Discussion
We have shown that S-GP-TS enjoys the same regret order as exact GP-TS but with greatly reduced O M)) computation per step f, compared to the O(N2) cost of the standard sampling. However, the discretization size N, is exponential in the dimension d of the search space and so remains a limiting computational factor when optimizing over high dimensional search spaces. Hence, while S-GP-TS with decoupled sampling rule allows orders of magnitude larger optimization budgets compared vanilla GP-TS, it still suffers from the curse of dimensionality. Intuitively, this seems inevitable due to NP-Hardness of non-convex optimization problems [see, e.g., 63] as required to find the maximizer of the GP-UCB acquisition function [see, e.g., 30], or even in the application of UCB to linear bandits [64]. In particular, the computational cost of the ofe adaptive sketching method for implementing GP-UCB [30] was reported as O ONr where dett, referred to as the effective dimension of the problem, is upper bounded by YT.

---

0.45
BNN-TS
0.40
S-GP-TS m=2000
0.35
BNN-Greedy
10% Recall
0.30
S-GP-TS m=500
Random Search
0.25
0.20
0.15
0.10
0.05
0.00
20000
40000
60000
80000
100000
120000
# Evaluations (B=500)
Figure 2: Proportion of the top 10% of molecules found by each of the search routines. S-GP-TS is able to process substantial data volumes and effectively allocates large batches, matching the performance of the well-established BNN baselines.
An important practical consideration when using S-GP-TS in practice is how choose its inducing points. The performance improvement provided by choosing inducing points by k-means rather than greedy variance selection, as demonstrated in our experiments, raises the possibility that BO-specific routines for choosing inducing points could allow even better performance. This is an important avenue for future work.
References
[1] William Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 1933.
[2] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas. Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 2016.
[3] Sayak Ray Chowdhury and Aditya Gopalan. On kernelized multi-armed bandits. In International Conference Machine Learning, 2017.
[4] Ricardo Baptista and Mathhis Poloczek. Bayesian optimization of combinatorial structures. In International Conference on Machine Learning. 2018.
[5] David Eriksson and Matthias Poloczek. Scalable constrained bayesian optimization. In Interna- tional Conference on Artificial Intelligncce and Statistics, 2021.
[6] Samuel Daulton, Shaun Singh, Vashist Avadhanula, Drew Dimmcry, and Eytan Bakshy. Thomp- son sampling for contextual bandit problems with auxiliary safety constraints. arXiv preprint arXiv: 1900633, 2019.
[7] Clément Chevalier and David Ginsbourger. Fast computation of the multi multi-points expected improvement with applications in batch selection. In International Conference on Learning and Intelligent Optimization, 2013.
[8] Javier González, Zhenwen Dai, Philipp Hennig, and Neil Lawrence. Batch Bayesian optimiza- tion via local penalization. In Artificial intelligence and statistics, 2016.
[9] Jian Wu and Peter I Frazier. The parallel knowledge gradient method for batch Bayesian optimization. In Advances Neural Information Processing Systems, 2016.
[10] Henry B Moss, David S Leslie, Javier Gonzalez, and Paul Rayson. Gibbon: General-purpose information-bart based bayesian optimisation. arXiv preprint arXiv:2102.03324, 2021.
[11] Hamed Jalali, Inneke Van Nieuwenhuyse, and Victor Picheny. Comparison of kriging-based based algorithms for simulation optimization with heterogeneous noise. European Journal of Opera- tional Research, 2017.
10

---

[12] Mickaël Binois, Jiangeng Huang, Robert B Gramacy, and Mike Ludkovski. Replication or exploration? sequential design for stochastic simulation experiments. Technometrica, 2019.
[13] José Miguel Hernández-Lobato, James Requeima, Edward O Pyzer-Knapp. and Alán Aspuru- Guzik. Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space. In International Conference on Machine Leurning, 2017.
[14] Kirthevasan Kandasamy, Akshay Krishnamurthy, Jeff Schneider, and Barnabás Póczos. Paral- lelised Bayesian optimisation via Thompson sampling. In International Conference on Artificial Intelligence and Statistics, 2018.
[15] Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable Bayesian optimization using deep neural networks. In International conference on machine learning, pages 2171-2180, 2015.
[16] Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale Bayesian optimization in high-dimensional spaces. In International Conference on Artificial Intelligence and Statistics, 2018.
[17] Carl E Rasmussen and Christopher KI Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
[18] Peter Diggle, Jonathan A Tawn, and Rana A Moyeed. Model-based geostatistics. Journal of the Royal Statistical Society: Series C, 1998.
[19] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, 2019.
[20] M. K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics, 2009.
[21] Mitchell McIntire, Daniel Ratner, and Stefano Ermon. Sparse gaussian processes for bayesian optimization. In Associateion for Uncertainty in Artificial Intelligence, 2016.
[22] Ryan-Rhys Griffiths and José Miguel Hernández-Lobato. Constrained Bayesian optimization for automatic chemical design. arXiv preprint arXiv:1709.05001 2017.
[23] Ang Yang, Cheng Li, Santu Rana, Sunil Gupta, and Svetha Venkatesh. Sparse spectrum Gaussian process for Bayesian optimization. arXiv preprint arXiv: /006::990008888 08898, 2019.
[24] James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Peter Deisenroth. Efficiently sampling functions from Gaussian process posteriors. International Conference on Machine Learning, 2020.
[25] My Phan, Yasin Abbasi Yadkori, and Justin Domke. Thompson sampling and approximate inference. In Advances in Neural Information Processing Systems, 2019.
[26] David Burt, Carl Edward Rasmussen, and Mark Van Der Wilk. Rates of convergence for sparse variational Gaussian process regression. In International Conference on Machine Learning, 2019.
[27] Daniel J. Russo, Benjamin Van Roy, Abbas Kazerouni, Ian Osband, and Zheng Wen, A tutorial on thompson sampling, Foundational Trends in Machine Learning, 2018
[28] D. Russo and B. Van Roy. Learning to optimize via posterior sampling. Mathematics of Operations Research, 2014.
[29] Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger. Gaussian process opti- mization in the bandit setting: no regret and experimental design. In International Conference on Machine Learning, 2010.
[30] Daniele Calandriello, Luigi Carratino, Alessandro Lazaric, Michal Valko, and Lorenzo Rosasco. Gaussian process optimization with adaptive sketching: scalable and no regret. In Conference on Learning Theory, 2019.
11

---

[31] Alain Berlinet and Christine Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011.
[32] James Hensman, Nicoló Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence (UAI 2013), 2013.
[33] James Hensman, Nicolas Durrande, Amo Solin, et al. Variational Fourier features for Gaussian processes. Journal of Machine Learning Research, 2017.
[34] Vincent Dutordoir, Nicolas Durrande, and James Hensman. Sparse Gaussian Processes with Spherical Harmonic Features. In Proceedings of the 37th International Conference on Machine Learning, 2020.
[35] Miguel Lázaro-Gredille and Anibal Figueiras- Vidal. Inter-domain Gaussian processes for sparse inference using inducing features. In Advances in Neural Information Processing Systems, 2009.
[36] Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaus- sian processes and kernel methods: A review on connections and equivalences. arXiv preprint 1807.08807.02588 2018.
[37] Viacheslav Borovitskiy, Alexander Teremin. Peter Mostowsky, and Marc P. Deisenroth. Matern Gaussian processes on Riemannian manifolds. In Advances Neural Information Processing Systems, 2020.
[38] Huaiyu Zhu, Christopher KI Williams, Richard Rohwer, and Michal Moreiniec. Gaussian regression and optimal finite dimensional linear models, 1997.
[39] Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 2020.
[40] José Miguel Hernández-Lobato. Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Advances in Neural Information Processing Systems, 2014.
[41] José Miguel entropy search for efficient global optimization of black-box functions. In Advances in Neural Information Processing Systems, 2014.
Matthew W Hoffman, and Zoubin Ghahramani. Predictive
[42] Salomon Bochner et al. Lectures on Fourier integrals. Princcton University Press, 1959.
[43] Mojmir Mutny and Andreas Krause. Efficient high dimensional Bayesian optimization with additivity and quadrature fourier features. In Advances in Neural Information Processing Systems 31, 2018
[44] Thomas M Cover. Elements of information theory. John Wiley Sons, 1999.
[45] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. In Advances Neural Information Processing Systems 25, 2012.
[46] Gabriele Santin and Robert Schaback. Approximation of cigenfunctions in kernel spaces. Advances in Computational Mathematics, 2016.
[47] Mikhail Belkin, Approximation beats concentration? an approximation view on inference with smooth radial kernels. In Conference On Learning Theory, 2018.
[48] Gabriel Rayyol,-Mayoll, Paul-Christian Burkner, Michael R. Andersen, Amo Solin, and Aki Vehtari. Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. arXiv preprint 2rXiv: /408, 2020.
[49] Sattar Vakili, Kia Khezeli, and Victor Picheny, On information gain and regret bounds in Gaussian process bandits, In International Conference on Artificial Intelligence and Statistics, 1.21.
12

---

[50] Joel Berkeley, Henry B. Moss, Artem Artemev, Sergio Pascual-Diaz, Uri Granta, Hrvoje Stojic, Ivo Couckuyt, Jixiang Quing, Loka Satrio, and Victor Picheny. Trieste, 2021.
[51] Alexander G de G Matthews, Mark van der Wilk, Tom Nickson, Keisuke Fujii, Alexis Bouk- ouvalas, Pablo León-Villagee Villagrá, Zoubin Ghahramani, and James Hensman. Gpflow: A gaussian process library using tensorflow. Journal of Machine Learning Research, 2017.
[52] Vincent Dutordoir, Hugh Salimbeni, Eric Hambro, John McLeod, Felix Leibfried, Artem Artemev, Mark van der Wilk, James Hensman, Mare P Deisenroth, and ST John. Gpflux: A library for deep Gaussian processes. arXiv preprint arXiv:2
2021.
[53] Donald R Jones, Matthias Schonlau, and William Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 1998.
[54] Deng Huang, Theodore T Allen, William Notz, and Ning Zeng. Global optimization of stochastic black-box systems via sequential kriging meta-models. Journal of global optimization, 2006.
[55] Philipp Hennig and Christian J Schuler. Entropy search for information-efffiient global opti- mization. Journal of Machine Learning Research, 2012.
[56] Rafael
Jorge Aguilera-Iparraguirre, Timothy D Hirzel, David Duvenaud, Dougal Maclaurin, Martin A Blood--Foryythe, Hyun Sik Chac, Markus Einzinger, Dong-Gwang Ha, Tony Wu, al. Design of efficient molecular organic fight-emiting diodes by a high- throughput virtual screening and experimental approach. Nature Materials, 2016.
[57] Henry B Moss, Daniel Beck, Javier González, David S Leslie, and Paul Rayson. Boss: Bayesian optimization over string spaces. Advances in Neural Information Processing Systems, 2020.
[58] Johannes Hachmann, Roberto Olivares Amaya, Sule Mahann-vrenne Carlos Amador- Bedool,, Roel S Sánchez-Carrera, Arych Gold-Parker Leslic Vogt, Anna M Brockway, and Alán Aspuru- Guzik. The Harvard clean energy project: large-scale computational screening and design of organic photovoltaics on the world community grid. The Journal of Physical Chemistry Letters, 2011.
[59] David Rogers and Mathew Hahn. Extended-commectivity fingerprints. Journal of chemical information and modeling, 2010.
[60] Henry Moss and Ryan-Rhys Griffiths. Gaussian process molecule property prediction with flowmo. Advances in Neural Information Processing Systems: Workshop on Machine Learning for Molecules., 2020.
[61] Youngmin Cho. Kernel methods for deep learning. PhD thesis, UC San Diego, 2012.
[62] Kurt Cutajar, Edwin V Bonilla, Pietro Michiardi, and Maurizio Filippone. Random feature expansions for deep Gaussian processes. In International Conference on Machine Learning, 2017.
[63] Prateek Jain and Purushottam Kar. Non-convex optimization for machine learning, Foundational Trends in Machine Learning, 2017.
[64] V. Dani, T. P. Hayes, and M. Kakade. Stochastic linear optimization under bandit feedback. In Conference on Learning Theory, 2008
[65] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[66] Rodolphe Le Riche and Victor Picheny. Revisiting Bayesian optimization in the light of the coco benchmark. arXiv preprint 1rXiv:2103.16644 2021.
[67] Dong Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical programming, 1989.
13


Reviews:
[{"review_id": "ylZ41X_V8QH", "reviewer": "Reviewer_AeT2", "summary": "This work considers Bayesian optimization associated with a Thompson sampling scheme, where the posterior over the function is represented as a Gaussian Process (GPs). Due to the poor scaling of GPs, a pseudo-input approximation of the GP is proposed and incorporated into the bandit optimization procedure. Regret bounds are provided for the proposed scheme which establish accuracy/complexity tradeoffs. Experiments demonstrate the merits fo the proposed approach.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "\nThe paper is well-written, puts forth a novel algorithm, and contributes interesting convergence theory to substantiate the proposed technique. This is a worthy contribution to the intersection of sparse Gaussian Processes and bandit optimization. \n\nIn discussing the use of sparse Gaussian Processes in Bayesian settings, it should be noted that methods other than Thompson sampling have been previously studied, which are very similar to this work, based upon, e.g., Upper-confidence bound or expected improvement:\n\nCalandriello, D., Carratino, L., Lazaric, A., Valko, M., & Rosasco, L. (2020, November). Near-linear time Gaussian process optimization with adaptive batching and resparsification. In International Conference on Machine Learning (pp. 1295-1305). PMLR.\n\nBedi, A. S., Peddireddy, D., Aggarwal, V., & Koppel, A. (2020). Sublinear Regret and Belief Complexity in Gaussian Process Bandits via Information Thresholding. arXiv preprint arXiv:2003.10550.\n\nWhere explicit tradeoffs between regret and the complexity of the GP are provided. I'd like to see a quantitative contrast of Theorem 1 with the results presented in the above works, especially because the later reference incorporates the selection of inducing inputs directly via the notion of information-theoretic regret, whereas in the work under review, the inducing input selection must be done from a determinant point process that is untethered from the underlying optimization problem. Put another way, is there any relationship between the selected pseudo-inputs in this work and the optimal information gain defined over the collection of $m_t$ points?\n\nThe analysis offers a novel incorporation of the approximation error incurred by inducing inputs into Thompson sampling techniques for Gaussian Processes. While this is a nice confluence of ideas, I am a little unclear on what is completely conceptually innovative here in terms of the algorithm design/convergence analysis.\n\nExperimentally, it is also conspicuous that there is no comparison to reference [30]. How does Thompson sampling alleviate limitations of existing works based on other acquisition functions? The experiments do not demonstrate the proposed approach yields a state of the art tradeoff in regret/posterior complexity. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The submission is sufficient along this axis.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "yIYhkXC6Jf2", "reviewer": "Reviewer_Fhfj", "summary": "This paper studies the use of Sparse Gaussian Process regression in combination with Thompson sampling (S-GP-TS) for sequential (batch) decision processes.\n\nVia theoretical analysis, the authors show that batch S-GP-TS can achieve the same regret order as an exact GP-TS algorithm, as long as the GP approximation technique chosen meets certain criteria.\n\nThey argue that approximate S-GP-TS provides computational benefits, with no regret performance loss over exact GP-TS.\n\nThe proposed S-GP-TS is shown to provide good regret results in both synthetic and a real-world Bayesian optimization tasks.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The main contribution of this work is to show that approximate S-GP-TS provides computational benefits, with no regret performance loss, when compared to exact GP-TS.\n\nEven though the originality is limited (S-GP-TS has been already proposed in [24] and the proof technique leverages work in [26]), the significance of this work lies on providing a theoretical regret bound for the general GP approximation case (Theorem 1), as well as a more specific bound (Theorem 2) for the approximation choices of S-GP-TS (as presented in Section 4).\n\nThe paper is clearly written, with a concise introduction and problem formulation, a very clear description of GPs with their corresponding approximations in Section 3 and the S-GP-TS approach in Section 4.\n\nSection 5 is where the main contributions lies (Theorems 1 and 2), with experiments provided in Section 6.\n\nThe theoretical work in Section 5 leverages already presented rates of convergence for sparse variational Gaussian process regression in [26], but accounts for\n1. the errors introduced by S-GP-TS decoupled sampling rules, and\n2. the sequential selection of the observation points $x_{t,b}$ of S-GP-TS.\n\nThe provided bounds depend on:\n1. the approximation error of Assumptions 3 and 4,\n2. the discretization grid of Assumption 5.\n\nThese errors are made specific for S-GP-TS, based on the spectrum properties of selected GP kernels and the number of prior features and inducing variables required. The authors argue that their bound shows no order regret of S-GP-TS when compared to the exact case if such assumptions hold.\n\nTo counteract the approximation error of Assumptions 3 and 4, the number of features $m_t$ and $M_T$ have an exponential dependency over dimensionality (as presented in Table 1):\n- Exponentially many points are needed to provide an approximation to the exact GP with no regret.\n- The authors argue that the discretization size being exponential in the dimension $d$ is inherent when optimizing over high dimensional spaces, and not specific to the approximation.\n\nImportantly, the authors show that the approximation techniques (with enough inducing points/features) allow for a computational cost reduction (which is significant in many practical Bayesian optimization techniques).\n- As a limitation, they acknowledge that studying the selection mechanism of inducing points/features is an open avenue of work.\n\nSome clarifications on the applicability of their results would help the unfamiliar reader:\n- $\\tau$ first appears in Equation (2), as the Gaussian observation noise variance or precision.\n- However, it also appears in the constants of Theorem 1 and Proposition 1.\n- Therefore, are these results only applicable to Gaussian noise observation models? I.e., for exact Gaussian regression models?\n    * [Updated] I encourage the authors to incorporate the informative discussion along these lines from their rebuttal to the main manuscript\n\nFinally, they provide empirical results for optimization problems with large batch sizes in Section 6.\n- The significance of this section is a bit unclear as they focus on evaluating regret for large batch sizes, but don't provide computational complexity comparisons with exact TS-GP.\n- For their synthetic benchmarks, they state that S-GP-TS has access to 5000 observations (t=50 times B=100) while the baselines are run only for t=750 steps. Why are the authors comparing algorithms with different number of observations? If the alternatives can not deal with these many observations, computational results (i.e., run-time) for the same number of observations could be shown to better illustrate the claimed benefits (same reward, less computational burden).\n    * [Updated] The explanations provided in the author responses helps understand the subtleties of the experimental section better, which if incorporated into the main manuscript, would clarify the significance of the presented results.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discuss the general limitations (e.g., dimensionality dependence) of the presented work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "PDgom7vjlx7", "reviewer": "Reviewer_MHRh", "summary": "This paper provides a regret bound for Thompson sampling of SVGP posterior. Explicitly, this paper first derives a generic regret bound for any arbitrary approximate posterior of a certain uniform closeness (i.e., approximation quality) to the full posterior. This bound is then applied to two different sampling schemes (for two respective variants of SVGP). \n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "\nI find the paper very well written and clearly organized. There might be a small typo in Lemma 1, where the wrong symbol of B is used (either that or the typo is in the proof of Lemma 1).\n\nI think the theoretical analysis is interesting and quite different from the theory of GP-TS. It has also provided good insights to the high level challenges of extending the regret bound of GP-TS to S-GP-TS. As such, the resulting theory is a nice complement to S-GP-TS. On the other hand, while I appreciate the extensive empirical study on many real-world/synthetic benchmarks, it feels somewhat loosely connected to the theory:\n\n1. The entire point of S-GP-TS is to scale up GP-TS, but there is no runtime plot to showcase the difference. I don't think performance vs #query and performance vs #iterations are sufficient to demonstrate this point.\n2. There seems to be no experiment to backup the analysis of SVGP with inducing features.\n3. Molecular search experiment uses an entirely different kernel which has little relevance to the theory developed.\n4. No ablation study to verify other aspects of the theory (e.g., approximation quality in Proposition 1, validity of various assumptions).\n\nAlso, regarding proof of Theorem 1: It is claimed that the approximation error in the standard dev might result in an under exploration which the original GP-TS analysis cannot address. Does this mean the same proof technique will not go through, or will it get a significantly worse regret bound because of this term, and if so, how much worse?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have not discussed the broader impact of this work, but I believe there is no potential negative impact to society.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "GMi3mmKyIf6", "reviewer": "Reviewer_MheU", "summary": "This paper analyzes Thompson Sampling for sequential optimization where the GP is replaced by a scalable approximation. In particular, leveraging prior works on Sparse Variational GPs and ways of sampling from and updating Sparse Variational GPs, the authors show that the resulting sparse-GP-TS algorithm (which makes use of these approximations) has $\\tilde O(\\sqrt(T))$ regret. In particular, the authors first prove a general regret bound for Thompson Sampling with approximate GPs that satisfy certain error bounds on the means and covariances. They then show that SVGPs satisfy these bounds. The authors then use their algorithm and show performance on synthetic and real problems when compared to existing BO methods.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "This paper uses Sparse Variational GPs in Thompson Sampling to allow GP-TS to scale to larger dimensional problems. The proposed algorithm makes extensive use of recent work on Sparse Variational GPs and sampling from Sparse Variational GPs. The analysis makes use of recent work on characterizing the approximation power of SVGPs which---importantly-- the authors show to result in a computationally tractable algorithm.\n\nThis paper is well written with relatively few and addresses an important limitation of GP-TS, namely its scalablility and computational tractability. Furthermore, it provides frequentist regret bounds for TS with a class of approximations to GPs (which is shown to include SVGPs), and extensive numerical simulations showing the impressive performance of their algorithm against benchmarks. While this paper really draws on results from prior works on SVGPs, it extends upon them and applies the findings to the important problem of Thompson sampling, so I think the contribution of this paper is good enough for acceptance.\n\nMy main question is on the cost of choosing inducing points. Choosing inducing points is said to be efficient, but shouldn't this computational cost be included in the computational cost of the TS algorithm for a proper apples-to-apples comparison?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do discuss the limitations of their work, and do a good job of comparing to prior works on GP-TS. I'm curious about how the computation of the inducing points factors into the computational complexity of the algorithm, since solving e.g., K-means at every step may be computationally expensive. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This was a somewhat borderline paper. The reviewers found the theory to be slightly derivative from earlier work, but that it was nontrivial and ultimately had merit. The reviewers also find the empirical results to still be somewhat limited. The reviewers certainly did appreciate the dialogue with the authors during the review process. Authors are strongly encouraged to read the reviews and responses, and ensure that all feedback is taken into account when preparing the camera ready submission.
Decision: Accept (Poster)