Fast, Precise Thompson Sampling for Bayesian Optimization

Authors: david  sweet
Venue: NeurIPS_cc_2024_Workshop_BDU
Published: 2024-08-25 02:49:30

Abstract: Thompson sampling (TS) has optimal regret and excellent empirical performance
in multi-armed bandit problems. Yet, in Bayesian optimization, TS underperforms
popular acquisition functions (e.g., EI, UCB). TS samples arms according to the
probability that they are optimal. A recent algorithm, P-Star Sampler (PSS), per-
forms such a sampling via Hit-and-Run. We present an improved version, Stagger
Thompson Sampler (STS). STS more precisely locates the maximizer than does
TS using less computation time. We demonstrate that STS outperforms TS, PSS,
and other acquisition methods in numerical experiments of optimizations of sev-
eral test functions across a broad range of dimension. Additionally, since PSS
was originally presented not as a standalone acquisition method but as an input to
a batching algorithm called Minimal Terminal Variance (MTV), we also demon-
strate that STS matches PSS performance when used as the input to MTV.

Main paper: 
 introduction: 
 Bayesian optimization (BO) is applied to many experiment-and simulation-based optimization problems in science and engineering [17] . The aim of BO methods is to minimize the number of measurements needed to find a good system configuration. Measurements are taken in a sequence of batches of one or more arms -where an arm is one system configuration. System performance is measured for each arm in a batch, then a new batch is produced, performance is measured, and so on.Thompson sampling (TS) samples arms according to the probability that they maximize system performance [22] . Let's denote an arm as x a , its performance as f (x a ), and the probability that an arm is best, i.e., that x a = argmax x f (x), as p * (x). TS draws an arm as: x a ∼ p * (x). TS has optimal or near-optimal regret [1, 3] for multi-armed bandits (MAB) and is also used in BO [9] . Application of TS to BO is not straightforward. Since BO arms are typically continuous -e.g., x ∈ [0, 1] d -sampling from p * (x) is non-trivial. The usual approach is to first sample many candidate arms uniformly, x i ∼ U([0, 1] d ), then draw a value, y i , from a model distribution of f (x) at each x i . The x i that yields the highest-valued draw, i.e., x a = argmax xi y(x i ), is taken as the arm. BO methods usually model y(x) with a Gaussian process, GP [16] , i.e. y(x) ∼ N (µ(x), σ 2 (x)). While appealing for its simplicity, TS, implemented as described, tends to underperform popular acquisition functions such as Expected Improvement (EI) [10] or Upper Confidence bound (UCB) [5] (see Subsection 3.2). 
stagger thompson sampling: 
 The TS arm candidates, being uniform in [0, 1] d , are unlikely to fall where p * (x) has high density. If we imagine that the bulk of p * (x) lies in a hypercube of side ε < 1, then the probability that a randomly-chosen x ∈ [0, 1] d falls in the hypercube is just the hypercube volume, v = ε d . Note that (i) v decreases exponentially with dimension, and (ii) ε, thus v, decreases with each additional measurement added to the model of p * (x) (see figure 5 (c) in appendix A). Effect (ii) is the aim of BO -to localize the maximizer. Effect (i) is the "curse of dimensionality". Thus, the number of candidates required to find the bulk of p * (x) (i) increases exponentially in dimension, and (ii) increases with each additional measurement (albeit in a non-obvious way).The P-Star Sampler (PSS) [17] is a Hit-and-Run [19] sampler with a Metropolis filter [7] . Our algorithm, Stagger Thompson Sampler (STS), is, also, but differs in several details, discussed below algorithm 1.1: if no measurements yet then 2: return x i ∼ U([0, 1] d ) ▷ Take p * (x) prior as uniform in x 3: x * = argmax x µ(x) ▷ µ(x)x t = U([0, 1] d ) ▷ Perturbation target 7:s = e -kU ([0,1]) ▷ A "stagger" perturbation length 8:x′ a = x a + s(x t -x a ) ▷ Perturbations 9: [y, y ′ ] ∼ GP([x a , x ′ a ]) ▷ Joint sample 10: x a ← x ′ a if y ′ > y ▷ MH acceptance 11: return x a ▷ A sample from p * (x)Algorithm 1 modifies vanilla Hit-and-Run in two ways: (i) Instead of initializing randomly, we initialize an arm candidate, x a , at x * . (ii) Instead of perturbing uniformly along some directionsince the scale of p * (x) is unknown and may be small -we choose the length of the perturbation uniformly in its exponent, i.e. as ∼ e -kU , a log-uniform random variable. (k is a hyperparameter, which we choose to be k = ln 10 -6 ). Numerical ablation studies in appendix B show that these modifications improve performance in optimization. We refer to the log-uniform perturbation as a "stagger" proposal, following [21] . Besides being empirically effective, a stagger proposal also obviates the need to adapt the scale of the proposal distribution, as is done in PSS (which uses a Gaussian proposal). We see this as a valuable, practical simplification.Since a log-uniform distribution is a symmetric proposal we expect the Markov chain generated by algorithm 1 to converge to p * (x). Appendix A provides some numerical support for this.Perturbations are made along a line from x a to a target point, x t , which is chosen uniformly inside the bounding box, [0, 1] d . This ensures that the final perturbation [a convex combination of points in the (convex) bounding box] will lie inside the bounding box. It also simplifies the implementation somewhat, since boundary detection is unnecessary. PSS performs a bisection search to find the boundary of the box along a randomly-oriented line passing through x a .We accept a perturbation of x, x ′ , with Metropolis acceptance probability p acc = min{1, p * (x ′ )/p * (x)}. As a coarse (and fast) approximation to p acc , we follow PSS and just take a single joint sample from GP and accept whichever point, x or x ′ , has a larger sample value. Note that this is a Thompson sample from the set {x, x ′ }, so we might describe STS as iterated Thompson sampling.Appendix A offers numerical evidence that (i) samples from STS are nearer the true maximizer than are samples from TS, and (ii) STS produces samples more quickly than standard TS while better approximating p * (x).Previous work applying MCMC methods to Thompson sampling include random-walk Metropolis algorithms constrained to a trust region [25] and a sequential Monte Carlo algorithm [2] .  
numerical experiments: 
 To evaluate STS, we optimize various test functions, tracking the maximum measured function value at each round and comparing the values to those found by other methods. We use the term round to refer to the generation of one or more arms followed by the measurement of them.To introduce our comparison methodology, we compare STS to a few other optimization methods, in particular to TuRBO, a trust-region-enhanced Thompson sampling method [6] . In [6] , the authors optimize the Ackley function in 200 dimensions with 100 arms/round on a restricted subspace of parameters. Figure 2 optimizes the same function on the standard parameter space using STS, TuRBO, and other methods. STS finds higher values of y more quickly than the other methods. The error areas are twice the standard error over 10 runs. STS (sts) finds higher values more quickly than other optimization methods: turbo-1 -TuRBO [6] with one trust region. cma -CMA-ES [8] , an evolution strategy. random -Choose arms uniformly randomly (serving as a baseline). [6] .We can summarize each method's performance in Figure 2 with a single number, which we'll call the score. At each round, i, find the maximum measured values so far for each method, m: y i,m . Rank these values across m and scale: r i,m = [rank(y i,m ) -1]/(M -1), where M is the number of methods. Repeat this for every round, i, then average over all R rounds to get the score: s m = R i r i,m /R. The scores in figure 2 are s sts = 1, s turbo-1 = 2/3, s cma = 1/2, and s random = 0.Using an normalized score enables us to average over runs on different functions (which, in general, have different scales for y). Using a rank-based score prevents a dramatic result, like the one in figure 2 , from dominating the average.In our experiments below we optimize over nine common functions. To add variety to the function set and to avoid an artifact where an optimization method might coincidentally prefer to select points near a function's optimum (e.g., at the center of the parameter space), we randomly distort each function as in [17] , repeating the optimization 30 times with different random distortions.We compare STS to other BO methods in dimensions 3 through 300, all generating one arm per round. For each dimension, each method's score is averaged over all test functions. See Figure 3 . STS has the highest score in each dimension, and its advantage appears to increase with dimension. Data from dimensions 1, 10, and 100 (unpublished for space) follow the same pattern. We optimize for max(30, num dim) rounds with num arms / round over the functions ackley, dixonprice, griewank, levy, michalewicz, rastrigin, rosenbrock, sphere, and stybtang [20] with random distortions (see section 3.1). Error bars are two standard errors over all functions and 30 runs/function. Figure represents a total of 874,800 function evaluations. (We were not able to run pss for num dim=300 due to long computation times.The various optimization methods are: sts -Stagger Thompson Sampling (Algorithm 1), random -Uniformly-random arm, sobol -Uniform, space-filling arms [18, Chapter 5], sr -Simple Regret, µ(x), ts -Thompson Sampling [9] , ucb -Upper Confidence Bound [5, 24] , ei -Expected Improvement [12, 24] , gibbon -GIBBON, an entropy-based method [13] , and optuna -an open-source optimizer [15] employing a tree-structured Parzen estimator [23] . For the methods ei, lei, ucb, sr, gibbon, and turbo-1, we initialize by taking a Sobol' sample for the first round's arm. sts and ts do not require initialization.STS makes no explicit accommodations for higher-dimensional problems yet performs well on them. Of the methods evaluated, only turbo-1 specifically targets higher-dimensional problems [6] , so it may be valuable for future work to compare STS to other methods specifically designed for such problems. (See references to methods in [6] .)Thompson sampling can be extended to batches of more than one arm simply by taking multiple samples from p * (x), e.g., by running Algorithm 1 multiple times. However, this approach can be inefficient [14] because some samples -since they are independently generated -may lie very near each other and, thus, provide less information about f (x) than if they were to lie farther apart. This problem, that of generating effective batches of arms, is not unique to TS but exists for all approaches to acquisition, and there are various methods for dealing with it [14] [24] [13] .One method, Minimal Terminal Variance (MTV) [17] , minimizes the post-measurement, average variance of the GP, weighted by p * (x):EQUATIONapproximated by i σ 2 (x i |x a ), where x i are drawn from x i ∼ p * (x) with P-Star Sampler (PSS).MTV is interesting, in part, because it can design experiments both when prior measurements are available and ab initio (e.g., at initialization time). It not only outperforms acquisition functions (like EI or UCB) but the same formulation also outperforms common initialization methods, such as Sobol' sampling [17] .We modify MTV to draw x i ∼ p * (x) using STS instead of PSS. Note, also, that the arms, x a , that minimize M T V (x a ) are not drawn from the set x i but are chosen by a continuous minimization algorithm (specifically, scipy.minimize, as implemented in [11] ), such that x a = argmin x ′ a i σ 2 (x i |x ′ a ). Figure 4 compares MTV, with P-Star Sampler replaced by STS (mtv+sts), to other methods using various dimensions (num dim) and batch sizes (num arms): mtv -MTV, as in [17] , lei -q-Log EI, an improved EI [4] , and dpp -DPP-TS [14] , a diversified-batching TS. For the methods ei, ucb, sr, gibbon, and dpp, we initialize by taking Sobol' samples for the first round. sts, mtv, and mtv+sts do not require initialization.Figure 4 roughly reproduces figure 3 of [17] , adding more methods and extending to higher dimensions. Additionally, we include pure PSS and STS sampling, where arms are simply independent draws from p * (x), to highlight the positive impact of MTV on batch design. When MTV's input samples come from STS (mtv+sts), performance is similar to the original MTV (mtv). 
conclusion: 
 We presented Stagger Thompson Sampler, which is novel in several ways:• It outperforms not only the standard approach to Thompson sampling but, also, popular acquisition functions, a trust region Thompson sampling method, and an evolution strategy. • It is simpler and more effective at acquisition than an earlier sampler (P-Star Sampler).• It works on high-dimensional problems without modification.Additionally, the combination MTV+STS is unique in that it applies to so broad a range of optimization problems: It solves problems with zero or more pre-existing measurements, with one or more arms/batch, and in dimensions ranging from low to high.  
 
 Appendix:  



Reviews: review: title: title: Review of Submission 6
main_review: main_review: **Summary:** This work presents a novel approach for fast and exact Thompson sampling in Gaussian process-based Bayesian optimization. The authors introduce an MCMC scheme that enables efficient sampling of points based on their probability of being the global optimum. The proposed method shows significant performance gains in various numerical experiments.

**Strengths:** The method is straightforward to implement, and the empirical results are impressive, demonstrating notable improvements over existing approaches.

**Weaknesses:** My primary concern is the lack of clarity regarding the "exactness" of the proposed sampling scheme. It is not immediately clear that the introduced MCMC approach samples from the posterior distribution over the optimum. A detailed discussion addressing this aspect is necessary.
title: title: Review of the submission
main_review: main_review: The article presents an MCMC-inspired method for sampling from the Thompson sampling distribution in Gaussian process bandits, which is also a credible (if slightly misguided) approach to Bayesian optimisation - in pure Bayesian optimisation we should not care about incurring poor performance on the way through, whereas bandits and TS do need to care about just how bad some arms might be. Nevertheless, I am reviewing purely whether this is a good article on how to perform Thompson sampling in GPs.

The approach builds on the p* sampler of Ren and Sweet. What is proposed is not daft, and seems to work well in the presented experiments. However I have significant misgivings about the clarity of the paper, and think it is not ready for presentation.

My three key concerns are:
1. The algorithm is presented in a reproducible way, but there is no analysis of what the stationary distribution of this Markov chain is, and therefore no justification for why it might be sampling from approximately the Thompson sampling distribution.
2. The interfacing of MTV with the approach in section 3.3 is not at all well explained. Is it the case that too many Thompson samples are produced then MTV thins these down?
3. It is claimed that the new method is significantly computationally faster than the p* sampler on which it builds. But both are simply MCMC schemes with an arbitrary run length. I cannot understand how a reasonable comparison of computational effort can be made, and certainly no attempt is made to do so.

More minor comments, to assist a future revision:
- Use of "ex. " instead of "for example" or "eg" is non-standard and quite confusing
- In the abstract, the use of "arm dimension" is a little jarring for some reason. Perhaps "action space dimension"?
- In l32, v is the volume of a fixed sized box. It doesn't change with data. So this statement is incorrect.
- In Alg 1, what is p_max?
- It's unfortunate that the key prior work cannot be implemented in the first example (l59). It would be much better to start with an example where the prior work can be implemented then to move on.

The topic is kind of interesting for the workshop, and it's possible that a poster would be interesting. But I suspect there will be better submissions.