Fast, Precise Thompson Sampling for Bayesian Optimization

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}=\operatorname*{argmax}_{x}f(x)$, as $p_{*}(x)$. TS draws an arm as: $x_{a}\sim 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\in[0,1]^{d}$ – sampling from $p_{*}(x)$ is non-trivial. The usual approach is to first sample many candidate arms uniformly, $x_{i}\sim\mathcal{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}=\operatorname*{argmax}_{x_{i}}y(x_{i})$, is taken as the arm. BO methods usually model $y(x)$ with a Gaussian process, $\mathcal{GP}$ [16], i.e. $y(x)\sim\mathcal{N}(\mu(x),\sigma^{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).

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 $\varepsilon<1$, then the probability that a randomly-chosen $x\in[0,1]^{d}$ falls in the hypercube is just the hypercube volume, $v=\varepsilon^{d}$. Note that (i) $v$ decreases exponentially with dimension, and (ii) $\varepsilon$, 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.

Algorithm 1 modifies vanilla Hit-and-Run in two ways: (i) Instead of initializing randomly, we initialize an arm candidate, $x_{a}$, at $\tilde{x}_{*}$. (ii) Instead of perturbing uniformly along some direction – since 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 $\sim 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^{\prime}$, with Metropolis acceptance probability $p_{acc}=\min\{1,p_{*}(x^{\prime})/p_{*}(x)\}$. As a coarse (and fast) approximation to $p_{acc}$, we follow PSS and just take a single joint sample from $\mathcal{GP}$ and accept whichever point, $x$ or $x^{\prime}$, has a larger sample value. Note that this is a Thompson sample from the set $\{x,x^{\prime}\}$, 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].

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.

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.

This work was carried out in affiliation with Yeshiva University. The author is additionally affiliated with DRW Holdings, LLC. This work was supported, in part, by a grant from Modal.