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 sightthoogghout 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
R(T, B: E
-
"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 ( 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(s,,,,. 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, = - + @ 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, 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 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()) - -t 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 a. S cfoo 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 (ffl:aa;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 20c(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
G. = /op a, V1- 13ms 21.6mn
ā. /1 + BBke and e = VCmm,ōM, where C₁ the same constant
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) = YTBBB 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 nT 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 of
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 Dutordorr, 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": "iEGWqfMwjU7", "reviewer": "Reviewer_eTmV", "summary": "The paper present a new approximate Bandit/Bayesian optimization algorithm named S-GP-TS based on Thompson sampling with Gaussian process (GP-TS). GP-TS is known to achieve good regret, it is also highly computationally inefficient. In particular it has two bottlenecks:\n- computing the posterior on T historical points requires at least $O(T^2)$ and is not scalable\n- even given the posterior, drawing a sample and evaluating it on N test points requires $O(N^2)$. This is true even if $T=0$, i.e. even evaluating an exact GP draw from prior on N points is computationally hard.\n\nThe paper leverages some old and some new results on speeding up these two bottlenecks to speed up the overall S-GP-TS algorithm:\n- it replaces slow exact inference of the posterior with a sparse variational approximation, in two variants based on random features or inducing points\n- it replaces slow exact sampling from the posterior with a recently proposed technique called decomposed sampling\n\nBoth of these are known techniques, but the authors take care in correctly matching the right variational approximation form and parameters with the right decomposed sampling form and parameters. This is not difficult but it is important and well explained in the paper. They also add batching in the most straightforward way for TS (draw a batch of B samples).\n\nOnce they have derived their algorithm, they proceed to provide sufficient conditions (mostly on the number and distribution of the random feature/inducing points) to guarantee that the algorithm achieve low regret. This analysis is essentially a direct combination of three existing analysis:\n- \"[3] Chowdhury and Gopalan\" for the overall regret bounding framework, which adds the necessity of introducing a discretization of the action space\n- \"[26] Burt et al.\" for the accuracy bound of the sparse variational posterior. Normally this results could not be used as they do not hold for adaptive sequences, but since the authors are already operating on a fixed discretization following [3], they smartly take advantage of this fact to be able to apply the results of [24]\n- \"[24] Wilson et al.\" for the decomposed sampling, here once again thanks to the discretization the application is straightforward\n\nThe only extra technical issue that the authors resolve is handling additive error in their posterior estimate that was not considered often in the literature. However thanks to the sufficiently fast converging rate of the approximation this additive error does not require new tools, only a careful bookkeeping.\n\nOverall the algorithm achieves the same rate as GP-TS. The only downside is the dependency on the batch size which is essentially linear (i.e. it would be almost equivalent to B instances operating in parallel and not sharing information).\n\nThe runtime is greatly reduced although the final complexity table (Table 1) could be accompanied by more explainations (more in the main review).\n\nFinally the authors evaluate the methods empirically on large scale artificial and real world datasets. However the tuning of all the S-GP-TS is completely decoupled from theory, making it a less convincing evaluation. The authors also report simple regret instead of cumulative regret (the focus of the paper) and only report iterations and not wall-clock runtime.", "questions": "", "limitations": "", "rating": 4, "confidence": 5, "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": "Overall I like the idea of the paper as combining two solid and popular techniques (GP-TS and SVGP) with the newly introduced decoupled sampling. However the theoretical and algorithmic part seems almost completely incremental, with no clearly highlighted contribution over the secondary one of handling of additive error.\n\nThis could have been compensated with a good experimental evaluation, but although the authors positively focused on the latest algorithm and real-world, large scale dataset they also diverge significantly from the main premise of the paper: reducing wall-clock time (at least for the GP part) while guaranteeing small cumulative regret. The experiments do not report runtime nor cumulative regret, and are run with parameters that do not guarantee low regret.\n\nOverall this looks like a highly polished, but incremental paper. The theoretical impact is medium, but since it is disconnected by the impact evaluated by the experiments it falls under the threshold.\n\nMore concrete ways in which the theoretical part of the paper could be improved:\n1.  if I understood correctly, the discretization is not only necessary to guarantee correct maximization of the TS draw. Even in a domain where optimization is easy (e.g. linear kernel or 1d) the discretization needs to be there to make the selected points non-adaptive. Since the discretization scales as $t^{2d}$, the linear $N_t$ cost by itself would effectively negate any improvement due to the sparse approximation. It would be maybe better to explicitely target a finite action space (e.g. the library of molecules used in 6.2 is finite and fixed). Would this be enough to avoid having to concretely construct the discretization to run the algorithm?\n2. Thm. 1 establishes a required value for $\\alpha_t$, Proposition 1 is analyzed under $\\alpha_t = 1$. This is not rigorous and should be corrected. It is also not clear how the authors plan to compute a sufficiently accurate value of $\\gamma_t$ to be used in $\\alpha$ with theoretical guarantees.\n3. The whole paper is extremely focused on GP-TS, and almost completely ignores other important baselines in BO. Theoretically, at least UCB based algorithms, which actually achieve a lower regret and potentially runtime than S-GP-TS, should be mentioned and compared in terms of regret and runtime. Currently, they are only mentioned once in the conclusions. The experiments at least have EI, AEI and GIBBON, but again UCB is at least as popular as TS and should be included. Scalability is not an issue since as the authors comment, \"[30] Calandriello et al.\" seems as scalable as S-GP-TS, and there seems to be a [follow-up version capable of batching available](https://github.com/luigicarratino/batch-bkb).\n4. Table 1 is quite unreadable, and again it is not compared to any other scalable GP optimization method. Even not considering the lack of baselines, in most of the paper the authors use maximum information gain to express complexities, while here everything is expanded (e.g. the first $m_t$ seems to be $\\gamma_t^2$ but it is not unclear. It would be good to 4.A) uniform somehow the various complexity so that they can be compared (e.g. it is hard to compare conditions on M for inducing point and features for the matern kernel), 4.B) add remarks explaining rates and pros/cons of each alternative, and 4.C) add a comparison with existing BO baselines.\n\nThe experiments also could be greatly improved, both to more closely match the theory and in general to be made a more rigorous evaluation:\n1) $\\alpha_t$ in the experiment should be set according to Thm. 1 and not set to 1. $m_t$ and $M$ should also try to satisfy the values dictated by theory. If not, they should be cross validated and all curves (also the sub-optimal ones) reported. Cross validation is not usually possible in real-world BO, but toy experiments are an ideal place for this kind of explorations. \n2) Cumulative regret is the goal of the paper and it is not reported.\n3) Wall clock or some form of similar evaluation should be reported. In general the golden rule would be to report your own computational resources and wall-clock (as recommended by the checklist), and take the baselines numbers from the original papers. In case the original papers use different dataset or do not report numbers, it is worth to reimplement and run at least one for comparison. If all else fails, one can substitute per-step theoretical complexities when wall-clock is not available.\n4) Why do you compare BO iterations and function evaluations without including at least another batched baseline? Usually these quantities are compared to show that a batched algorithm does not have much more regret than a sequential one. If instead the goal is to compare convergence speed then the baselines should also be batched (as they are in the appendix). The procedure used to choose the batch size has also not been reported (cross-validated?).\n5) In the main paper experiments S-GP-TS has access to almost 5k evaluation, and baselines to 750. Without wall-clock constraint this seems completely arbitrary, and with enough evaluations the baselines might outperform S-GP-TS.\n\nMinor comments:\n- Thm. 1 is more precise, and a $O(\\sqrt{d})$ factor should also be present in Thm. 2 ($\\gamma_T$ which is not a constant depends on it)\n- recent $k$-DPP sampling algorithms (e.g. [the intermediate samplers from DPPy](https://dppy.readthedocs.io/en/latest/finite_dpps/exact_sampling.html#id26)) are not slow, and they are guaranteed to be roughly as efficient as S-GP-TS, at least for moderate $k$. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors discuss a bit the main limitations of the algorithms in the end, but these limitations should have been considered more throughout the paper (the necessity of the discretization mostly). The open question on how to choose inducing point is also not very in-depth (it misses developments in $k$-DPP sampling from 2019).\n\nBO can have significant societal impact, but as this is mostly a theory paper it is ok to not discuss them.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "9"}, "extra_text": {}}, {"review_id": "dzpWQpuAZJw", "reviewer": "Reviewer_sgAN", "summary": "The paper provides regret bounds for Thompson sampling based on scalable Gaussian processes. Under certain conditions, scalable GP Thompson sampling enjoys the same regret bound as the exact GP Thompson sampling.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "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 work presents regret bounds for Thompson sampling based on scalable Gaussian processes (i.e. with the inducing point approximation) that are sampled with decoupled sampling (i.e. incurring yet another approximation). There is a proof for the given regret bounds that seems solid. The work is original and cites related work adequately.\n\nThe paper is clearly written, the mathematical details being hidden in the appendix. I have not checked all the details. \n\nThe results are important because the main conclusion is that scalable GP-based Thompson sampling enjoys, under some conditions, the same regret bound as exact GP Thompson sampling. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors address the limitations of the scalable GP Thompson sampling, the main of those being the curse of dimensionality. The authors do not address potential societal impact. As the work is largely a theoretical analysis, the societal impact will depend on practitioners applying scalable Thomspon sampling to real-world problems. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "BjwpLuw-Wgw", "reviewer": "Reviewer_Psad", "summary": "The work is dedicated to a regret vs computational complexity analysis of scalable Thompson sampling technique based on the ideas of Wilson et al, ICML 2020. Theoretical bounds of regret for the widely used Matern and squared exponential kernels are derived and the computational complexity is estimated. The conclusion is drawn that the scalable Thompson sampling technique achieves the same order of regret asymptotics as exact Thomspon sampling with substantially reduced complexity. Empirical tests on several benchmark functions and one real world problem are performed where the scalable Thompson sampling technique turns out to be on par with existing baselines.", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "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": "Theoretical results of the paper are interesting. They support the new scalable approximate Thompson sampling technique by confirming that in some particular regime (i.e. number of features, number of inducing variables) it attains the regret bound of the exact analog. The empirical study is rather scarce. Considering more benchmark functions would be beneficial as well as benchmarking runtimes and quantifying uncertainty for the real world problem. Although the delicacy of these is understandable because of implementation dependency of runtimes and the reported high computational demand of the experiment that hinders uncertainty quantification. Another point to mention is that theoretical analysis and the studied empirical setting do not agree: the former uses Mercer theorem expansions and minimizes acquisitions over a discrete set while the latter uses random Fourier features and minimizes acquisitions by a (repeated) gradient descent. This illustrates the fact that the theoretical setting studied differs from the setting that actually usually arises in applications and somewhat impairs the significance of this work.\n\n**Warning**: due to time constraints and large number of papers to review I was unable to check the proofs and calculations of the paper. Since those are the main contributions I write this review assuming those are correct but retaining a very low level of confidence.\n\n# Originality\nThe paper presents original theoretical bounds of regret for the new scalable approximate Thompson sampling technique and the basic supporting empirical study.\n\n# Quality\nSee warning above on the theoretical side. The extended empirical study would be beneficial.\n\n# Clarity\nThe paper is well-written. \n\n# Significance\nThe results seem novel and useful for applications. As a byproduct of the theoretical study authors identify the regime of the scalable approximate Thompson sampling that is efficient and attains the regret asymptotic order of the exact analog. Probably the main problem is that the theoretical study does not align with the empirical one (uses a different setting).\n\n# Remarks\n- Lines 78–79. It is worth mentioning that the usual definition of sub-Gaussianity allows for nonzero mean and contains a multiplicative constant before the exponential function. Furthermore, it seems unnecessary and a bit confusing to introduce the sub-Gaussianity assumption for noise and then working with Gaussian noise. In my opinion, just assuming the Gaussianity of noise would be fairer and clearer.\n- Line 94. This line refers the reader to the appendix for details on RKHS. After a quick look over the appendix I could not find a discussion on RKHS. If it does actually exist, a more explicit reference would help.\n- Equation (2). $\\tau$ is not defined/explained in text.\n- Line 120. Probably “is given” -> “are given”.\n- Line 130. I suggest doing something like “we can decompose our chosen kernel”->“we can usually decompose our chosen kernel”. After all, we can decompose only sufficiently regular kernels.\n- Lines 170 and 178. It would be useful to elaborate on why the formulas are what they are, i.e. why do we have $\\alpha_t (u_t - m_t) + m_t$ instead of just $\\alpha_t u_t.$ Additionally, it would be useful to expand on why you consider two different sampling rules (given by (4) and (5)) and what are their benefits relative to each other.\n- Line 194. “assumption”->”assumptions”.\n- Equation after line 196. The notation for constant $\\varepsilon_t$ overloads the notation for noise.\n- Line 205. “$\\tilde{u}_t (\\tilde{\\sigma}_t(x) + \\varepsilon_t)$” -> “$\\tilde{u}_t(\\delta) (\\tilde{\\sigma}_t(x) + \\varepsilon_t)$” because otherwise it looks like $\\tilde{u}_t$ is evaluated at $\\tilde{\\sigma}_t(x) + \\varepsilon_t$ rather than multiplied by it.\n- Theorem 1. It would be useful to elaborate on the fact that the seemingly linear term in this regret will be subsequently removed by the very small value of $\\varepsilon$.\n- Line 222. Probably “Decomposed”->”Decoupled”.\n- Line 350. Probably “$m_t$” -> “$m$”.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Some limitations are addressed. The societal impact is adequately addressed.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}, {"review_id": "2tacer_RkOn", "reviewer": "Reviewer_52bF", "summary": "The manuscript provides a theoretical analysis of Thompson sampling with the sparse Gaussian processes (S-GP-TS) for Bayesian optimization (BO) as sparse GPs have been successfully used to improve the scalability of Thompson sampling with the exact GP (GP-TS). The manuscript shows that TS with any (sparse) approximation GP can achieve the same regret order as the (exact) GP-TS with certain conditions. It also includes numerical experiments which shows empirical performance (i.e., the regret bound over iterations) of S-GP-TS. ", "questions": "", "limitations": "", "rating": 6, "confidence": 1, "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 manuscript does not proposes new methodologies but I think it has certain contributions. The manuscript provides theoretical background and guarantees for sparse GPs for Thompson sampling. \n\nI think that only one real-world dataset is not enough to make readers convinced about the effectiveness of the proposed method in practice. I also think that the authors could include the computational complexity of each method (BNN-TS or S-GP-TS) per each evaluation. \n\nLines 136 and 137: the subscript t is omitted (i.e., \\phy_m(x) and \\Lambda_m)\n\nLine 306: “BO iterations” is used without its definition. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes. The authors addressed the limitations in Discussion section.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: I am ultimately recommending acceptance for this paper on the basis of the reviewer consensus that the theoretical contributions made in this paper are strong. Namely, this is the first paper I am aware of to (1) present a regret bound for Bayesian optimization with sparse GP models using any acquisition function at all, and (2) moreover, the theoretical results are specifically for a recent scalable approximation of Thompson sampling using sparse GPs, noting that Thompson sampling is one of those rare GP operations not made more computationally tractable by simply using sparse GPs alone.

However, I want to echo Reviewer eTmV's concerns about the choice of \alpha_{t}=1 in the experiments. Keep in mind, the actual acquisition function that you've demonstrated convergence for is *not* Thompson sampling as we would all understand it, or even an approximation thereof: it's a modification to Thompson sampling in which the (co-)variance envelope samples are drawn from shrinks over time. This is fine, and indeed there is obvious precedent for such a modification for theoretical purposes, because as you point out people who use UCB don't actually typically anneal the trade-off parameter. That said, Srinivas et al. absolutely did evaluate their mechanism with the annealed trade-off parameter: they did not evaluate a constant \beta_{t} but rather scaled it down by a factor of five, leaving a schedule for \beta_{t} that still annealed over time. I don't agree that setting \alpha_{t}=1 is at all analogous to their experimental evaluation.

While I completely understand that in practice, people will simply just use Thompson sampling (which has no parameter \alpha_{t}), there's frankly just no good justification for not at least presenting some results using an \alpha_{t} annealing schedule within a constant factor of the one suggested by your theory. You should include these results, if not in the main body of the paper then in the supplementary materials.
Decision: Accept (Poster)