Optimizing Posterior Samples for Bayesian Optimization via Rootfinding

## Abstract

Abstract Bayesian optimization devolves the global optimization of a costly objective function
to the global optimization of a sequence of acquisition functions.
This inner-loop optimization can be catastrophically difficult if it involves posterior samples,
especially in higher dimensions.
We introduce an efficient global optimization strategy for posterior samples
based on global rootfinding.
It provides gradient-based optimizers with
judiciously selected starting points, designed to combine exploitation and exploration.
The algorithm scales practically linearly to high dimensions.
For posterior sample-based acquisition functions such as
Gaussian process Thompson sampling (GP-TS) and variants of entropy search,
we demonstrate remarkable improvement in both inner- and outer-loop optimization,
surprisingly outperforming alternatives like EI and GP-UCB in most cases.
We also propose a sample-average formulation of GP-TS,
which has a parameter to explicitly control exploitation
and can be computed at the cost of one posterior sample.
Our implementation is available at https://github.com/UQUH/TSRoots .

## 1 Introduction

Bayesian optimization (BO) is a highly successful approach to
the global optimization of expensive-to-evaluate black-box functions,
with applications ranging from hyper-parameter training of machine learning models
to scientific discovery and engineering design [, , , ].
Many BO strategies are also backed by strong theoretical guarantees
on their convergence to the global optimum [, , , ].

Consider the global optimization problem $\min_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x})$
where $\mathbf{x}\in\mathcal{X}\subset\mathbb{R}^{d}$ represents the vector of input variables and $f(\mathbf{x})\in\mathbb{R}$ the objective function which can be evaluated at a significant cost,
subject to observation noise.
At its core, BO is a sequential optimization algorithm
that uses a probabilistic model of the objective function to guide its evaluation decisions.
Starting with a prior probabilistic model and some initial data,
BO derives an acquisition function $\alpha(\mathbf{x})$ from the posterior model,
which is much easier to evaluate than the objective function
and often has easy-to-evaluate derivatives.
The acquisition function is then optimized globally, using off-the-shelf optimizers,
to provide a location to evaluate the objective function.
This process is iterated until some predefined stopping criteria are met.

Effectively there are two nested iterations in BO:
the outer-loop seeks to optimize the objective function $f(\mathbf{x})$,
and the inner-loop seeks to optimize the acquisition function $\alpha(\mathbf{x})$ at each BO iteration.
The premise of BO is that the inner-loop optimization can be solved accurately and efficiently,
so that the outer-loop optimization proceeds informatively with a negligible added cost.
In fact, the convergence guarantees of many BO strategies
assume exact global optimization of the acquisition function.
However, the efficient and accurate global optimization of acquisition functions
is less trivial than it is often assumed to be [].

Acquisition functions are, in general, highly non-convex and have many local optima.
In addition, many common acquisition functions
are mostly flat surfaces with a few peaks [], which take up an overwhelmingly large portion of the domain as the input dimension grows.
This creates a significant challenge for generic global optimization methods.

Some acquisition functions involve sample functions from the posterior model,
which need to be optimized globally.
Gaussian process Thompson sampling (GP-TS) [] uses posterior samples directly as random acquisition functions.
In many information-theoretic acquisition functions
such as predictive entropy search (PES) [], max-value entropy search (MES) [], and joint entropy search (JES) [, ],
multiple posterior samples are drawn and optimized to find their global optimum location and/or value.
Such acquisition functions are celebrated for their nice properties in BO:
TS has strong theoretical guarantees [, ]
and can be scaled to high dimensions [];
information-theoretic acquisition functions are grounded
in principles for optimal experimental design [];
and both types can be easily parallelized in synchronous batches
[, , ].
However, posterior samples are much more complex than other acquisition functions as they fluctuate throughout the design space, and are less smooth than the posterior mean and marginal variance.
The latter are the basis of many acquisition functions,
such as expected improvement (EI) [], probability of improvement (PI) [],
and upper confidence bound (GP-UCB) [].
As a consequence, posterior samples have many more local optima,
and the number scales exponentially with the input dimension.

While there is a rich literature on prior probabilistic models and acquisition functions for BO,
global optimization algorithms for acquisition functions have received little attention.
One class of global optimization methods is derivative-free,
such as the dividing rectangles (DIRECT) algorithm [],
covariance matrix adaptation evolution strategy (CMA-ES) algorithms [],
and genetic algorithms [].
Gradient-based multistart optimization, on the other hand,
is often seen as the best practice to reduce the risk of getting trapped in local minima [],
and enjoys the efficiency of being embarrassingly parallelizable.
For posterior samples, their global optimization may use
joint sampling on a finite set of points [],
or approximate sampling of function realizations followed by gradient-based optimization
[, ].
The selection of starting points is crucial for the success of gradient-based multistart optimization.
This selection can be deterministic (e.g., grid search), random [, ], or adaptive [].

In this paper, we propose an adaptive strategy for selecting starting points for
gradient-based multistart optimization of posterior samples.
This algorithm builds on the decomposition of posterior samples by pathwise conditioning,
taps into robust software in univariate function computation based on univariate global rootfinding,
and exploits the separability of multivariate GP priors.
Our key contributions include:

- •
A novel strategy for the global optimization of posterior sample-based acquisition functions.
The starting points are dependent on the posterior sample,
so that each is close to a local optimum that is a candidate for the global optimum.
The selection algorithm scales linearly to high dimensions.
- •
We give empirical results across a diverse set of problems
with input dimensions ranging from 2 to 16,
establishing the effectiveness of our optimization strategy.
Although our algorithm is proposed for the inner-loop optimization of posterior samples,
perhaps surprisingly, we see significant improvement in outer-loop optimization performance,
which often allows acquisition functions based on posterior samples to converge faster
than other common acquisition functions.
- •
A new acquisition function via the posterior sample average that explicitly controls the exploration–exploitation balance [],
and can be generated at the same cost as one posterior sample.

## 2 General Background

#### Gaussian Processes.

Consider an unknown function $f_{\text{true}}:\mathcal{X}\mapsto\mathbb{R}$, where domain $\mathcal{X}\subset\mathbb{R}^{d}$.
We can collect noisy observations of the function through the model $y^{i}=f_{\text{true}}(\mathbf{x}^{i})+\varepsilon^{i}$,
$i\in\{1,\cdots,n\}$, with $\boldsymbol{\varepsilon}\sim\mathcal{N}_{n}(0,\boldsymbol{\Sigma})$.
To model the function $f_{\text{true}}$, we use a Gaussian process (GP) as the prior probabilistic model: $f\sim\pi\in\mathcal{GP}$.
A GP is a random function $f$ such that for any finite set of points $X=\{\mathbf{x}^{i}\}_{i=1}^{n}$, $n\in\mathbb{N}$,
the values $\mathbf{f}_{n}=(f(\mathbf{x}^{i}))_{i=1}^{n}$ have a multivariate Gaussian distribution $\mathcal{N}_{n}(\boldsymbol{\mu}_{n},\mathbf{K}_{n,n})$,
with mean $\boldsymbol{\mu}_{n}=(\mu(\mathbf{x}^{i}))_{i=1}^{n}$
and covariance $\mathbf{K}_{n,n}=[\kappa(\mathbf{x}^{i},\mathbf{x}^{j})]_{i\in n}^{j\in n}$.
Here, $\mu(\mathbf{x})$ is the mean function and $\kappa(\mathbf{x},\mathbf{x}^{\prime})$ is the covariance function.

#### Decoupled Representation of GP Posteriors.

Given a data set $\mathcal{D}=\{(\mathbf{x}^{i},y^{i})\}_{i=1}^{n}$,
the posterior model $f|\mathcal{D}$ is also a GP.
Samples from the posterior have a decoupled representation called pathwise conditioning, originally proposed in []:

$$ $(f|\mathcal{D})(\mathbf{x})\overset{\text{d}}{=}f(\mathbf{x})+\boldsymbol{ \kappa}_{\cdot,n}(\mathbf{x})(\mathbf{K}_{n,n}+\boldsymbol{\Sigma})^{-1}( \mathbf{y}-\mathbf{f}_{n}-\boldsymbol{\varepsilon}),\quad f\sim\pi,\; \boldsymbol{\varepsilon}\sim\mathcal{N}_{n}(0,\boldsymbol{\Sigma}),$ (1) $$

where $\overset{\text{d}}{=}$ denotes an equality in distribution,
and $\boldsymbol{\kappa}_{\cdot,n}(\mathbf{x})=(\kappa(\mathbf{x},\mathbf{x}^{i}))_
{i=1}^{n}$ is the canonical basis.

## 3 Global Optimization of Posterior Samples

In this section, we introduce an efficient algorithm for the global optimization of posterior sample-based acquisition functions.
For this, we exploit the separability of prior samples and useful properties of posterior samples to judiciously select a set of starting points
for gradient-based multistart optimizers.

#### Assumptions.

Following Section 2, we impose a few common assumptions throughout this paper:
(1) the domain is a hypercube: $\mathcal{X}=\prod_{i=1}^{d}[\underline{x}_{i},\overline{x}_{i}]$;
(2) prior covariance is separable: $\kappa(\mathbf{x},\mathbf{x}^{\prime})=\prod_{i=1}^{d}\kappa_{i}(x_{i},x_{i}^{
\prime})$;
(3) prior samples are continuously differentiable: $f(\mathbf{x};\omega)\in C^{1}$.
Without loss of generality, we also assume that the prior mean $\mu(\mathbf{x})=0$:
any non-zero mean function can be subtracted from the data by replacing $f_{\text{true}}$ with $f_{\text{true}}-\mu$.
While additive and multiplicative compositions of univariate kernels can be used in the prior [],
assumption (2) is the most popular choice in BO.
It is possible to extend our method to generalized additive models.

### 3.1 TS-roots Algorithm

To find the global minimum $(\widetilde{\mathbf{x}}_{\widetilde{\omega}}^{\star},\widetilde{f}_{\widetilde
{\omega}}^{\star})$ of a posterior sample $\widetilde{f}_{\widetilde{\omega}}(\mathbf{x})$,
we observe that given the assumptions, Eq. 1 can be rewritten as:

$$ $\widetilde{f}(\mathbf{x};\widetilde{\omega})=f(\mathbf{x};\omega)+b(\mathbf{x} ;\widetilde{\omega}),\quad f(\mathbf{x};\omega)=\prod_{i=1}^{d}f_{i}(x_{i}; \omega_{i}),\quad b(\mathbf{x};\widetilde{\omega})=\sum_{j=1}^{n}v_{j}\kappa( \mathbf{x},\mathbf{x}^{j}).$ (2) $$

Here, the prior sample $f(\mathbf{x})$ is a separable function determined by the random bits $\omega$.
Data adjustment $b(\mathbf{x})$ is a sum of the canonical basis with coefficients
$\mathbf{v}=(\mathbf{K}_{n,n}+\boldsymbol{\Sigma})^{-1}(\mathbf{y}-\mathbf{f}_{
n}-\boldsymbol{\varepsilon})$.
Both the data adjustment and the posterior sample are determined by the random bits
$\widetilde{\omega}=(\omega,\boldsymbol{\varepsilon})$.
In the following, we denote the prior and the posterior samples as $f_{\omega}$ and $\widetilde{f}_{\widetilde{\omega}}$, respectively.

The global minimization of a generic function, in principle, requires finding all its local minima and then selecting the best among them.
However, computationally efficient approaches to this problem are lacking even in low dimensions
and, more pessimistically, the number of local minima grows exponentially as the domain dimension increases.
The core idea of this work is to use the prior sample $f_{\omega}$ as a surrogate of the posterior sample $\widetilde{f}_{\widetilde{\omega}}$ for global optimization.
Another key is to exploit the separability of the prior sample for
efficient representation and ordering of its local extrema.

We define TS-roots as a global optimization algorithm for GP posterior samples
(given the assumptions) via gradient-based multistart optimization,
where the starting points include
a subset of the local minima $\breve{X}_{\omega}$ of a corresponding prior sample, $S_{\text{e}}\subseteq\breve{X}_{\omega}$,
and a subset of the observed locations, $S_{\text{x}}\subseteq X$.
We call $S_{\text{e}}$ the exploration set and $S_{\text{x}}$ the exploitation set.
Let

$$ $S_{\text{o}}=\operatorname*{arg\,mink}_{\mathbf{x}\in\breve{X}_{\omega}}(f_{ \omega}(\mathbf{x}),n_{\text{o}})$ (3) $$

be the $n_{\text{o}}$ smallest local minima of the prior sample
and $\widetilde{\mu}(\mathbf{x})=\boldsymbol{\kappa}_{\cdot,n}(\mathbf{x})(\mathbf{
K}_{n,n}+\boldsymbol{\Sigma})^{-1}\mathbf{y}$
be the posterior mean function.
The starting points are defined as:

$$ $S=S_{\text{e}}\cup S_{\text{x}},\quad S_{\text{e}}=\operatorname*{arg\,mink}_{ \mathbf{x}\in S_{\text{o}}}(\widetilde{f}_{\widetilde{\omega}}(\mathbf{x}),n_{ \text{e}}),\quad S_{\text{x}}=\operatorname*{arg\,mink}_{\mathbf{x}\in X}( \widetilde{\mu}(\mathbf{x}),n_{\text{x}}).$ (4) $$

Here, $n_{\text{e}}$ and $n_{\text{x}}$ are imposed to control the cost of gradient-based multistart optimization,
and $n_{\text{o}}$ is set to limit the number of evaluations of $\widetilde{f}_{\widetilde{\omega}}$ in the determination of $S_{\text{e}}$.
Considering the cost difference, in general we can have $n_{\text{o}}\gg n_{\text{e}}$.

Figure: Figure 1: Illustrations of exploration and exploitation sets for the global optimization of GP-TS acquisition functions in one dimension (top row) and two dimensions (bottom row). Left column: When the global minimum $\widetilde{\mathbf{x}}_{\widetilde{\omega}}^{\star}$ of the GP-TS acquisition function lies outside the interpolation region, it is typically identified by starting the gradient-based multistart optimizer at a local minimum of the prior sample. Right column: When $\widetilde{\mathbf{x}}_{\widetilde{\omega}}^{\star}$ is within the interpolation region, it can be found by starting the optimizer at either an observed location or a local minimum of the prior sample.
Refer to caption: x1.png

### 3.2 Relations between the Local Minima of Prior and Posterior Samples

Figure 1 shows several posterior samples $\widetilde{f}_{\widetilde{\omega}}$ in one and two dimensions,
each marked with its local minima $\breve{\widetilde{X}}_{\widetilde{\omega}}$ and global minimum $\widetilde{\mathbf{x}}_{\widetilde{\omega}}^{\star}$.
Here the exploration set $S_{\text{e}}$ is the local minima $\breve{X}_{\omega}$ of $f_{\omega}$,
and the exploitation set $S_{\text{x}}$ is the observed locations $X$.
We make the following observations:

- (1)
The prior sample $f_{\omega}$ is more complex than the data adjustment $b$
in the sense that it is less smooth and has more critical points.
The comparison of smoothness can be made rigorous in various ways: for example,
for GPs with a Matérn covariance function where the smoothness parameter is finite,
$f_{\omega}$ is almost everywhere one time less differentiable than $b$
(see e.g., Garnett [] Sec. 10.2, Kanagawa et al. []).
- (2)
Item (1) implies that
when the prior sample $f_{\omega}$ is added to the smoother landscape of $b$,
each local minimum $\breve{\mathbf{x}}_{\omega}$ of $f_{\omega}$ will be located near a local minimum $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$ of
the posterior sample $\widetilde{f}_{\widetilde{\omega}}$.
Away from the observed locations $X$, each $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$
is closely associated with an $\breve{\mathbf{x}}_{\omega}$, with minimal change in location.
In the vicinity of $X$, an $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$ may have
both a data point $\mathbf{x}^{i}$ and an $\breve{\mathbf{x}}_{\omega}$ nearby,
but because of the smoothness difference of $f_{\omega}$ and $b$,
in most cases the one closest to $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$ is an $\breve{\mathbf{x}}_{\omega}$.
- (3)
It is possible that near $X$, sharp changes in $f_{\omega}$
may require sharp adjustments to the data,
which may move some $\breve{\mathbf{x}}_{\omega}$ by a significant distance,
or create new $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$
that are not near any $\breve{\mathbf{x}}_{\omega}$ or any $\mathbf{x}^{i}$.
- (4)
Searching from $\mathbf{x}^{i}$ with good observed values can discover
good $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$ in the vicinity of $X$,
which can pick up some local minima not readily discovered by $\breve{X}_{\omega}$.
This is especially true if $f_{\omega}$ is relatively flat near $\mathbf{x}^{i}$.
- (5)
Since the posterior sample $\widetilde{f}_{\widetilde{\omega}}$ adapts to the data set, searches from $\mathbf{x}^{i}$ will tend to converge to
a good $\breve{\widetilde{\mathbf{x}}}_{\widetilde{\omega}}$ among all the local optima near $\mathbf{x}^{i}$.
Even if the searches from $X$ cannot discover all the local minima in its vicinity,
they tend to discover a good subset of them.
Therefore, (4) can help address the issue in (3), if not fully eliminating it.
By combining subsets of $\breve{X}_{\omega}$ and $X$,
we can expect that the set of local minima of $\widetilde{f}_{\widetilde{\omega}}$ discovered with these starting points
includes the global minimum $\widetilde{\mathbf{x}}_{\widetilde{\omega}}^{\star}$ with a high probability with respect to $\widetilde{\omega}$.

### 3.3 Representation of Prior Sample Local Minima

For each component function $f_{i}(x_{i};\omega_{i})$ of the prior sample $f_{\omega}(\mathbf{x})$,
define $h_{i}(\underline{x}_{i})=f_{i}^{\prime}(\underline{x}_{i})$, $h_{i}(\overline{x}_{i})=-f_{i}^{\prime}(\overline{x}_{i})$, and $h_{i}(x_{i})=f_{i}^{\prime\prime}(x_{i})$ for $x_{i}\in(\underline{x}_{i},\overline{x}_{i})$.
We call a coordinate $\xi_{i}\in[\underline{x}_{i},\overline{x}_{i}]$ of mono type if $f_{i}(\xi_{i})h_{i}(\xi_{i})>0$
and call it of mixed type if $f_{i}(\xi_{i})h_{i}(\xi_{i})<0$.
Let $\mathring{\Xi}_{i}=\{\xi_{i,j}\}_{j=1}^{r_{i}}$ be the set of interior critical points of $f_{i}$
such that $\xi_{i,j}\in(\underline{x}_{i},\overline{x}_{i})$ and $f_{i}^{\prime}(\xi_{i,j})=0$, $j\in\{1,\cdots,r_{i}\}$.
Denote $\xi_{i,0}=\underline{x}_{i}$ and $\xi_{i,r_{i}+1}=\overline{x}_{i}$.
Partition the set of candidate coordinates $\Xi_{i}=\{\xi_{i,j}\}_{j=0}^{r_{i}+1}$
into mono type $\Xi_{i}^{(0)}$ and mixed type $\Xi_{i}^{(1)}$.
Proposition 1 gives a representation of the sets of strong local extrema of the prior sample.
Its proof and the set sizes therein are given in Appendix A.

###### Proposition 1 .

The set of strong local minima of the prior sample $f_{\omega}(\mathbf{x})$ can be written as:

$$ $\breve{X}_{\omega}=\breve{X}_{\omega}^{-}\sqcup\breve{X}_{\omega}^{+},\quad \breve{X}_{\omega}^{-}=\{\boldsymbol{\xi}\in\Xi^{(1)}:f_{\omega}(\boldsymbol{ \xi})<0\},\quad\breve{X}_{\omega}^{+}=\{\boldsymbol{\xi}\in\Xi^{(0)}:f_{\omega }(\boldsymbol{\xi})>0\},$ (5) $$

where tensor grids $\Xi^{(j)}=\prod_{i=1}^{d}\Xi_{i}^{(j)}$, $j\in\{0,1\}$.
The set $\widehat{X}_{\omega}$ of strong local maxima of $f_{\omega}(\mathbf{x})$ has an analogous representation,
and satisfies $\widehat{X}_{\omega}\sqcup\breve{X}_{\omega}=\Xi^{(0)}\sqcup\Xi^{(1)}$, where $\sqcup$ is the disjoint union.

###### Proposition 1 .

#### Critical Points of Univariate Functions via Global Rootfinding.

To compute the set $\mathring{\Xi}_{i}$ of all critical points of $f_{i}$
is to compute all the roots of the derivative $f^{\prime}_{i}$ on the interval $[\underline{x}_{i},\overline{x}_{i}]$.
Since $f^{\prime}_{i}$ is continuous, this can be solved robustly and efficiently
by approximating the function with a Chebyshev or Legendre polynomial
and solving a structured eigenvalue problem (see e.g., Trefethen []).
The roots algorithm for global rootfinding based on polynomial approximation
is given as Algorithm 2 in Appendix C.

### 3.4 Ordering of Prior Sample Local Minima

While the size of $\breve{X}_{\omega}$ grows exponentially in domain dimension $d$,
its representation in Eq. 5 enables an efficient algorithm
to compute the best subset $S_{\text{o}}$ (eq. 3) without enumerating its elements.
With Eq. 5,
we see that $\breve{X}_{\omega}^{-}$ consists of all the local minima of $f_{\omega}$ with negative function values.
Consider the case where $\breve{X}_{\omega}^{-}$ has at least $n_{\text{o}}$ elements
so that in the definition of $S_{\text{o}}$ we can replace $\breve{X}_{\omega}$ with $\breve{X}_{\omega}^{-}$,
which in turn can be replaced with $\Xi^{(1)}$.
As we will show later, the problem of finding the largest elements of $|f_{\omega}(\mathbf{x})|$ within $\Xi^{(1)}$
is easier than finding the smallest negative elements of $f_{\omega}(\mathbf{x})$.
Once the former is solved, we can solve the latter simply by removing the positive elements.
Therefore, we convert the problem of Eq. 3 into three steps:

- 1.
$S^{(1)}=\operatorname*{arg\,maxk}_{\mathbf{x}\in\Xi^{(1)}}(|f_{\omega}(\mathbf
{x})|,\alpha n_{\text{o}})$,
with buffer coefficient $\alpha\geq 1$;
- 2.
$\breve{S}^{-}=\{\mathbf{x}\in S^{(1)}:f_{\omega}(\mathbf{x})<0\}$,
so that $\breve{S}^{-}\subseteq\breve{X}_{\omega}^{-}$;
- 3.
$S_{\text{o}}=\operatorname*{arg\,mink}_{\mathbf{x}\in\breve{S}^{-}}(f_{\omega}
(\mathbf{x}),n_{\text{o}})$,
assuming that $|\breve{S}^{-}|\geq n_{\text{o}}$.

The last two steps are by enumeration and straightforward.
The first step can be solved efficiently using min-heaps,
with a time complexity that scales linearly in $\sum_{i=1}^{d}|\Xi_{i}^{(1)}|$ rather than $\prod_{i=1}^{d}|\Xi_{i}^{(1)}|$,
see Appendix B.
The coefficient $\alpha$ is chosen so that $|\breve{S}^{-}|\geq n_{\text{o}}$.
The case when $|\breve{X}_{\omega}^{-}|<n_{\text{o}}<|\breve{X}_{\omega}|$ can be handled similarly.
If $n_{\text{o}}\geq|\breve{X}_{\omega}|$, no subsetting is needed.
The overall procedure to compute $S_{\text{o}}$ is given in Algorithm 3 in Appendix C.

## 4 Related Works

#### Sampling from Gaussian Process Posteriors.

A prevalent method to sample GP posteriors with stationary covariance functions is via weight-space approximations based on Bayesian linear models of random Fourier features [].
This method, unfortunately, is subject to the variance starvation problem [, ] which can be mitigated using more accurate feature representations (see e.g., [, ]).
An alternative is pathwise conditioning [] that draws GP posterior samples by updating the corresponding prior samples.
The decoupled representation of the pathwise conditioning can be further reformulated as two stochastic optimization problems for the posterior mean and an uncertainty reduction term, which are then efficiently solved using stochastic gradient descent [].

#### Optimization of Acquisition Functions.

While their global optima guarantee the Bayes’ decision rule, BO acquisition functions are highly non-convex and difficult to optimize [].
Nevertheless, less attention has been given to the development of robust algorithms for optimizing these acquisition functions.
For this inner-loop optimization, gradient-based optimizers are often selected because of their fast convergence and robust performance [].
The implementation of such optimizers is facilitated by Monte Carlo (MC) acquisition functions whose derivatives are easy to evaluate [].
Gradient-based optimizers also allow multistart settings that use a set of starting points which can be, for example,
midpoints of data points [], uniformly distributed samples over the input variable space [, ],
or random points from a Latin hypercube design [].
However, multistart-based methods with random search may have difficulty determining the non-flat regions of acquisition functions, especially in high dimensions [].

#### Posterior Sample-Based Acquisition Functions.

As discussed in Section 1, the family of posterior sample-based acquisition functions is determined from samples of the posterior.
GP-TS [] is a notable member that extends the classical TS for finite-armed bandits to continuous settings of BO
(see algorithms in Appendix D).
GP-TS prefers exploration by the mechanism that iteratively samples a function from the GP posterior of the objective function, optimizes this function, and selects the resulting solution as the next candidate for objective evaluation.
To further improve the exploitation of GP-TS, the sample mean of MC acquisition functions can be defined from multiple samples of the posterior [, ].
Different types of MC acquisition functions can also be used to inject beliefs about functions into the prior [].

## 5 Results

We assess the performance of TS-roots in optimizing benchmark functions. We then compare the quality of solutions to the inner-loop optimization of GP-TS acquisition functions obtained from our proposed method, a gradient-based multistart optimizer with uniformly random starting points, and a genetic algorithm.
We also show how TS-roots can improve the performance of MES.
Finally, we propose a new sample-average posterior function and show how it affects the performance of GP-TS.
The experimental details for the presented results are in Appendix F.

#### Optimizing Benchmark Functions.

We test the empirical performance of TS-roots on challenging minimization problems of five benchmark functions: the 2D Schwefel, 4D Rosenbrock, 10D Levy, 16D Ackley, and 16D Powell functions [].
The analytical expressions for these functions and their global minimum are given in Appendix E.

In each optimization iteration, we record the best observed value of the error $\log(y_{\min}-f^{\star})$
and the distance $\log\left(\left\|\mathbf{x}_{\min}-\mathbf{x}^{\star}\right\|\right)$,
where $y_{\min}$, $\mathbf{x}_{\min}$, $f^{\star}$, and $\mathbf{x}^{\star}$
are the best observation of the objective function in each iteration, the corresponding location of the observation,
the true minimum value of the objective function, and the true minimum location, respectively.
We compare the optimization results obtained from TS-roots and other BO methods, including GP-TS using decoupling sampling with random Fourier features (TS-DSRF), GP-TS with random Fourier features (TS-RF), expected improvement (EI) [], and lower confidence bound (LCB)—the version of GP-UCB [] for minimization problems.

Figure 2 shows the medians and interquartile ranges of solution values obtained from 20 runs of each of the considered BO methods.
The corresponding histories of solution locations are in Figure 7 of Appendix G.
With a fair comparison of outer-loop results (detailed in Appendix F),
TS-roots surprisingly shows the best performance on the 2D Schwefel, 16D Ackley, and 16D Powell functions,
and gives competitive results in the 4D Rosenbrock and 10D Levy problems.
Notably, TS-roots recommends better solutions than its counterparts, TS-DSRF and TS-RF, in high-dimensional problems
and offers competitive performance in low-dimensional problems.
Across all the examples, EI and LCB tend to perform better in the initial stages, while TS-roots shows fast improvement in later stages. This is because GP-TS favors exploration, which delays rewards.
The exploration phase, in general, takes longer for higher-dimensional problems.

Figure: Figure 2: Outer-loop optimization results for the (a) 2D Schwefel, (b) 4D Rosenbrock, (c) 10D Levy, (d) 16D Ackley, and (e) 16D Powell functions. The plots are histories of medians and interquartile ranges of solution values from 20 runs of TS-roots, TS-DSRF (i.e., TS using decoupled sampling with random Fourier features), TS-RF (i.e., TS using random Fourier features), EI, and LCB.
Refer to caption: x2.png

Figure: Figure 3: Inner-loop optimization results by rootfinding, a gradient-based multistart optimizer with random starting points (random multistart), and a genetic algorithm for (a) the 2D Schwefel and (b) 4D Rosenbrock functions. The plots are cumulative values of optimized GP-TS acquisition functions $\alpha_{k}^{\star}$, cumulative distances between new solution points ${\bf x}_{k}^{\star}$ and the true global minima ${\bf x}_{k}^{\text{t}}$ of the acquisition functions, and cumulative CPU times $t_{k}$ for optimizing the acquisition functions.
Refer to caption: x3.png

Figure: Figure 4: Inner-loop optimization results by rootfinding, a gradient-based multistart optimizer with random starting points (random multistart), and a genetic algorithm for (a) the 10D Levy, (b) 16D Ackley, and (c) 16D Powell functions. The plots are cumulative values of optimized GP-TS acquisition functions $\alpha_{k}^{\star}$ and cumulative CPU times $t_{k}$ for optimizing the acquisition functions.
Refer to caption: x4.png

#### Optimizing GP-TS Acquisition Functions via Rootfinding.

We assess the quality of solutions and computational cost for the inner-loop optimization of GP-TS acquisition functions by the proposed global optimization algorithm, termed rootfinding hereafter.
We do so by computing the optimized values $\alpha_{k}^{\star}$ of the GP-TS acquisition functions, the corresponding solution points $\mathbf{x}_{k}^{\star}$, and the CPU times $t_{k}$ required for optimizing the acquisition functions during the optimization process.
For low-dimensional problems of the 2D Schwefel and 4D Rosenbrock functions, we also compute the exact global solution points $\mathbf{x}_{k}^{\mathrm{t}}$ of the GP-TS acquisition functions by starting the gradient-based optimizer at a large number of initial points (set as $10^{4}$), which is much greater than the maximum possible number of starting points set for TS-roots.
For comparison, we extend the same GP-TS acquisition functions to inner-loop optimization using a gradient-based multistart optimizer with random starting points (i.e., random multistart) and a genetic algorithm.
In each outer-loop optimization iteration, the number of starting points for the random multistart and the population size of the genetic algorithm are equal to the number of starting points recommended for rootfinding. The same termination conditions are used for the three algorithms.

Figure 3 shows the comparative performance of the inner-loop optimization for low-dimensional cases:
the 2D Schwefel and 4D Rosenbrock functions.
We see that the optimized acquisition function values and the optimization runtimes by rootfinding and the random multistart algorithm are almost identical, both of which are much better than those by the genetic algorithm. Rootfinding gives the best quality of the new solution points in both cases, while the genetic algorithm gives the worst. Notably in higher-dimensional settings of the 10D Levy, 16D Ackley, and 16D Powell functions shown in Figure 4, rootfinding performs much better than the random multistart and genetic algorithm in terms of optimized acquisition values and optimization runtimes, which verifies the importance of the judicious selection of starting points for global optimization of the GP-TS acquisition functions and the efficiency of rootfinding in high dimensions.
The random multistart becomes worse in high dimensions.

#### TS-roots to Information-Theoretic Acquisition Functions.

We show how TS-roots can enhance the performance of MES [], which uses information about the maximum function value $f^{\star}$ for processing BO.
One approach to computing MES generates a set of GP posterior samples using TS-RF
and subsequently optimizes the generated functions for samples of $f^{\star}$
using a gradient-based multistart optimizer with a large number of random starting points [].
We hypothesize that high-quality $f^{\star}$ samples can improve the performance of MES.
Thus, we assign both TS-roots and TS-RF as the inner workings of MES and then compare the resulting optimal solutions.

Specifically, we minimize the 4D Rosenbrock, 6d Hartmann, and 10D Levy functions using four versions of MES, namely MES-R 10, MES-R 50, MES-TS-roots 10, and MES-TS-roots 50.
Here, MES-R [] and MES-TS-roots correspond to TS-RF and TS-roots, respectively, while 10 and 50 represent the number of random samples $f^{\star}$ generated for computing the MES acquisition function in each iteration.

Figure: Figure 5: Performance of MES-R 10 and MES-R 50 for (a) the 4D Rosenbrock function, (b) the 6d Hartmann function, and (c) the 10D Levy function when TS-RF and TS-roots are used for generating random samples from $f^{\star}|\mathcal{D}$. The plots are histories of medians and interquartile ranges of solutions from ten runs of each method.
Refer to caption: x5.png

Figure 5 shows the optimization histories for ten independent runs of each MES method.
On the 4D Rosenbrock and 6d Hartmann functions, MES with TS-roots demonstrates superior optimization performance and faster convergence compared to MES with TS-RF, especially when 50 samples of $f^{\star}$ are generated. For the 10D Levy function, TS-roots outperforms TS-RF when using 10 samples of $f^{\star}$, while their performances are comparable when 50 samples are used.

#### Sample-Average Posterior Function.

We finally propose a sample-average posterior function that explicitly controls exploration-exploitation balance
and, notably, can be generated at the cost of generating one posterior sample. For noiseless observations with $\widetilde{\omega}=\omega$,
we can rewrite the GP posterior function in Eq. 2 as
$\widetilde{f}_{\omega}(\mathbf{x})=f_{\omega}(\mathbf{x})+\widetilde{\mu}(
\mathbf{x})+\xi_{\omega}(\mathbf{x})$,
where $\xi_{\omega}(\mathbf{x})=-\boldsymbol{\kappa}_{\cdot,n}(\mathbf{x})\mathbf{K}_
{n,n}^{-1}\mathbf{f}_{n}$.
Define $\alpha_{\text{aTS}}(\mathbf{x})=\frac{1}{N_{\text{c}}}\sum_{j=1}^{N_{\text{c}}
}\widetilde{f}^{j}_{\omega}(\mathbf{x})$
as the sample-average posterior function, where $\widetilde{f}^{j}_{\omega}(\mathbf{x})$ are samples generated from the GP posterior
and $N_{\text{c}}\in\mathbb{N}_{>0}$. Since $\widetilde{\mu}(\mathbf{x})$ is deterministic,
and the scaled prior sample $\frac{1}{\sqrt{N_{\text{c}}}}f^{j}_{\omega}(\mathbf{x})$ can be written as
$\frac{1}{\sqrt{N_{\text{c}}}}f^{j}_{\omega}(\mathbf{x})\overset{\text{iid}}{
\sim}\mathcal{GP}(0,\frac{1}{N_{\text{c}}}\kappa(\mathbf{x},\mathbf{x}^{\prime
}))$,
we have $\alpha_{\text{aTS}}(\mathbf{x})=\mu(\mathbf{x})+\frac{1}{\sqrt{N_{\text{c}}}}
\left(f_{\omega}(\mathbf{x})+\xi_{\omega}(\mathbf{x})\right)$, where the first and second terms favor exploitation and exploration, respectively.
Thus, we can consider $N_{\text{c}}$ an exploration-exploitation control parameter that, at large values, prioritizes exploitation
by concentrating the conditional distribution of the global minimum location, i.e., $p(\mathbf{x}^{\star}|\mathcal{D})$,
at the minimum location of $\widetilde{\mu}(\mathbf{x})$, see Figure 8 in LABEL:{apd:AdditionalResults}.
With $\alpha_{\text{aTS}}(\mathbf{x})$, we can reproduce $\widetilde{f}_{\omega}(\mathbf{x})$ and the GP mean function $\widetilde{\mu}(\mathbf{x})$
by setting $N_{\text{c}}=1$ and $N_{\text{c}}=\infty$, respectively.

We investigate how $\alpha_{\text{aTS}}(\mathbf{x})$ influences the outer-loop optimization results.
For this, we set $N_{\text{c}}\in\{1,10,50,100\}$ for TS-roots to optimize the 2D Schwefel, 4D Rosenbrock, and 6d Ackley functions.
We observe that increasing $N_{\text{c}}$ from $1$ to $10$ improves TS-roots performance on the 2D Schwefel, 4D Rosenbrock, and 6d Ackley functions (see Figure 9 in Appendix G).
However, further increases in $N_{\text{c}}$ from 10 to 50 and 100 result in slight declines in solution quality as TS-roots transitions to exploitation.
These observations indicate that there is an optimal value of $N_{\text{c}}$ for each problem at which TS-roots achieves its best performance by balancing exploitation and exploration.
However, identifying such an optimal value to maximize the performance of $\alpha_{\text{aTS}}(\mathbf{x})$ for a particular optimization problem is an open issue.

## 6 Conclusion and Future Work

We presented TS-roots, a global optimization strategy for posterior sample paths.
It comprises an adaptive selection of starting points for gradient-based multistart optimizers, combining exploration with exploitation.
This strategy breaks the curse of dimensionality by exploiting the separability of Gaussian process priors.
Compared to random multistart and a genetic algorithm, TS-roots consistently yields higher-quality solutions
in optimizing posterior sample-based acquisition functions for Bayesian optimization, both in low- and high-dimensional settings.
It also improves the outer-loop optimization performance of GP-TS and information-theoretic acquisition functions such as MES.
For future work, we extend TS-roots to other spectral representations per Bochner’s theorem [, , ].
We also plan to study the modes and the probability when TS-roots fails to find the global optimum,
and the impact of subset sizes.

## Appendix A Characterizing the Local Minima of a Separable Function

### A.1 Proof of Proposition 1 : A Representation of the Set of Local Minima

Proposition 1 broadly applies to separable functions on a hypercube.
Consider a separable function $f(\mathbf{x})=\prod_{i=1}^{d}f_{i}(x_{i})$ with domain $\mathcal{X}=\prod_{i=1}^{d}[\underline{x}_{i},\overline{x}_{i}]$,
where $f_{i}\in C^{1}([\underline{x}_{i},\overline{x}_{i}];\mathbb{R})$.
To simplify the discussion, we further assume that $f_{i}$ is twice differentiable at its interior critical points $\mathring{\Xi}_{i}$.
The gradient of $f$ can be written as:

$$ $\nabla f(\mathbf{x})=\Big{(}f^{\prime}_{i}(x_{i})\cdot\prod_{j\neq i}f_{j}(x_{ j})\Big{)}_{i=1}^{d}=\bigg{(}f(\mathbf{x})\cdot\frac{f^{\prime}_{i}(x_{i})}{f_ {i}(x_{i})}\bigg{)}_{i=1}^{d}=f(\mathbf{x})\cdot\mathbf{v}(\mathbf{x}),$ (6) $$

where $\mathbf{v}(\mathbf{x})=\big{(}f^{\prime}_{i}/f_{i}\big{)}_{i=1}^{d}=\big{(}
\frac{\,\mathrm{d}}{\,\mathrm{d}x_{i}}\log f_{i}\big{)}_{i=1}^{d}$.
The Hessian of $f$ can be written as:

$$ $\nabla^{2}f(\mathbf{x})=\operatorname{diag}\Big{\{}f^{\prime\prime}_{i}(x_{i}) \prod_{j\neq i}f_{j}(x_{j})\Big{\}}_{i=1}^{d}+\Big{[}f^{\prime}_{i}(x_{i})f^{ \prime}_{j}(x_{j})\prod_{k\neq i,j}f_{k}(x_{k})\Big{]}_{i\in d}^{j\neq i}=f( \mathbf{x})\operatorname{diag}(\mathbf{s}+\mathbf{v}\mathbf{v}^{\intercal}),$ (7) $$

where $\mathbf{s}(\mathbf{x})=\big{(}f^{\prime\prime}_{i}/f_{i}-(f^{\prime}_{i}/f_{i}
)^{2}\big{)}_{i=1}^{d}=\big{(}\frac{\,\mathrm{d}^{2}}{\,\mathrm{d}x_{i}^{2}}
\log f_{i}\big{)}_{i=1}^{d}$.

An interior point $\mathbf{x}\in\operatorname*{int}\mathcal{X}:=\prod_{i=1}^{d}(\underline{x}_{i}
,\overline{x}_{i})$ is a strong local minimum of $f$
if and only if $\nabla f(\mathbf{x})=0$ and $\nabla^{2}f(\mathbf{x})>0$.
From Eq. 6, the first condition is satisfied in any of the following three cases:
(1) $f_{i}(x_{i})\neq 0$and $f^{\prime}_{i}(x_{i})=0$ for all $i\in\{1,\cdots,d\}$;
(2) $f_{i}(x_{i})=0$for exactly one $i\in\{1,\cdots,d\}$ and $f^{\prime}_{i}(x_{i})=0$; or
(3) $f_{i}(x_{i})=0$for all $i\in I\subseteq\{1,\cdots,d\}$ where $|I|\geq 2$.

In case (1), the Hessian Eq. 7 reduces to
$\nabla^{2}f(\mathbf{x})=f(\mathbf{x})\cdot\operatorname{diag}\{f^{\prime\prime
}_{i}(x_{i})/f_{i}(x_{i})\}_{i=1}^{d}$,
which is positive definite if and only if one of the following holds:
(i) $f(\mathbf{x})>0$and $f_{i}(x_{i})f^{\prime\prime}_{i}(x_{i})>0$, for all $i\in\{1,\cdots,d\}$; or
(ii) $f(\mathbf{x})<0$and $f_{i}(x_{i})f^{\prime\prime}_{i}(x_{i})<0$, for all $i\in\{1,\cdots,d\}$.

In case (2), the Hessian reduces to an all-zero matrix except for the $i$th diagonal entry:
$[\nabla^{2}f(\mathbf{x})]_{i,i}=f^{\prime\prime}_{i}(x_{i})\prod_{j\neq i}f_{j
}(x_{j})$.
Even if this entry is positive, the Hessian is still positive semi-definite,
which means that there is a continuum of weak local minima: $\{x_{i}\}\times\prod_{j\neq i}[\underline{x}_{j},\overline{x}_{j}]$.
Besides, this case requires $f_{i}$ and $f^{\prime}_{i}$ to have an identical root, which an event with probability zero.

In case (3), let $g_{i}(r_{i}):=f_{i}(x_{i}+r_{i})$ be a shifted version of $f_{i}$, $i\in\{1,\cdots,d\}$.
Taylor expansion at $\mathbf{r}=0$ gives $g_{i}(r_{i})=0+g^{\prime}_{i}(0)\,r_{i}+o(r_{i})$ for all $i\in I$
and $g_{j}(r_{j})=g_{j}(0)+O(r_{j})$ for all $j\notin I$.
We have $g(\mathbf{r}):=\prod_{i=1}^{d}g_{i}(r_{i})=c\prod_{i\in I}r_{i}+o(\prod_{i\in I
}r_{i})\cdot O(\prod_{j\notin I}r_{j})$,
where $c=\prod_{i\in I}g^{\prime}_{i}(0)\cdot\prod_{j\notin I}g_{j}(0)\neq 0$.
This means that there is a continuum of saddle points: $\{x_{i}\}_{i\in I}\times\prod_{j\notin I}[\underline{x}_{j},\overline{x}_{j}]$.

For a boundary point $\mathbf{x}\in\partial\mathcal{X}:=\mathcal{X}\setminus\operatorname*{int}
\mathcal{X}$,
we partition the index set $\{1,\cdots,d\}$ into $L,R$, and $I$ such that
$x_{i}=\underline{x}_{i}$ for all $i\in L$, $x_{i}=\overline{x}_{i}$ for all $i\in R$, and $x_{i}\in(\underline{x}_{i},\overline{x}_{i})$ for all $i\in I$.
Define $\nabla_{J}:=(\partial_{j})_{j\in J}$ for any subset $J$ of the indices.
Then $\mathbf{x}$ is a strong local minimum of $f$ if and only if the following conditions hold:
(a) $\mathbf{x}$is a strong local minimum in $\{x_{j}\}_{j\notin I}\times\prod_{i\in I}[\underline{x}_{j},\overline{x}_{j}]$;
(b) $\nabla_{L}f(\mathbf{x})>0$; and
(c) $\nabla_{R}f(\mathbf{x})<0$.

Condition (a) holds if any only if $\nabla_{I}f(\mathbf{x})=0$ and $\nabla_{I}^{2}f(\mathbf{x})>0$.
Based on the previous discussion on interior local minima, it is equivalent to:
(i) $f(\mathbf{x})>0$and $f_{i}(x_{i})f^{\prime\prime}_{i}(x_{i})>0$, for all $i\in I$; or
(ii) $f(\mathbf{x})<0$and $f_{i}(x_{i})f^{\prime\prime}_{i}(x_{i})<0$, for all $i\in I$.

From Eq. 6, condition (b) is equivalent to:
(i) $f(\mathbf{x})>0$and $f_{i}(x_{i})f^{\prime}_{i}(x_{i})>0$, for all $i\in L$; or
(ii) $f(\mathbf{x})<0$and $f_{i}(x_{i})f^{\prime}_{i}(x_{i})<0$, for all $i\in L$.

Similarly, condition (c) is equivalent to:
(i) $f(\mathbf{x})>0$and $-f_{i}(x_{i})f^{\prime}_{i}(x_{i})>0$, for all $i\in R$; or
(ii) $f(\mathbf{x})<0$and $-f_{i}(x_{i})f^{\prime}_{i}(x_{i})<0$, for all $i\in R$.

Summarizing the above discussions, we see that there is a unified way to identify the set $\breve{X}$
of all strong local minima of $f$, which is stated in Proposition 1.
The discussion for the set $\widehat{X}$ of local maxima is the exactly the same, except that the signs are flipped.
This also means that $\widehat{X}$ and $\breve{X}$ form a partition of the union $\Xi^{(0)}\sqcup\Xi^{(1)}$
of the two tensor grids.

If $f_{i}$ is not twice differentiable at some interior critical point $x_{i}$,
we may replace $f^{\prime\prime}_{i}(x_{i})>0$ with the statement that $x_{i}$ is a strong local minimum of $f_{i}$,
and replace $f^{\prime\prime}_{i}(x_{i})<0$ with the statement that $x_{i}$ is a strong local maximum of $f_{i}$.
The rest of the discussion still follows.
In practice, the differentiability of the prior sample is not an issue,
because it is almost always approximated by a finite sum of analytic functions, which is again analytic.

### A.2 Number of Local Minima of a Separable Function

In Proposition 1, each set of candidate coordinates $\Xi_{i}$ is partitioned into
mono type and mixed type:

$$ $\Xi_{i}^{(1)}=\{\xi_{i,j}\in\Xi_{i}:f_{i}(\xi_{i,j})h_{i}(\xi_{i,j})<0\},\quad \Xi_{i}^{(0)}=\{\xi_{i,j}\in\Xi_{i}:f_{i}(\xi_{i,j})h_{i}(\xi_{i,j})>0\}.$ $$

Another partition of $\Xi_{i}$ is by the sign of the corresponding component function value:

$$ $\Xi_{i}^{-}=\{\xi_{i,j}\in\Xi_{i}:f_{i}(\xi_{i,j})<0\},\quad\Xi_{i}^{+}=\{\xi_ {i,j}\in\Xi_{i}:f_{i}(\xi_{i,j})>0\}.$ $$

These two partitions create a finer partition of $\Xi_{i}$ into four subsets:

$$ $\Xi_{i}^{-(1)}=\Xi_{i}^{-}\cap\Xi_{i}^{(1)},\quad\Xi_{i}^{-(0)}=\Xi_{i}^{-} \cap\Xi_{i}^{(0)},\quad\Xi_{i}^{+(1)}=\Xi_{i}^{+}\cap\Xi_{i}^{(1)},\quad\Xi_{i }^{+(0)}=\Xi_{i}^{+}\cap\Xi_{i}^{(0)}.$ $$

Denote the sizes of mixed and mono type candidate coordinates as
$n_{i}^{(1)}=|\Xi_{i}^{(1)}|$ and $n_{i}^{(0)}=|\Xi_{i}^{(0)}|$,
then the sizes of the two tensor grids $\Xi^{(1)}$ and $\Xi^{(0)}$ can be written as:

$$ $N^{(1)}:=|\Xi^{(1)}|=\prod_{i=1}^{d}n_{i}^{(1)},\quad N^{(0)}:=|\Xi^{(0)}|= \prod_{i=1}^{d}n_{i}^{(0)}.$ $$

Define signed sums as the sums of signs of function values on the two tensor grids:

$$ $S^{(1)}:=\sum_{\boldsymbol{\xi}\in\Xi^{(1)}}\operatorname{sign}(f(\boldsymbol{ \xi})),\quad S^{(0)}:=\sum_{\boldsymbol{\xi}\in\Xi^{(0)}}\operatorname{sign}(f (\boldsymbol{\xi})).$ $$

We now derive efficient formulas to calculate these signed sums, using $S^{(1)}$ as an example.
Denote each coordinate in $\Xi_{i}^{(1)}$ as $\xi_{i,j}^{(1)}$.
Denote each point in $\Xi^{(1)}$ as $\boldsymbol{\xi}_{J}^{(1)}=(\xi_{i,J_{i}}^{(1)})_{i=1}^{d}$,
where multi-index $J=(J_{i})_{i=1}^{d}\in\Pi^{(1)}:=\prod_{i=1}^{d}\{1,\cdots,n_{i}^{(1)}\}$.
The signed sum $S^{(1)}$ can be written as:

$$ $\displaystyle S^{(1)}$ $\displaystyle=\sum_{J\in\Pi^{(1)}}\operatorname{sign}(f(\boldsymbol{\xi}_{J}^{ (1)}))=\sum_{J\in\Pi^{(1)}}\operatorname{sign}\bigg{(}\prod_{i=1}^{d}f_{i}(\xi _{i,J_{i}}^{(1)})\bigg{)}$ $\displaystyle=\sum_{J\in\Pi^{(1)}}\prod_{i=1}^{d}\operatorname{sign}(f_{i}(\xi _{i,J_{i}}^{(1)}))=\prod_{i=1}^{d}\sum_{j=1}^{n_{i}^{(1)}}\operatorname{sign}( f_{i}(\xi_{i,j}^{(1)}))$ $\displaystyle=\prod_{i=1}^{d}\left[\sum_{j=1}^{n_{i}^{(1)}}1(f_{i}(\xi_{i,j}^{ (1)})>0)-\sum_{j=1}^{n_{i}^{(1)}}1(f_{i}(\xi_{i,j}^{(1)})<0)\right]=\prod_{i=1 }^{d}\left[|\Xi_{i}^{+(1)}|-|\Xi_{i}^{-(1)}|\right].$ $$

A formula for $S^{(0)}$ can be derived analogously.
Denote set sizes:

$$ $n_{i}^{-(1)}=|\Xi_{i}^{-(1)}|,\quad n_{i}^{-(0)}=|\Xi_{i}^{-(0)}|,\quad n_{i}^ {+(1)}=|\Xi_{i}^{+(1)}|,\quad n_{i}^{+(0)}=|\Xi_{i}^{+(0)}|,$ $$

then the signed sums can be calculated as:

$$ $S^{(1)}=\prod_{i=1}^{d}(n_{i}^{+(1)}-n_{i}^{-(1)}),\quad S^{(0)}=\prod_{i=1}^{ d}(n_{i}^{+(0)}-n_{i}^{-(0)}).$ $$

The sizes of negative and positive strong local minima of a separable function can be written as:

$$ $\displaystyle\breve{N}^{-}:=|\breve{X}^{-}|=\sum_{\boldsymbol{\xi}\in\Xi^{(1)} }1(f(\boldsymbol{\xi})<0)=\frac{1}{2}(N^{(1)}-S^{(1)}),$ (8) $\displaystyle\breve{N}^{+}:=|\breve{X}^{+}|=\sum_{\boldsymbol{\xi}\in\Xi^{(0)} }1(f(\boldsymbol{\xi})>0)=\frac{1}{2}(N^{(0)}+S^{(0)}).$ $$

Therefore, the size of the strong local minima of a separable function can be written as:

$$ $\breve{N}:=|\breve{X}|=|\breve{X}^{-}|+|\breve{X}^{+}|=\frac{1}{2}(N^{(1)}+N^{ (0)}-S^{(1)}+S^{(0)}).$ (9) $$

## Appendix B Ordering the Local Minima of a Separable Function

### B.1 Filtering a Tensor Grid for High Absolute Values of a Separable Function

The step one in Section 3.4 is equivalent to the following problem:
given coordinates $Z_{i}=\{\zeta_{i,1},\cdots,\zeta_{i,t_{i}}\}$
and components values $F_{i}=\{f_{i,1},\cdots,f_{i,t_{i}}\}$, $i\in\{1,\cdots,d\}$,
of a separable function $f(\mathbf{x})=\prod_{i=1}^{d}f_{i}(x_{i})$,
find points $\boldsymbol{\zeta}$ such that $|f(\boldsymbol{\zeta})|$ are the $k$ largest
in the tensor grid $Z=\prod_{i=1}^{d}Z_{i}$.

Because $\log|f(\mathbf{x})|=\log|\prod_{i=1}^{d}f_{i}(x_{i})|=\sum_{i=1}^{d}\log|f_{i}
(x_{i})|$,
we can solve this problem as follows:
define two-dimensional arrays $F=[F_{1},\cdots,F_{d}]$ and $A=\log|F|$,
solve $S=\texttt{maxk_sum}(A,k)$,
and return $\{\boldsymbol{\zeta}=(\zeta_{1,I_{1}},\cdots,\zeta_{d,I_{d}}):I\in S\}$.
Here the maxk_sum algorithms finds the combinations from $A$ that gives the $k$ largest sums,
which is described next.

### B.2 Top Combinations with the Largest Sums

Consider this problem: given a two-dimensional array $A=[\mathbf{a}_{1},\cdots,\mathbf{a}_{d}]$,
$\mathbf{a}_{i}=[a_{i,1},\cdots,a_{i,t_{i}}]$, with $a_{i,1}\geq\cdots\geq a_{i,t_{i}}$, $i\in\{1,\cdots,d\}$,
find $k$ multi-indices of the form $I=[I_{1},\cdots,I_{d}]$ such that the sums $s_{I}:=\sum_{i=1}^{d}a_{i,I_{i}}$
are the $k$ largest among all combinations $I\in\prod_{i=1}^{d}\{1,\cdots,t_{i}\}$.

An exhaustive search is intractable because the number of all possible combinations
grows exponentially as $\prod_{i=1}^{d}t_{i}$.
Instead, we use a min-heap to efficiently keep track of the top $k$ combinations.
A min-heap is a complete binary tree, where each node is no greater than its children.
The operations of inserting an element and removing the smallest element from a min-heap can be done in logarithmic time.
Algorithm 1 gives a procedure to solve the above problem using min-heaps.

Figure: Algorithm 1 maxk_sum: Combinations with the $k$ largest sums

This algorithm has time complexity $\mathcal{O}(tk\log k)$, where $t=\sum_{i=1}^{d}t_{i}\ll\prod_{i=1}^{d}t_{i}$,
and space complexity $\mathcal{O}(dk)$.
In TS-roots, the cost of maxk_sum is small compared with the gradient-based multistart optimization of the posterior sample.

## Appendix C Algorithms for TS-roots

### C.1 Spectral Sampling of Separable Gaussian Process Priors

Per Mercer’s theorem on probability spaces (see e.g., [], Sec 4.3),
any positive definite covariance function that is essentially bounded
with respect to some probability measure $\mu$ on a compact domain $\mathcal{X}$
has a spectral representation $\kappa(\mathbf{x},\mathbf{x}^{\prime})=\sum_{k=0}^{\infty}\lambda_{k}\phi_{k}(
\mathbf{x})\phi_{k}(\mathbf{x}^{\prime})$,
where $(\lambda_{k},\phi_{k}(\mathbf{x}))$ is a pair of eigenvalue and eigenfunction of the kernel integral operator.
The corresponding GP prior can be written as
$f_{\omega}(\mathbf{x})=\sum_{k=0}^{\infty}w_{k}\sqrt{\lambda_{k}}\phi_{k}(
\mathbf{x})$,
where $w_{k}\overset{\text{iid}}{\sim}\mathcal{N}(0,1)$ are independent standard Gaussian random variables.
Similar spectral representations exist per Bochner’s theorem, which may have efficient discretizations [, ].

Given spectral representations of the univariate component functions of a separable Gaussian Process prior,
we can accurately approximate the prior sample as

$$ $f_{\omega}(\mathbf{x})=\prod_{i=1}^{d}f_{i}(x_{i};\omega_{i}),\quad f_{i}(x_{i };\omega_{i})\approx\sum_{k=0}^{N_{i}-1}w_{i,k}\sqrt{\lambda_{i,k}}\phi_{i,k}( x_{i}).$ (10) $$

Here $N_{i}$ is selected for each variate such that $\lambda_{i,N_{i}-1}/\lambda_{i,1}\leq\eta_{i}$,
where $\eta_{i}$ is sufficiently small (see Appendix F for the value used in this study).
Using spectral representations of the univariate components as in Eq. 10
is much more efficient than directly using a spectral representation of the separable GP prior,
because the former uses $\sum_{i=1}^{d}N_{i}$ univariate terms to
exactly represent $\prod_{i=1}^{d}N_{i}$ multivariate terms in the latter.

#### Spectrum of the Squared Exponential Covariance Function.

The univariate squared exponential (SE) covariance function can be written as
$\kappa(x,x^{\prime};l)=\exp(-\frac{1}{2}s^{2})$, where the relative distance $s=|x-x^{\prime}|/l$.
The spectral representation of such a covariance function per Mercer’s theorem is $\kappa(x,x^{\prime})=\sum_{k=0}^{\infty}\lambda_{k}\phi_{k}(x)\phi_{k}(x^{
\prime})$.
With a Gaussian measure $\mu=\mathcal{N}(0,\sigma^{2})$ over the domain $\mathcal{X}=\mathbb{R}$,
we can write the eigenvalues $\lambda_{k}$ and eigenfunctions $\phi_{k}(x)$ of the kernel integral operator as follows.
(See e.g., [] Sec. 4 and [] 7.374 eq. 8.)

Define constants $a=(2\sigma^{2})^{-1}$, $b=(2l)^{-1}$, $c=\sqrt{a^{2}+4ab}$, and
$A=\frac{1}{2}a+b+\frac{1}{2}c$.
For $k\in\mathbb{N}$, the $k$th eigenvalue is $\lambda_{k}=\sqrt{\frac{a}{A}}\left(\frac{b}{A}\right)^{k}$ and the corresponding eigenfunction is
$\phi_{k}(x)=\left(\frac{\pi c}{a}\right)^{1/4}\psi_{k}(\sqrt{c}x)\exp\left(
\frac{1}{2}ax^{2}\right)$,
where $\psi_{k}(x)=\left(\pi^{1/2}2^{k}k!\right)^{-1/2}H_{k}(x)\exp\left(-\frac{1}{2}
x^{2}\right)$ and $H_{k}(x)$ the $k$th-order Hermite polynomial defined by $H_{k}(x)=(-1)^{k}\exp(x^{2})\frac{d^{k}}{dx^{k}}\exp(-x^{2})$.

Figure 6 shows approximations to the SE covariance function
by truncated spectral representations with the first $N$ eigenpairs
and by random Fourier features [] with $N$ basis functions.
The spectral representation per Mercer’s theorem converges quickly to the true covariance function,
while the random Fourier features representation requires a large number of basis functions and is inaccurate for $N<1000$.

Figure: Figure 6: Approximate SE covariance functions from (a) the spectral representation per Mercer’s theorem with the first $N$ eigenpairs and (b) the random Fourier features representation with $N$ basis functions. The plots are generated for $l=1$.
Refer to caption: x6.png

### C.2 Univariate Global Rootfinding

Algorithm 2 outlines a method for univariate global rootfinding on an interval by solving an eigenvalue problem.
When the orthogonal polynomial basis is the Chebyshev polynomials, the corresponding comrade matrix is called a colleague matrix,
and we have the following theorem:

Figure: Algorithm 2 roots: Univariate global rootfinding on an interval.

###### Theorem 1 .

Let $p(x)=\sum_{k=0}^{m}a_{k}T_{k}(x)$, $a_{m}\neq 0$, be a polynomial of degree $m$,
where $T_{k}$ is the $k$th Chebyshev polynomial and $a_{k}$ is the corresponding weight.
The roots of $p(x)$ are the eigenvalues of the following $m\times m$ colleague matrix:

$$ $C=\left(\begin{matrix}0&1&&&&\\ 1/2&0&1/2&&&\\ &1/2&0&1/2&&\\ &&\ddots&\ddots&\ddots&\\ &&&&&1/2\\ &&&&1/2&0\end{matrix}\right)-\frac{1}{2a_{m}}\left(\begin{matrix}&&&&&\\ &&&&&\\ &&&&&\\ &&&&&\\ &&&&&\\ a_{0}&a_{1}&a_{2}&\cdots&\cdots&a_{m-1}\end{matrix}\right),$ (11) $$

where the elements not displayed are zero.

A proof of Theorem 1 is provided in [], Chapter 18.
A classical formula to compute the weights $\{a_{k}\}$ requires $\mathcal{O}{(m^{2})}$ floating point operations,
which can be reduced to $\mathcal{O}{(m\log m)}$ using a fast Fourier transform.
Since the colleague matrix is tridiagonal except in the final row,
the complexity of computing its eigenvalues can be improved from $\mathcal{O}{(m^{3})}$ to $\mathcal{O}{(m^{2})}$ operations,
which can be further improved to $\mathcal{O}{(m)}$ via recursive subdivision of intervals (see Trefethen []).
The roots algorithm is implemented in the Chebfun package in MATLAB []
and the chebpy package in Python [];
both packages also implement other related programs such as chebfun for Chebyshev polynomial approximation
and diff for differentiation.

Figure: Algorithm 3 minsort: Best local minima of a separable function

###### Theorem 1 .

### C.3 Best Local Minima of a Separable Function

Given the univariate component functions of a separable function,
Algorithm 3 finds the subset $S_{\text{o}}$ of the local minima of the function with the $n_{o}$ smallest function values.
This procedure requires the maxk_sum algorithm in Algorithm 1,
the roots algorithm in Algorithm 2 and the related programs chebfun and diff,
see also Appendix F.

In Algorithm 3, $\boldsymbol{\xi},\mathbf{f},\mathbf{h},J,P$ are two-dimensional arrays,
while $\mathbf{I},\boldsymbol{\Pi}^{(1)},\boldsymbol{\Pi}^{(0)}$ are matrices.
Function evaluations at
10, 11, 15, 20, 21 and 25
are only notational: the sign and value of the function can be computed efficiently
by multiplying the signs and values of its components at the selected coordinates.
For example, the statement $f(\boldsymbol{\xi}^{(1)}(\mathbf{I}))<0$ at 10 can be evaluated as
$\texttt{rowXor}(P^{(1)}(\mathbf{I}))$,
where $P^{(1)}$ is a two-dimensional array with $P_{i}^{(1)}=P_{i}(\neg J_{i})$,
$P^{(1)}(\mathbf{I})$ is a matrix with $d$ columns,
and rowXor is row-wise exclusive or operation.
Similarly, the statement $f(\boldsymbol{\xi}^{(1)}(\mathbf{I}))$ at 11 can be evaluated as
$\texttt{rowProd}(\mathbf{f}^{(1)}(\mathbf{I}))$, where rowProd is row products.

### C.4 Decoupled Sampling from Gaussian Process Posteriors

The decoupled sampling method for GP posteriors [],
together with the spectral sampling of separable GP priors, is outlined in Algorithm 4.

Figure: Algorithm 4 Decoupled sampling of Gaussian process posterior

Figure: Algorithm 5 Sequential optimization []

Figure: Algorithm 6 Bayesian optimization policy

Figure: Algorithm 7 Gaussian process Thompson sampling acquisition function

## Appendix D Bayesian Optimization via Thompson Sampling

A general procedure for sequential optimization is given in Algorithm 5.
The initial dataset $\mathcal{D}^{0}$ can either be empty or contain some observations.
In the latter case we can write $\mathcal{D}^{0}=\{(\mathbf{x}^{i},y^{i})\}_{i=1}^{n_{0}}$, where $n_{0}\in\mathbb{N}_{>0}$.
Three components of this algorithm can be customized:
the observation model $\mathrm{Observe}(\mathbf{x})$, the optimization policy $\mathrm{Policy}(\mathcal{D})$,
and the termination condition.

BO can be seen as an optimization policy for sequential optimization.
A formal procedure is given in Algorithm 6.
Three components of this algorithm can be customized: the prior probabilistic model $f$,
the acquisition function $\alpha$, and the global optimization algorithm.
Any probabilistic model of the objective function $f_{\text{true}}$ can be seen as a probability distribution on a function space,
and the prior $f$ is usually specified as a stochastic process such as a GP.
The acquisition function $\alpha$ derived from the posterior $f|\mathcal{D}$
can be either deterministic—such as EI and LCB—or stochastic, such as GP-TS.
To simplify notation, we state the global optimization problem of
$\alpha(\mathbf{x})$ as minimization rather than maximization.
The two problems are the same with a change of sign to the objective.

When applied to BO, GP-TS generates a random acquisition function simply by sampling the posterior model.
That is, given the posterior $f^{k}$ at the $k$th BO iteration,
the GP-TS acquisition function is a random function:
$\alpha^{k}(\mathbf{x})\sim f^{k}$.

## Appendix E Benchmark Functions

The analytical expressions for the benchmark functions used in Section 5 are given below. The global solutions of these functions are detailed in [].

#### Schwefel Function:

$$ $f(\mathbf{x})=418.9829d-\sum_{i=1}^{d}x_{i}\sin\left(\sqrt{|x_{i}|}\right).$ (12) $$

This function is evaluated on $\mathcal{X}=[-500,500]^{d}$
and has a global minimum $f^{\star}:=f(\mathbf{x}^{\star})=0$ at $\mathbf{x}^{\star}=[420.9687,\cdots,420.9687]^{\intercal}$.
This function is $C^{1}$ at $\mathbf{x}=0$.

#### Rosenbrock Function:

$$ $f(\mathbf{x})=\sum_{i=1}^{d-1}\left[100(x_{i+1}-x_{i}^{2})^{2}+(x_{i}-1)^{2} \right].$ (13) $$

This function is evaluated on $\mathcal{X}=[-5,10]^{d}$
and has a global minimum $f^{\star}=0$ at $\mathbf{x}^{\star}=[1,\cdots,1]^{\intercal}$.

#### Levy Function:

$$ $f(\mathbf{x})=\sin^{2}(\pi w_{1})+\sum_{i=1}^{d-1}(w_{i}-1)^{2}\left[1+10\sin^ {2}(\pi w_{i}+1)\right]+(w_{d}-1)^{2}\left[1+\sin^{2}(2\pi w_{d})\right],$ (14) $$

where $w_{i}=1+\frac{x_{i}-1}{4}$, $i=1,\cdots,d$.
This function is evaluated on $\mathcal{X}=[-10,10]^{d}$
and has a global minimum $f^{\star}=0$ at $\mathbf{x}^{\star}=[1,\cdots,1]^{\intercal}$.

#### Ackley Function:

$$ $f(\mathbf{x})=-a\exp\left(-b\sqrt{\frac{1}{d}\sum_{i=1}^{d}x^{2}_{i}}\right)- \exp\left(\frac{1}{d}\sum_{i=1}^{d}\cos(cx_{i})\right)+a+\exp(1),$ (15) $$

where $a=20$, $b=0.2$, and $c=2\pi$.
This function is evaluated on $\mathcal{X}=[-10,10]^{d}$
and has a global minimum $f^{\star}=0$ at $\mathbf{x}^{\star}=[0,\cdots,0]^{\intercal}$.
This function not differentiable at $\mathbf{x}^{\star}$.

#### Powell Function:

$$ $f(\mathbf{x})=\sum_{i=1}^{d/4}\left[\left(x_{4i-3}+10x_{4i-2}\right)^{2}+5 \left(x_{4i-1}-x_{4i}\right)^{2}+\left(x_{4i-2}-2x_{4i-1}\right)^{4}+10\left(x _{4i-3}-x_{4i}\right)^{4}\right].$ (16) $$

This function is evaluated on $\mathcal{X}=[-4,5]^{d}$
and has a global minimum $f^{\star}=0$ at $\mathbf{x}^{\star}=[0,\cdots,0]^{\intercal}$.

#### 6d Hartmann Function:

$$ $f(\mathbf{x})=-\sum_{i=1}^{4}a_{i}\exp\left(-\sum_{j=1}^{6}A_{ij}(x_{j}-P_{ij} )^{2}\right),$ (17) $$

where

$$ $\mathbf{a}=[1,1.2,3,3.2]^{\intercal},$ (18a) $\mathbf{A}=\begin{bmatrix}10&3&17&3.5&1.7&8\\ 0.05&10&17&0.1&8&14\\ 3&3.5&1.7&10&17&8\\ 17&8&0.05&10&0.1&14\end{bmatrix},$ (18b) $\mathbf{P}=10^{-4}\begin{bmatrix}1312&1696&5569&124&8283&5886\\ 2329&4135&8307&3736&1004&9991\\ 2348&1451&3522&2883&3047&6650\\ 4047&8828&8732&5743&1091&381\end{bmatrix}$ (18c) $$

This function is evaluated on $\mathcal{X}=[0,1]^{6}$
and has a global minimum $f^{\star}=-3.32237$ at $\mathbf{x}^{\star}=[0.20169,0.150011,0.476874,0.275332,0.311625,0.6573]^{\intercal}$.
The rescaled version $\widetilde{f}(\mathbf{x})=\frac{f(\mathbf{x})-2.58}{1.94}$ [] is used in the experiments.

## Appendix F Experimental Details

#### Data Generation.

We generate 20 initial datasets for each problem.
The input observations are randomly generated using the Latin hypercube sampling [] within $[-1,1]^{d}$,
where $d$ represents the number of input variables.
The normalized input observations are transformed into their real spaces to evaluate the corresponding objective function values
which are then standardized using the $z$-score for processing optimization.
Each BO method in comparison starts from each of the generated datasets.

#### Key Parameters for TS-roots and other BO Methods.

We use squared exponential (SE) covariance functions for our experiments.
The spectra of univariate SE covariance functions for all problems (see Section C.1)
are determined using the Gaussian measure $\mu=\mathcal{N}(0,1)$.
The number of terms $N_{i}$, $i\in\{1,\cdots,d\}$, of each truncated univariate spectrum is determined such that
$\lambda_{i,N_{i}-1}/\lambda_{i,1}\leq\eta_{i}$, where $\eta_{i}=10^{-16}$.
If $N_{i}>1000$, we set $N_{i}=1000$ to trade off between the accuracy of truncated spectra and computational cost.
The maximum size of the exploration set is $n_{\text{e}}=250$.
The maximum size of the exploitation set is $n_{\text{x}}=200$.

The number of initial observations is $10d$ for all problems.
The standard deviation of observation noise $\sigma_{\text{n}}=10^{-6}$ is applied for standardized output observations.
The number of BO iterations for the 2D Schwefel and 4D Rosenbrock functions is 200, while that for the 10D Levy, 16D Ackley, and 16D Powell functions is 800. Other GP-TS methods for optimization of benchmark test functions including TS-DSRF (i.e., TS using decoupled sampling with random Fourier
features) and TS-RF (i.e., TS using random Fourier features) are characterized by a total of 2000 random Fourier features.

To ensure a fair comparison of outer-optimization results, we first implement TS-roots and record the number of starting points used in each optimization iteration. We then apply other BO methods, each employing a gradient-based multistart optimizer with the same number of random starting points and identical termination criteria as those used for TS-roots in each iteration.

For the comparative inner-loop optimization performance of the proposed method via rootfinding with the random multistart and genetic algorithm approaches, we set the same termination tolerance on the objective function value as the stopping criterion for the methods. In addition, the number of starting points for the random multistart and the population size of the genetic algorithm are the same as the number of points in both the exploration and exploitation sets of rootfinding in each optimization iteration.

#### Computational Tools.

We carry out all experiments, except those for inner-loop optimization, using a designated cluster at our host institution. This cluster hosts 9984 Intel CPU cores and 327680 Nvidia GPU cores integrated within 188 compute and
20 GPU nodes.
The inner-loop optimization is implemented on a PC with an Intel^® Core^TM i7-1165G7 @ 2.80 GHz and 16 GB memory.

For the univariate global rootfinding via Chebyshev polynomials, we use MATLAB’s Chebfun package []
and its corresponding implementation in Python, called chebpy [].

## Appendix G Additional Results

#### Distance to Global Minimum.

Figure 7 shows the solution locations from 20 runs of TS-roots, TS-DSRF, TS-RF, EI, and LCB for the 2D Schwefel, 4D Rosenbrock, 10D Levy, 16D Ackley, 16D Powell functions.

#### Sample-average Posterior Function.

Figure 8 shows how we can improve the exploitation of GP-TS when increasing the exploration-exploitation control parameter $N_{\text{c}}$.

#### Performance of Sample-average TS-roots.

Figure 9 shows the performance of sample-average TS-roots with different exploration-exploitation control parameters $N_{\text{c}}$ for the 2D Schwefel, 4D Rosenbrock, and 6d Ackley functions.

Figure: Figure 7: Outer-loop optimization results for (a) the 2D Schwefel, (b) 4D Rosenbrock, (c) 10D Levy, (d) 16D Ackley, (e) 16D Powell functions. The plots are histories of medians and interquartile ranges of solution locations from 20 runs of TS-roots, TS-DSRF (i.e., TS using decoupled sampling with random Fourier features), TS-RF (i.e., TS using random Fourier features), EI, and LCB.
Refer to caption: x7.png

Figure: Figure 8: Sample-average posterior function for different values of $N_{\text{c}}$. The posterior function approaches the GP mean and the conditional distribution of the solution location $p(x^{\star}|\mathcal{D})$ is more concentrated when we increase $N_{\text{c}}$.
Refer to caption: x8.png

Figure: Figure 9: Performance of sample-average TS-roots with different control values $N_{\text{c}}$ for (a) the 2D Schwefel, 4D Rosenbrock, and (b) 6d Ackley functions. The plots are histories of medians and interquartile ranges of solution values and solution locations from 20 runs of TS-roots for each $N_{\text{c}}$ value.
Refer to caption: x9.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: -	There are also other works addressing the issue on mostly-flat acquisition functions (stated by authors in line 49). One of the notable approaches is to compute log of acquisition values, which has been shown to mitigate the flatness problem. One example is a recent work [1], which combines previous similar attempts into a systematic solution, the logEI acquisition function. Can the authors comment on this log approach? For example, how about computing log acquisition values + random multi-start for gradient solver?

[1] Ament, Sebastian, et al. "Unexpected improvements to expected improvement for bayesian optimization." Advances in Neural Information Processing Systems 36 (2023): 20577-20612.
soundness: 3
strengths: -	The paper addresses the problems of finding the global optimum of GP posterior samples efficiently, which is an important topic for many BO acquisition functions.
-	There are many empirical analyses for the proposed methods.
-	The method significantly reduces the CPU times to optimize GP-TS compared against genetic algorithm and random multi-start.
-	The method can be a robust choice for general problem – TS-root performance is the most robust among other baselines, despite only being on par with the best baselines in each problem.
summary: The paper proposes a method for optimizing the sample paths generated by GP posterior. The method, TS-root, globally optimizes the posterior samples via gradient-based multi-start optimization. 
Specifically, TS-root can provide the gradient-based solver with a set of starting points accounting for exploration and exploitation, thereby finding the global optimum of posterior sample more efficiently. This starting point set S is discussed in section 3.1. Basically, S consists of exploration set S_e and exploitation set S_x. As S_x is a subset of observed data points satisfying eq (4), the number of points n_x can be small (hence easy to find, I think?). However, the set S_e is a n_e subset from the set S_0 of n_0 prior sample local optimum, and because n_0 can be larger, especially in high dimensions, finding S_e is non-trivial, which is discussed in Section 3.4 and 3.5. The main focus of the paper is to compute set S_e.

Another important component of TS-root is that it does not follow the conventional representation for GP posterior, and instead adopt the Decoupled form of GP posterior. This technique was presented in Wilson et al., 2020. The main idea of this decoupled form is to reduce to computational cost of conventional GP posterior, while maintaining similar predictive performance.
weaknesses: -	The contribution is fair, as the main contribution seems to be computing the starting points for gradient-based optimization process. The decoupled representation for GP posterior was proposed by previous works.
-	The idea seems to not very effective in improving MES acquisition function. Fig. 5 shows that TS-Random Fourier and TS-root are statistically similar.
-	No real-world benchmark problems.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: 1. The author describes the computational complexity of the algorithm in the appendix. It seems that finding local optimizers of the prior sample path is of O(m^2). It seems that this additional computational cost may be a bit of problematic in high(or even moderate) dimension, where you need more initial locations to optimize the inner acquisition function. Can you provide some computational efficiency comparison with other existing methods (based on LHS or uniform grid) to perform posterior sample path based acquisition function optimization?

2. Rather subtle question I have is, the robustness of the separability assumption on Kernel. It seems that if one is only interested in characterizing the optimum, utilizing the separable kernel instead of standard ones (like standard Matern Kernel) is sufficient for describing the geometry near the optima. Can you comment a bit on the necessity of the separability assumption for Kernel, besides the implementation or computaion aspect?
soundness: 3
strengths: The paper delves into a problem that's important in Bayesian optimization but less understood, namely the effect or efficient strategy of optimizing the acquisition function. As the paper said, many (almost all) Bayesian optimization methods assume the exact inner optimization solution, which can be very difficult in practice. The authors provide a systematic way to come up with a starting point to facilitate existing multi-start gradient-based optimization strategies in acquisition functions. The effectiveness of this strategy and improving the existing sample path based acquisition function strategy looks impressive.
summary: The paper proposes a new methodology to optimize acquisition functions based on posterior sample paths. By judiciously selecting initial starting points, the proposed method makes a gradient-based multi-start strategy an appealing technique for optimizing the acquisition function. The effectiveness of the proposed methods has been demonstrated in multiple examples and comparisons with other existing methods are provided.
weaknesses: No major weakness. The paper is clearly written, except for the part where it describes the starting point for the inner acquisition optimization task. It seems that there are some inconsistencies between descriptions and notations between lines 142-161 on page 3. 1) I suggest authors revise this part to enhance the clarity of the method. 2) Furthermore, having an algorithmic box description would help readers gain the gist of the proposed method.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. There are some hyperparameters, e.g., n_o, n_e, n_x, that control the balance between exploration and exploitation. Given a specific optimization problem to be resolved, how to specify appropriate values for these hyperparameters? Is there any thoretical or empirical guidance?

2. The proposed method is claimed to be able to find the global optimum of the aquisition functions in a scalable way. Is there a theoretical guarantee?  It seems that the performance of the posterior sampling procedure plays a key role for the follow-up root finding based global optimum searching procedure. Is the posterior sampling procedure scalable to higher dimensions? As far as we know, posterior sampling of a high dimensional parameter space is much more challenging issue to be resolved. If this procedure only yields low quality samples, how to guarantee that the proposed method here can work well?
soundness: 2
strengths: 1. This paper selects an uncommon topic of BO, namely global optimization of acquisition functions, and conduct a systematic research on it.
2. The paper proposes a novel strategy that uses root finding to search the global optimum of the acquisition functions, and evaluates it by numerical experiments.
summary: This paper focuses on global optimization of acquisition functions in the framework of Bayesian optimization (BO). The common practice of BO is to utilize routine algorithms, such as dividing rectangles (DIRECT),  covariance matrix adaptation evolution strategy (CMA-ES), and Thompson sampling (TS), to solve the acquisition function optimization problem. The authors here propose to leverage root finding procedures to improve/accelerate optimization of posterior sample-based acquisition functions. The basic idea is to select a subset of posterior samples as the starting points to conduct gradient-based root finding for searching local minima of the acquisition function. The global optimum is obtained by ranking the local optima. The performance of the proposed methods is evaluated by comparison with some standard BO algorithms on several benchmark testing functions.
weaknesses: 1. The premise for that this paper is an important contribution to the community is that global optimization of acquisition functions is really necessary for guaranteeing BO's performance. Unfortunately this is not proved or even discussed in the paper. Is it worth to develop more complex algorithms which consume more computations to optimize the acquisition function? Is this better than we using the same resource to build a better surrogate model to approximate the objective function? Why the community paid more attention to the latter issue? It's likely because it plays a more important role for BO's performance.   

2. As shown in Figure 2, for some cases, the proposed method, TS-roots, performs better than othe competitors; while for the other cases, it just performs on par with or even worse than the standard commonly used baseline BO methods. In addition, only several simple benchmark test objective functions are involved in the experiments. Combining the above two points, the experiment results don't provide enough evidence for that the proposed solution can bring remarkable benefits for solving challenging real-life problems.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: * Is the method also applicable to Entropy Search [1]? The paragraph (l. 54-63) in the introduction mentions many information-theoretic acquisition functions besides this one and Entropy Search is arguably one of the main/first such acquisition functions. 

[1] Entropy Search for Information-Efficient Global Optimization (Hennig and Schuler, 2012)
soundness: 3
strengths: * The method is motivated well (Section 3.2)
* The method performs better than (or at least on-par with ) the baselines empirically. 
* I think it is good that also the outer-loop performance is tested.  At the first sight, I was wondering whether a potentially slight improvement in the inner loop would really matter for the final performance, but now I am more convinced that it does.
* Figure 1 is very useful
summary: The paper is about improving the inner-loop in Bayesian Optimization methods that rely on optimization of posterior samples. This is done by selecting better starting points for the gradient based optimization of the posterior samples. The set of starting points is made up from a subset of the local minima of the prior and the local minimal of the posterior mean function (the subset where the posterior values are highest, respectively). The empirical performance of the method is assessed on common benchmark functions.
weaknesses: * For me, a bit more explanation on the decouple representation of GP Posteriors would have been useful. There is a reference to Wilson et al. 2020, but if one is not familiar with that particular piece of work, a bit more explanation on the decouple representation of GP Posteriors would have been useful. In particular, because this representation later on is key in the construction of the method. 
* I think it would be nicer to either have Figure 8/Figure 9 from the Appendix near the "Sample-Average Posterior Function" Section. Or maybe even have the "Sample-Average Posterior Function" Section in the methods part instead of the results part? Otherwise, I reads more like an afterthought.

Meta-review: This study proposes a novel method for the global optimization of acquisition functions by optimizing sample paths generated from the posterior distribution of Gaussian processes (GPs). The key feature of this method is its utilization of the decoupled representation of the GP posterior. By leveraging this decoupled framework, the proposed method provides a systematic algorithm based on gradient optimization that promotes the use of multiple initial points.　A potential weakness is the assumption of separability. However, as the authors point out, this assumption is satisfied by the squared exponential kernel and can be met for the Matérn kernel with appropriate adjustments, making it unlikely to pose a significant practical limitation.
All reviewers have acknowledged the significant contributions of this paper, and numerical experiments further demonstrate its utility. Therefore, I believe this work is valuable to the community and recommend acceptance.

Decision: Accept (Poster)