Optimizing Posterior Samples for Bayesian Optimization via Rootfinding
Authors: Taiwo Adebiyi
Bach Do
Ruda Zhang
Venue: ICLR 2025 Poster (2025)
Keywords: Bayesian optimization
global optimization
acquisition function
Thompson sampling
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

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 sample paths, especially in higher dimensions. We introduce an efficient global optimization strategy for posterior samples based on global rootfinding. It provides gradient-based optimizers with two sets of judiciously selected starting points, designed to combine exploration and exploitation. The number of starting points can be kept small without sacrificing optimization quality. Remarkably, even with just one point from each set, the global optimum is discovered most of the time. The algorithm scales practically linearly to high dimensions, breaking the curse of dimensionality. For Gaussian process Thompson sampling (GP-TS), we demonstrate remarkable improvement in both inner- and outer-loop optimization, surprisingly outperforming alternatives like EI and GP-UCB in most cases. Our approach also improves the performance of other posterior sample-based acquisition functions, such as variants of entropy search. Furthermore, we 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.

Published as a conference paper at ICLR 2025
OPTIMIZING POSTERIOR SAMPLES FOR BAYESIAN OPTIMIZATION VIA ROOTFINDING
Taiwo A. Adebiyi University of Houston taadebi2e8 20uh. edu
Bach Do
Ruda Zhang
University of Houston bdo3@uh.
University of Houston
rudaz@uh. edu
ABSTRACT
Bayesian optimization devolves the global optimization of a costly objective func- tion to the global optimization of a sequence of acquisition functions. This inner- loop optimization can be catastrophically difficult if it involves posterior sample paths, especially in higher dimensions. We introduce an efficient global opti- mization strategy for posterior samples based on global rootfinding. It provides gradient-basee optimizers with two sets of judiciously selected starting points, de- signed to combine exploration and exploitation. The number of starting points can be kept small without sacrificing optimization quality. Remarkably, even with just one point from each set, the global optimum is discovered most of the time. The algorithm scales practically linearly to high dimensions, breaking the curse of dimensionality. For Gaussian process Thompson sampling (GP-TS), we demonstrate remarkable improvement in both inner- and outer-loop optimization, surprisingly outperforming alternatives like El and GP-UCB in most cases. Our approach also improves the performance of other posterior sample-based acqui- sition functions, such as variants of entropy search. Furthermore, we propose a sample-average formulation of GP-TS, which has a parameter to explicitly con- trol exploitation and can be computed the cost of one posterior sample.
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 ma-
chine learning models to scientific discovery and engineering design (Jones et al. 1998] Snock
ctt al. 2012 Frazier 2018 Garnett 2023). Many BO strategies are also backed by strong theoretical
guarantees their convergence the global optimum (Srinivas al. 2010 Bull 2011 Russo &
Van Roy] 2014 Chowdhury Gopalan 2017).
Consider the global optimization problem minxex f(x) where X' Rd represents the vector of input variables and [xx E IR the objective function which can be evaluated at significant cost, subject to observation noise. At its core, BO is sequential optimization algorithm that uses probabilistic model of the objective function to guide its evaluation decisions. Starting with prior probabilistic model and some initial data, BO derives an acquisition function axx) from the posterior model, which is much easier to evaluate than the objective function and often has casy-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 func- tion [(x), and the inner-loop seeks to optimize the acquisition function ox(x) at each BO iteration. be solved accurately and efficiently, so
The premise of BO is that the inner-loop optimization 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 (Wilson al. 2018).
Acquisition functions are, in general, highly non-convex and have many local optima. In addition, many common acquisition functions are mostly flat surfaces with few peaks (Rana et al. 2017)

---

Published as a conference paper at ICLR 2025
which take up overwhelmingly large portion of the domain 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) (Chowdhury & Gopalan 2017) uses posterior sample paths directly as random acquisition functions. In many information-theoreeee acquisition functions such as entropy search (ES) (Hennig & Schuler (012).
predictive entropy
search (PES) (Hernández-Lobato
al., 2014), max-value entropy search (MES) (Wang & Jegelkal
2017), and joint entropy search (JES) Tu et al., 2022] Hvarfner et al. 2022). multiple posterior
samples are drawn and optimized find their global optimum location and/or value. Such acqui-
sition functions are celebrated for their nice properties in BO: TS has strong theoretical guarantees
(Russo & Van Roy 2014, 2016) and can be scaled to high dimensions (Mutny & Krause 2018); information-theoretie acquisition functions are grounded in principles for optimal experimental de- sign (MacKay 2003); and both types can be easily parallelized in synchronous batches (Shah & Ghahramani] 2015 Hernández-Lbbal et al., 2017 Kandasamy et al., 2018). However, posterior sample paths 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) (Jones et al. 1998), probability of improvement (PI) (Kushner 1964). and upper confidence bound (GP-UCB) (Srini- vas et al.] 2010). As a consequence, posterior sample paths 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 (Jones al. 1993). covariance matrix adaptation evolution strategy (CMA-ES) algorithms (Hansen al. 2003 and genetic algorithms [Mitchell 1998). Gradient-based multistart optimization, on the other hand, often seen as the best practice to reduce the risk of getting trapped local minima (Kim & Choi 202T), and enjoys the efficiency of being embarrassingly parallelizable. For posterior samples, their global optimization may use joint sampling on a finite set of points (Kandasamy et al. (2018), or approximate sampling of function realizations followed by gradient-basee optimization (Hernández- Lobato et al., 2014] Mutny & Krause, 2018). The selection of starting points is crucial for the success of gradient-basee based multistart optimization. This selection can be deterministic (e.g.. grid search), random (Bergstra & Bengio 2012] Balandat et al.) 2020), or adaptive (Feo & Resende 1995).
In this paper, propose an adaptive strategy for selecting starting points for gradient-based mul- tistart optimization of posterior samples. This algorithm builds on the decomposition of poste- rior 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 paths. The starting points are dependent on the posterior sample, so that each is close to local optimum that is a candidate for the global optimum. The selection algorithm scales linearly to high dimensions.
We give empirical results across diverse set of problems with input dimensions ranging from to 16, establishing the effectiveness of our optimization strategy. Although our algorithm is pro- posed 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 cetiioitttelllliiiiinno..oo balance (Chapelle Li 2011), and can be generated at the same cost one posterior sample.
2 GENERAL BACKGROUND
Gaussian Processes. Consider an unknown function fine x R, where domain X' C Rdd We can collect noisy observations of the function through the model y = fx'( + + with ~ N,, (0, E). To model the function fune, use Gaussian process (GP) the prior probabilistic model: f ~ π E GP. A GP is a random function ff such that for any finite set of points X {x+] 1 121110 n. E N. the values f,- (f(x)) have multivariate Gaussian distribution
(1. n).
2

---

Published as a conference paper at ICLR 2025
Na (phe Kn,n), with mean My (u(x')) ") and covariance Kn,,. [K(x', , the mean function and n(x, x') is the covariance function.
xx) lien Here, μ((x) is
Decoupled Representation of GP Posteriors. Given dataset D = {(x', y')] 1== the posterior model S/D also a GP. Samples from the posterior have a decoupled representation called pathwise conditioning, originally proposed in (Wilson et al. 2020 2021):
(x) = f(x) + K.,n K., (x) + 2)
2) 1y (y - fn - E), f~T, E (0,
(I)
where f(x) is sample path from the GP prior, (x,n(x)) = x1 is the canonical basis, f, (x(x') (x))))) and is sample of the noise. We may interpret f, + as sample from the prior distribution of the observations y = "= This representation has its roots Matheron's update rule that transforms joint distribution of Gaussian variables into conditional one (see e.g.,
Hoffman & Ribak 1191)). This formula is exact, in that & denotes equality in distribution, and it preserves the differentiability of the prior sample. is also computationally efficient for posterior sampling: the cost is independent of the input dimension d, linear in the data size n at evaluation time, and the weight vector for K can be solved accurately using an iterative algorithm that scales linearly with (Lin et al. 2023
3 GLOBAL OPTIMIZATION OF POSTERIOR SAMPLE PATHS
In this section, we introduce an efficient algorithm for the global optimization of posterior sample paths. For this, we exploit the separability of prior samples and useful properties of posterior samples to judiciously select set of starting points for gradient-based multistart optimizers.
Assumptions. Following Section 2 we impose few common assumptions throughout this paper:
()) the domain is a hypercube: X = II [Z, (2) prior covariance is separable: x(x, x') =
IT (4( (3) prior samples are continuously differentiable: f(x:w) C².. Without loss of generality, we also assume that the prior mean μ(x) = 0: any non-zero mean function can be subtracted from the data by replacing fine with fuue μ.. While additive and multiplicative compositions of univariate kernels can be used in the prior (Duvenaud et al. 2013), 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
We observe that, given the assumptions, a posterior sample in eq. T can be rewritten as:
n f(x:0) f(x:w) + ()x f(x:W) II ((:ww))) b(x:ū) Σ xy)(xxx)) x²).
(2)
Here, the prior sample f(x) is approximated as a separable function determined by the random bits w (see Appendix C. Data adjustment b(x) is sum of the canonical basis with coefficients == (Kn,, + 1)-1)(1(- I, e). Both the data adjustment and posterior sample are determined by the random bits W (w) In the following, we denote the prior and the posterior samples as for and for respectively. Our goal is find the global minimum (x ft) of the posterior sample (x).
The global minimization of 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 for surrogate of the posterior sample fo 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 roors as a global optimization algorithm for GP posterior samples (given the assump- tions) via gradient-based multistart optimization. The starting points include: (1) subset S, of the
3

---

Published as a conference paper at ICLR 2025
Algorithm roooo:::: Optimizing posterior samples via rootfinding
Input: prior samples fxx). E; prior covariances x(x.x') Σ; dataset D; set sizes no: Re: Rx-
1: So f minsort (f(x),
& Smallest n,, local minima of the prior sample (Algorithm Smallest ne values and indices of the posterior sample in S. b Smallest ny values and indices of the posterior mean in X >> Starting points: exploration set and exploitation set Gradient-basse multistart optimization > Global minimum of the posterior sample
2: [e, le] 1 mink (f.)) ne)
3: [f,, [ f- minn((((X),))
4: Se f Sole,ii,, Sx t [
5: [*) <- minimim ze (f 8), S = S.S.S
Output: Thompson sample x"
local minima X of corresponding prior sample; and (2) subset S, of the observed locations X. We call S. the exploration set and S, the exploitation set. Specifically, let
S. arg mfn ma)
(3)
xEX,
be the no smallest local minima the prior sample. The set of starting points, S, is defined as:
= Sx, S. = arg
ne), S = arg mink (aaink(f((6 (x), nx).
(4)
XES,
xEX
Algorithm Nutlines the procedure for TS-roots. Here, ng and n, are imposed to control the cost
of gradient- based multistart optimization, and is set to limit the number of evaluations of fa in the determination of S., Considering the cost difference, we can have no >> ne. We observe that
ne and n, can be set to small values, and n, to a moderate value, without sacrificing the quality of optimization, see Appendix D The TS-roots algorithm scales linearly in see Appendix C.5
3.2 RELATIONS BETWEEN THE LOCAL MINIMA OF PRIOR AND POSTERIOR SAMPLES
Figure shows several posterior samples Sa in one and two dimensions, each marked with its local minima X and global minimum X Here the exploration S. is the local minima X_, of fax and the exploitation set S₂ is the observed locations X. We make the following observations:
(1) The prior sample fa is more complex than the data adjustment à 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, fo is almost everywhere one time less differentiable than b (see e.g., Garnett (2023) Sec. 10.2,
Kanagawa et al. (2018)).
(2) Item (1) implies that when the prior sample fw is added to the smoother landscape of b, each local minimum X_, of fa will be located near local minimum Xx of the posterior sample fo. Away from the observed locations X, each N is closely associated with an X with minimal change in location. In the vicinity of X, an & may have both data point x' and an x. nearby, but because of the smoothness difference of fw and b, in most cases the one closest to X is an X...
(3) It is possible that near X, sharp changes in fa may require sharp adjustments to the data, which may move some X by significant distance, or create new X that are not near any X or any xf.
(4) Searching from x" with good observed values can discover good X in the vicinity of X, which can pick up some local minima not readily discovered by X This is especially true if fw is relatively flat x'.
(5) Since the posterior sample f adapts to the dataset, searches from x will tend to converge to a good x² among all the local optima near x'. 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
By combining subsets of X_,, and X, we can expect that the set of local minima of fo discovered with these starting points includes the global minimum X with high probability with respect to W.
4

---

Published as a conference paper at ICLR 2025
prior sample posterior mean
GP-TS acquisition
local minimum
exploitation set
x critical/ point
global minimum
exploration set
3
-1
-1
-2
-2
-2
-1
1
2
-2
-1
1
2
15
::
GP-TS acquisition
Figure 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 X of the GP-TS acquisition function lies outside the interpolation gion, it is typically identified by starting the gradient-based based multistart optimizer at local minimum of the prior sample. Right column: When X within the interpolation region, can be found by starting the optimizer at either observed location or local minimum of the prior sample.
3.3 A REPRESENTATION OF PRIOR SAMPLE LOCAL MINIMA
For each component function fo (wx wi) of the prior sample fu (x), define h, h( (T.), and h₂ hx,) (x,)) for 2; E (I, F.). We call coordinate E [. of mono type (E))))) 0 and call of mixed type (((E))))(EE) of interior critical points f, such that Ex,, (x,, and f((E.)) 0, [1,-.. . Denote E,o = I; and Eir++1 = z. Partition the set of candidate coordinates E, {G.)) into mono type = and mixed type E(") Proposition gives a representation of the sets of strong local extrema of the prior sample. Its proof and the sizes therein are given in Appendix A
0. Let E, [ be the set
15:
Proposition 1 The set of strong local minima of the prior sample (x) can be written as:
X = 8-- = = (E) < <0,, X = (E E fo (E)
0},
(5)
where tensor grids EO) -II (0, The set X., of strong local maxima of for (x) has an analogous representation, and satisfies 8_, U X__ = = UE where the disjoint union.
Critical Points of Univariate Functions via Global Rootfinding, To compute the set Ê of all critical points of is compute all the roots of the derivative f on the interval [.). Since
5

---

Published as a conference paper at ICLR 2025
f is continuous, this can be solved robustly and efficiently by approximating the function with Chebyshey or Legendre polynomial and solving a structured eigenvalue problem (see e.g., (Trefethen 2219)). The roots algorithm for global rootfinding based polynomial approximation is given as Algorithm
3
Appendix C
3.4 ORDERING OF PRIOR SAMPLE LOCAL MINIMA
While the size of X grows exponentially in domain dimension d, its representation in eq. 5 en- ables an efficient algorithm to compute the best subset S. (eq. (3) without enumerating its elements.
With eq. 5 we see that X consists of all the local minima of fw with negative function values. Consider the case where _ has at least no elements so that in the definition of S. we can replace X with X which the largest elements of (x)] within E()) is easier than finding the smallest negative elements of for (x). Once the former is solved, we solve the latter simply by removing the positive elements. Therefore, we convert the problem of eq. 33 into three steps:
which in turn can be replaced with E(1) As we will show later, the problem of finding
1. S(1) arg
(IJ2 am.), with buffer coefficient 2 1;
2. 5 (x E 5(), (x) < 0), so that 5- S X
3. S. arg mink xe8 (fw no), assuming that | 2 nh.
The last two steps are by enumeration and straightforward. The first step can solved efficiently
using min-heaps, with time complexity that scales linearly in Σ | rather than I[ |=
Appendix B The coefficient is chosen that | 2 no- The case when | no < X= can be handled similarly. If n. 2 |. no subsetting needed. The overall procedure to compute S, is given in Algorithm 4 in Appendix C
SAMPLE-AVERAGE POSTERIOR FUNCTION
We finally propose a sample-average posterior function that explicitly controls the exploration- exploitation balance and, notably, can be generated at the cost of generating one posterior sample. Let (ix) = K., (X) + Σ) - y be the posterior mean function. For noiseless observations with /- = w, we can rewrite the GP posterior function in eq. 2 as [x( = f=(x + 7i(x) + (x), where (x)
Kn,n (x) Define a (x) ΣS (x( the sample-average poste- rior function, where (x) are samples generated from the GP posterior and N. E N Since ī(x) is deterministic, and the scaled prior sample VN; fii can be written as fl) iid GP(O, N. xx x'),, we have a,TS (x) = μ(x) + VN. (fw + (x)), where the first and second terms favor exploitation and exploration, respectively. Thus, we can consider N, as exploration exploitation control parameter that, at large values, prioritizes exploitation by concentrating the conditional distribution of the global minimum location, i.c., p2x) at the minimum location of (x),, see Figure 13 in Appendix With aars (x), we can reproduce fa and the GP mean function xx() by setting N 1 and N. oc, respectively.
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 mod- els of random Fourier features (Rahimi & Rechi 2007). This method, unfortunately, subject to the variance starvation problem (Mutny & Krause 2018 Wilson et al., (0200 which can be miti- gated using more accurate feature representations (see e.g., Hensman et al. (2018); Solin & Särkkä 2020)). An alternative is pathwise conditioning (Wilson et al. 2020) that draws GP posterior sam- ples by updating the corresponding prior samples. The decoupled representation the pathwise conditioning can be further reformulated two stochastic optimization problems for the posterior

---

Published as a conference paper at ICLR 2025
mean and uncertainty reduction term, which are then efficiently solved using stochastic gradient descent (Lin et al. 2023).
Optimization of Acquisition Functions. While their global optima guarantee the Bayes' decision rule, BO acquisition functions are highly non-convex and difficult to optimize (Wilson et al. 2018 Nevertheless, less attention has been given to the development of robust algorithms for optimizing these acquisition functions. For this inner-loop optimization, gradient- gradient-based optimizers are often selected because of their fast convergence and robust performance (Daulton et al. 2020). The im- plementation of such optimizers is facilitated by Monte Carlo (MC) acquisition functions whose derivatives are easy to evaluate (Wilson et al. 2018). Gradient-based optimizerrs also allow multistart settings that use a set of starting points which can be, for example, midpoints of data points (Jones) 2001,, uniformly distributed samples over the input variable space (Frazier 2018 Ament et al., 2023), or random points from a Latin hypercube design (Wang et al., 2020). However, multistart- based methods with random search may have difficulty determining the non- flat regions of acquisi- tion functions, especially in high dimensions (Rana et al.) a (017). The log reformulation approach is a good solution to the numerical pathology of flat acquisition surface over large regions of the input variable space (Ament et al. 2023). While this approach works for acquisition functions prone to the flat surface issue such as the family of El-based acquisition functions, its performance has yet to be evaluated for acquisition functions with many local minima like those based on posterior samples.
Posterior Sample- Based Acquisition Functions. As discussed in Section I the family of posterior sample-based acquisition functions is determined from samples of the posterior. GP-TS (Chowd-) hury & Gopalan 2017) is a notable member that extends the classical TS for finite-armed bandits to continuous settings of BO (see algorithms in Appendix E.. GP-TS prefers exploration by the mech- anism that iteratively samples a function from the GP posterior of the objective function, optimizes this function, and selects the resulting solution 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 (Wilson et al., 2018 Balandat et al., (020). Different types of MC acquisition functions can also be used to inject beliefs about functions into the prior (Hvarfner et al. 2024).
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 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 H Our implementation of the TS-roots algorithm is available ttthtps: /gtp//gttthubh//t com/ JQUH/ TSRoots
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 (Surjanovic & Bingham (2013). The analytical expressions for these functions and their global minimum are given in Appendix F
In each optimization iteration, we record the best observed value of the error log(Umin - -) and the distance log (I|xmin - I)), where Umins Xmin- f", and x* 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- (TS-DSRF), GP-TS with random Fourier features (TS-RF), expected improvement (EI) (Jones et al., 1998), lower confidence bound (LCB) the version of GP-UCB Srinivas et al. 2010) for minimization problems, and LogEl (Ament et al. 2023).
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 IO of Appendix[] |I With a fair comparison of outer-loop results (detailed in Appendix H). TS-roots sur- prisingly 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
7

---

Published as a conference paper at ICLR 2025
TS roots ET TS-DSAP TCB
4
TS-maut EI
TS roots - Et TS-DSRP 1CB TS-8F
DS-DSR LCB
log loglghan MM)
Ing
log Shat
3
-2
200
50
100
150
200
600
Iteration
Iteration
Iteration
(b)
(c)
2.5
TS-sola
DS-D
LCB
E1
TS-
T8-06RF
10B
-2
TS-RP
-4
7.7
2-
-6-
1.5
pooba
18 DSIU
LCB
-10
TS-R5
800
200
600
800
100
200
300
350
Iteeattio
Iteration
Iteaaion
(d)
(e)
(f)
Figure 2: Outer-loop optimization results for the (a) 2D Schwefel function, (b) 4D Rosenbrock function, (c) 10D Levy function, (d) 16D Ackley function, (e) 16D Powell function, and (f) ten-bar truss problem. 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.
recommends better solutions than its counterparts, TS-aartt,TSSSSSSSS TS-DSRF and TS-RF, in high-dimensional prob- lems and offers competitive performance in low-dimensional problems. Across all the examples, El 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. TS-roots outperforms LogEl on the 2D Schwefel, 4D Rosenbrock, 10D Levy, and 16D Ackley functions (detailed in Figure of Appendix D.
Optimizing Real-world Problem. We implement TS-roots to optimize an engineered ten-bar truss structure (see Appendix G Ten design variables of the truss are the cross- cross-sectional areas of the truss members. The objective is to minimize a weighted sum of the scaled total cross-sectional area and the scaled vertical displacement a node of interest.
Figures ( and f) show the outer-loop optimization results for the truss obtained from 20 runs of each BO method, where f* is lower bound of the best objective function value we observed from all runs. TS-roots provides the best optimization result with rapid convergence.
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, referred as rootfinding hereafter. We do so by computing the opti- mized values ag of the GP-TS acquisition functions, the corresponding solution points xi. and the CPU times 4 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 x). of the GP-TS acquisition functions by starting the gradient-based optimizer at a large number of initial points (set as 10'), which is much larger than the maximum number of starting points set for TS-roots. For comparison, we extend the same GP-TS acquisi- tion functions to inner-loop optimization using gradient-based multistart optimizer with random starting points (i.c., 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 ge- netic algorithm are equal to the number of starting points recommended for rootfinding. The same termination conditions used for the three algorithms.
Figure 3 shows the comparative performance of the inner-loop optimization for dow-dimensional cases: the 2D Schwefel and 4D Rosenbrock functions. We that the optimized acquisition func-
8

---

Published as a conference paper at ICLR 2025
Rootflading
-100
Random maltistart algorithm
8
-200
-5
Random multistart
Rootrinding
0.6
.300
Genetic algorishm
Rootfinding
Al11
-400
Random multistatt
No.2
Genetic algoritham
-20
0
25
50
75
100
25
50
76
100
0
25
50
75
100
(8)
Iteration
Iteration
Iteration
Rootfiading
Rout
-20
Random moltistart
t4 100s))
0.6
Random mullistart
Gonotie algorithm
Genetic algorithon
40
0.4
-10
-60
log
-15
Random multistart
RoatEndinng
0.2
-80
-20
Genetic algoritiam
0
25
50
75
100
100
25
50
75
100
(b)
Iteration
Iteration
Iteration
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 oz cumulative distances between new solution points x and the true global minima x of the acquisition functions, and cumulative CPU times l for optimizing the acquisition functions.
-0.5
1001
-1
Rootdidding
RostSinding
RootEinding
-1.5
-andom muiltistart
Randorn multistart
Random multistart
algurithm
Genetic algoritlu
Genetie aigorishm
75
100
25
100
25
100
lteration
Iteration
Iteration
10
Rootfiading
15
Rootlinding
10
Rootfuniing
(x1's)
7.5
Random moltistart
Geuelio hum
d6 10%)
10
Randon malsistart
algorithm
100%
7.5
Random muiltattrt
Genetic algorithm
[
2.5
25
50
75
100
100
25
50
75
100
(a)
Iteration
(b)
Iteration
(c)
Iteration
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 acqui- sition functions at and cumulative CPU times t, for optimizing the acquisition functions.
tion values and the optimization runtimes by rootfinding and the random multistart algorithm are almost identical, both of which 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 higher-dimensional settings of the 10D Levy, 16D Ackley, and 16D Powell functions shown in Figure rootfinding performs much better than the random multistart and ge- netic 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 performance of in higher dimensions. Appendix ! provides additional results
random multistart becomes for gradient- based multistart optimization using two other initialization schemes: uniform grid and Latin hypercube sampling. Rootfinding outperforms both, especially in higher dimensions.

---

Published as a conference paper at ICLR 2025
12
MES-R
MES-R 10
MES-R
MES-R
0.5
MES-I 50
log Shatn -['
MES- 50
10
MES TS roots 10
MES-TS-rooie 10
MES- BES-TS-rools 10
log(stin
MES-TS-roota 50
logfgein
MPS. ro-tott 50
00.5
6
50
100
150
200
50
100
150
200
100
150
200
(8)
Iterution
(b)
Iteration
(c)
Iteration
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 [D) The plots are histories of medians and interquartile ranges of solutions from ten runs of each method.
TS-roots to Information- Theorretic Acquisition Functions. We show how TS-roots can enhance the performance of MES (Wang & Jegelka 2017). which uses information about the maximum func- tion value f* for conducting 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* using gradient-based multistart optimizer with large number of random starting points (Wang & Jegelka 2017). We hypothesize that high-quality f* 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 result- ing optimal solutions. Note that the inner-loop optimization of MES, which strongly influences the optimization results, is not addressed by TS-roots.
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-R (Wang& & Jegelka| 2017) and MES-TS-roots correspond to TS-RF and TS-roots, respec- tively, while 10 and 50 represent the number of random samples f* generated for computing the MES acquisition function in each iteration.
10, and MES-TS-roow 50. Here,
Figure 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* are generated. For the 10D Levy function, TS-roots outperforms TS-RF when using 10 samples of f*, while their performance is comparable when 50 samples are used.
Performance of Sample-Average Posterior Functions. We investigate how (aTs(x) influences the outer-loop optimization results. For this, we set No E (1. 10, 50, 100} for TS-roots to optimize the 2D Schwefel, 4D Rosenbrock, and 6D Ackley functions. We observe that increasing N: from to 10 improves TS-roots performance the 2D Schwefel, 4D Rosenbrock, and 6D Ackley functions (see Figure 14 in Appendix M. However, further increases in N. from 10 to 50 and 100 result in slight declines in solution quality TS-roots transitions exploitation. These observations indicate that there is an optimal value of Ne for each problem at which TS-roots achieves its best performance by balancing exploitation and exploration. However, identifying the optimal value to maximize the performance of aars (x) for particular optimization problem remains an open issue.
7 CONCLUSION AND FUTURE WORK
We presented TS-roots, global optimization strategy for posterior sample paths. It features an adaptive selection of starting points for gradient-based multistart optimizers, combining exploration and exploitation. This strategy breaks the curse of dimensionality by exploiting the separability of Gaussian process priors. Compared with random multistart and genetic algorithm, TS-roots con- sistently yields higher-quality solutions in optimizing posterior sample paths, across a range of input dimensions. It also improves the outer-loop optimization performance of GP-TS and information- theoretic acquisition functions such as MES for Bayesian optimization. For future work, we aim to extend TS TS-roots to other spectral representations per Bochner's theorem (Mutny & Krause 2018 Hensman et al. 2018 Solin & Särkkä (020). We also plan to study the ways and the probability of TS- failing to find the global optimum, as well as the impact of subset sizes.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
We would like thank Tian PAN for his help with the maxk_sum problem, and the organizers and the participants of the NeurIPS 2024 Workshop on Bayesian Decision-making and Uncertainty for their invaluable discussions and feedback on an earlier version of the paper. The authors are supported by the University of Houston through the SEED program no. 000189865

---

Published as a conference paper at ICLR 2025
REFERENCES
Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bak- shy. Unexpected improvements expected improvement for Bayesian optimization. In Advances in Neural Information Processing Systems. volume 36, pp. 20577-20612, 2023. URL |nttps://procedddngg fi File/419f  2l1//A87720066622111/23332222222 Confereeeeeeeeeeeeeecccc pdf
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wil- son, and Eytan Bakshy. BoTorch: A framework for efficient Monte-Carlo Bayesian opti- mization. In Advances in Neural Information Processing Systems, volume 33, pp. 21524 21538, 2020. URL http:////proceeddnnn.eeeerrrrr.s fallee/fblli 1W200F111/2511111882776788888111
ftteer/ppeeeeepeee//reea Pappedebeaaaaper..... Paper.
Zachary Battles and Lloyd N. Trefethen. An extension of MATLAB to continuous functions and operators. SIAM Journal on Scientific Computing, 25 2004. doi: 10.1137/ $1064827503450126.
James Bergstra and Yoshua Bengio. Random search for hyper-parameter optimization. Journal of Machine Learning Research, 13(10):281-305, 2012. URL http://jmmmt hwt.ii.errgppppss//////
v1 v13/be 113/erggstrala.hhee
Adam D. Bull. Convergence rates of efficient global optimization algorithms. Journal of Ma-
chine Learning Research, 12(88):2879-290) 2011. URL http: wwtp:i/or/...pprr..///m
#12/bulllia.. html
Olivier Chapelle and Lihong Li. An empirical evaluation of Thompson sampling. In
Advances in Neural Information Processing Systems, volume 24, pp. 2249-2257, 2011. URL https: /pwpess//ppppssppeeee /tttts////appprs.nnnnn http://pppppsssssccppppppaapeeep fww.e/iieeeeeeee www.le/ppapee//11//eeeeeee e
Sayak Ray Chowdhury and Aditya Gopalan. On kernelized multi-armed bandits. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 844-853, 06-11 11 Aug 2017.
mlr wwwppppwsssss.ep 7ttpa////70///hoohrr112222 html
Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hyper- volume improvement for parallel multi-objective bayesian optimization. In Advances in Neural Information Processing Systems, volume 33, pp. 9851-9864, 2020. URL file/
pdf
David Duvenaud, James Lloyd, Roger Grosse, Joshua Tenenbaum, and Ghahramani Zoubin. Struc- ture discovery in nonparametric regression through compositional kernel search. In Proceedings of the 30th International Conference on Machine Learning, volume 28, pp. 1166-1174, 2013. URL http::/
mml...r.....mmmmmm... www.srpp.ssssee.e.e..e...ee.eeee v28 dwv./auddn.aa....a .tmll
Thomas A. Feo and Mauricio G. C. Resende, Greedy randomized adaptive search procedures. Journal of Global Optimization, 602):109-133, 1995. doi: 10. 10.1007//001099663 htthttps:
//doi.org/Di 10. 1100/000100000000
Peter I. Frazier. Bayesian optimization. In Recent Advances Optimization and Modeling of Contemporary Problems, INFORMS TutORials in Operations Research, chapter 11, pp. 255- 278. October 2018 doi: 10. 10.1287/educ.201880188.
Roman Garnett. Bayesian Optimization. Cambridge University Press, Cambridge, UK, 2023. doi: 10.1017/9781108344889
I.S. Gradshteyn and I.M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, Boston, 8th edition, 2014. ISBN 978-0-12-3849333 doi: 10. 10.1066//211--0-448333
Nikolaus Hansen, Sibylle D. Müller, and Petros Koumoutsakos. Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES). Evolutionary Computation, 11(1)11-18, 2003. doi: 10.1162/066655633222828833
12

---

Published as a conference paper at ICLR 2025
Philipp Hennig and Christian J. Schuler. Entropy search for information-efficien global opti- mization. Journal of Machine Learning Research, 13(6):1809-115 2012. URL https:
www. jmir. wwtppppppp..vv....pp.pp haappe/s/vv33/emmllllaa.aaaaa
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian pro- cesses. Journal of Machine Learning Research, 18(151):1-52, 2018. URL http: jmlr..
prg/papers/vee v/8 579.9916669910 html
José Miguel Hernández-Lobato, James Requeima, Edward O. Pyzer-Knapp, and Alán Aspuru- Guzik. Parallel and distributed Thompson sampling for large-scale accelerated exploration of chemical space. In Proceedings of the 34th International Conference on Machine Learning, vol- ume 70, pp. 1470-1479, 1479, 06-11 Aug 2017. URL httpp:://pxoceedii iar.aiinnnsmm.m...m.mmmmmm.im.im press// v70/hernandez www.toll...aaooccooo.o.......... html
José Miguel Hermnddez-Lobato, Matthew W. Hoffman, and Zoubin Ghahramani. Predic- tive entropy search for efficient global optimization of black-box functions.
In Ad- vances in Neural Information Processing Systems, volume 27, pp. 918-920 2014. URL
06 (09530002cc3100055111f1600--00-2
Yehuda Hoffman and Erez Ribak. Constrained Realizations of Gaussian Fields: A Simple Algo-
rithm. Astrophysical Journal Letters, 380:15-L8, October 1991. doi: 10.1086/186160. URL
atttps://ule adsabs www.aarrddeooo. edu/ abs wttp:/a.laaaaa/aa/////// 19911pJ....................._.. 3BOL.
Carl Hvarfner, Frank Hutter, and Luigi Nardi. optimization.
Joint entropy search for In Advances in Neural Infor-
maximally-informed Bayesian mation Processing Systems, volume 35. pp. 11494-11506,
2022.
URL
Carl Hvarfner, Frank Hutter, and Luigi Nardi. A general framework for user-guided bayesian opti- mization. In The Twelfth International Conference on Learning Representations, pp. 9851-9864, 2024.
D. R. Jones, C. D. Perttunen, and B. E. Stuckman. Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Applicaations, 79(1):17-1111, BF00941892. BtLhttts:://ee
79pplcattion,790))5557-----3 1993. doi: 10.1007/
org/ 10 100.0100000011155555500
Donald R. Jones. A taxonomy of global optimization methods based response surfaces. Journal of Global Optimization, 21(4):345-355 -333, 2001. doi: 10. Nt.1022//:1112271055557
Donald R. Jones, Matthias Schonlau, and William Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4))455 492, 1998. doi:
10.10223/A110083364411444
URRp:///:::..../:::::://// ULLit::ss:/aol.org//D00 wttp://///....00011000000000.000
Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaussian processes and kernel methods: A review connections and equivalences, 2018. URL https:
URL /aaxxvv.rro/ab//880000 02582
Kirthevasan Kandasamy, Akshay Krishnamurthy, Jeff Schneider, and Barnabas Poczos. Parallelised
Bayesian optimisation via Thompson sampling. In Proceedings of the Twenty-Firss International Conference on Artificial Intelligence and Statistics, volume 84, pp. 133-142 142, 09-11 Apr 2018. html
Jungtaek Kim and Seungjin Choi. On local optimizers of acquisition functions in Bayesian opti- mization. In Machine Learning and Knowledge Discovery Databases, pp. 675-690, 2021. doi:
1111007/978-33000-66666111
Harold J. Kushner. new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise. Journal Basic Engineering, 86(1):97-106, 1964. doi:
13

---

Published as a conference paper at ICLR 2025
Jihao Andreas Lin, Javier Antorán, Shreyas Padhy, David Janz, José Miguel Hernández-Lobato, and Alexander Terenin. Sampling from Gaussian process posteriors using stochastic gradi- ent descent. In Advances in Neural Information Processing Systems, volume 36, pp. 36886- 36912. 2023. URL httpes///prrreeddnnnnn.aerraaaap 2023/F1 1tt3://66722222222a321oo22200111
David J. C. MacKay. Information theory, inference, and learning algorithms. Cambridge University Press, Cambridge, UK, 2003. URL httts:://wwet www.riiamaeieee http://..dg..oo//778001111233332
Melanie Mitchell. An introduction to genetic algorithms. The MIT Press, Massachusetts, USA, 1998.
Mojmir Mutny and Andreas Krause. Efficient high dimensional Bayesian optimization with addi- tivity and quadrature Fourier features. In Advances in Neural Information Processing Systems, volume 31, pp. 9005-9016, 2018. URL |ttp:t: .ew.:o.eeeeiinnn............p..s mwwrips wwtperissscccppppeee cc/paper/
html
Art B. Owen. Central Limit Theorem for Latin Hypercube Sampling Journal the cul Society; Series B
the Royal Statisti- 1992 doi: 101111110110..111111117716666111 tb01895.x. URL ottpa///////o/oo..oo........///0
Victor Picheny, Tobias Wagner, and David Ginsbourger. A benchmark of kriging- based infill criteria for noisy optimization. Structural and multidisciplinarry optimization, 48:607-626, 2013, doi: 10.
1007/s001188-013-099944 URL wtttt::///// stteps://dow //to://doo.oo. ot1/10.rg11... 1007 000.0000000000001000000-39333
Ali Rahimi and Benjamin Recht.
Random features for large-scale kerell machines. In Advances in Neural Information Processing Systems, volume 20, pp. 1177-1184, cwtpccpr_aaiciceeeeeeee fir saarrffles/ppaee /t11es/ppper/20077
2007. URL ELL2:03800000dd3322fffeeeaaaaaa3 Pap.s5---pppppppp..
Santu Rana, Cheng Li, Sunil Gupta, Vu Nguyen, and Svetha Venkatesh. Bayesian optimization with elastic Gaussian process. In Proceedings of the 34th Interna-
High dimensional
tional Conference on Machine Learning, volume 70, pp. 2883-2891, 06-11 Aug 2017. URL
https:: htttps://proceedings mlr ptt.ppres//70/raaaa
Carl Edward Rasmussen and Christopher K Williams. Gaussian processes for machine learning- The MIT Press, Massachusetts, USA, 2006. ISBN 9780521872508. doi: 10.755 / 001.0001. URL htttp::/ Tttt:s:://ooi.or/1/// wwog/oo/1..0. 75t1/m/tttemmes///2
001.00011
Mark Richardson. chebpy, a Python implementation of chebfun, 2016. URL https: com/ chebpy/chebpy
Daniel Russo and Benjamin Van Roy. Learning optimize via posterior sampling. Mathematics of Operations Research, 39(4):1211-1243, April 2014. doi: 10.1287//oor.2014.00555
Daniel Russo and Benjamin Van Roy. An information-tteeortic analysis of Thompson sampling.
The Journal of Machine Learning Research, 17(1):242-24711, 2016. URL hhtpps: / www.
jmir. wwwpo//ppppsss/////o.p
000.-087.hht 087.html
Amar Shah and Zoubin Ghahramani. Parallel predictive entropy search for batch global optimiza- tion of expensive objective functions. In Advances in Neural Information Processing Systems, volume 28, 2015. URL BRLLAtttss//ppccceiinn..eeeeeeee
httpp//scccppppe____eeeeee files/paper/ Patpeaaald-aaee.ppooom
Jasper Snoek, Hugo Larochelle, and Ryan Adams. Practical Bayesian optimization of ma- chine learning algorithms. In Advances in Neural Information Processing Systems, volume 25, pp. -951, 2012. URL hater//ee25S5111bb11aaaaaaa3aaa3
Amo Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. 30(p))) 30(2):419
Statistics and Computing,
2020. doi: 10. 1007/s 1222-019-09886-1
14

---

Published as a conference paper at ICLR 2025
Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In Proceedings of the 27th International Conference on Machine Learning, volume 13, pp. 1015-1022, 2010. URL pdf
https: / cmml...... httpcccccccccccccccc sece//oonfreennsee/////////eeeee
S Surjanovic and D Bingham. Virtual library of simulation expermmnns: Test functions and datasets, 2013. URLpttp: www. ca/ weaa/ssuuraaaaa https/ssaaaan/pttmmmaaaaaaaaaaaa html
Lloyd N. Trefethen. Approximation Theory and Approximation Practice, Extended Edition. So- ciety for Industrial and Applied Mathematics, Philadelphia, PA, 2019. ISBN 978-1-61197- 593-2. doi: 11113//197811119959999 Seffethen/ATA//
URL /ttpss://peolllle www.ll.........mmm ac.uk/
Ben Tu, Axel Gandy, Nikolas Kantas, and Behrang Shafei. Joint entropy search for multi-objective bayesian optimization. In Advances in
Neural In-
formation Processing Systems, volume 35, pp. 9922-9938, 2022.
URL
https://prccecinsss..iii..o.... hwt.::/iiincccccii........ooocaa f.les//aaer/222/111. 605E558c3887666680055554--2---=3
Jialei Wang, Scott C. Clark, Eric Liu, Peter Frazier. Parallel Bayesian global optimization of expensive functions. Operations Research, 6866)1850-1865 2020, doi: 10. 20.1229/
1966. URL https:: //cbss://do... wtt.aa.rrr.....00
Zi Wang and Stefanie Jegelka. Max-value entropy search for efficient Bayesian optimization. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 3627- 3635, 2017.
mlr wtt.s/mww........eemeemmmmm /
James Wilson, Frank Hutter, and Mare Deisenroth. Maximizing acquisition functions for Bayesian optimization. Advances Neural Information Processing Systems, volume 31, pp. 9884-9895, 2018. URL /t///hccppp//pprcceenng....oo
html
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisenroth.
Efficiently sampling functions from Gaussian process posteriors. In Proceedings of the 37th International Conference on Machine Learning, volume 119, pp. (0292-10302, 13-18 Jul 2020.
URL httphtttss//pprceeeeeet wtt.:nigs........ssssslosc. www.ppresssvlll www.wwlsonnooomm 20 html
James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of gaussian processes. Journal of Machine Learning Research, 22(105):1 2021. URL jmlr http://///....a//aaap///////////
html
Huaiyu Zhu, C. K. Williams, Rohwer, and M Morciniee. Gaussian regression and optimal finite dimensional linear models. In Neural Networks and Machine Learning, 1998. URL https:
//publication aston. ac, wwt..k///c//prirt///////3
15

---

Published as a conference paper at ICLR 2025
A CHARACTERIZING THE LOCAL MINIMA OF A SEPARABLE FUNCTION
A.1 PROOF OF PROPOSITION T A REPRESENTATION OF THE SET OF LOCAL MINIMA
Proposition I brooady applies to separable functions on a hypercube. Consider separable function f(x) I1:-1 fi(2rz) with domain X' -I=_ where fi E C) " [R). To simplify the
discussion, we further assume that fa is twice differentiable at its interior critical points E.. The gradient of f can be written as:
Vf(x) (f((at) I :m1 -
- (f(x)
f(xa) d f((2)) i=1
f(x) v(x).
(6)
where v(x) = (1//A)) = =( log fi) im1- The Hessian of f can be written as: (x;)
f(((x) = diag { (x;))
II [1tTt) }:m1 :=1 + [(()))(()) II (M(w))) iti
f(x) (
ied
3
kti
where s(x) (
(á log fi) fi)
An interior point x E int X' = I._ (c,, Fi) is a strong local minimum of f if and only if Vf(()) = O and V fxx >
0. From eq. (6). the first condition is satisfied in any of the following three cases: (1) (x:) + 0 and (x)) 0 for all E {1.... dd; (2) f1(2;) 0 for exactly one i E (1.... d)
and (x,))) = 0; or (3) f(()))) 00 for all i ETC C (1,-..,) where [ 2 2.
- which is positive definite if and only if one of the following holds: (i) [x) > 0 and ((a))) (14)) > 0, for all i E {1,...,,............ .d]: or (ii) f(x) < 0 and f((z))), < 0, for all i E {1,...,d}.
In case (I), the Hessian eq. 7 reduces to (²f(x) = f(x) ddig
In case (2), the Hessian reduces to an all-zero matrix except for the ith diagonal entry: fx)),, (x,) 11 f (xz)). Even if this entry is positive, the Hessian is still positive semi-definite, which means that there is a continuum of weak local minima: {2;} x 11 Fz]. Besides, this case requires f; and f to have an identical root, which an event with probability zero.
In case (3), let 9a (r:) == f((2, + rs) be shifted version of f.. E {1..... .d... Taylor expansion at
r == 0 gives 9: (r:) = 0 + (0) + o(r;) for all i E I and 93 (ry) = g,) O(r,) for all j & I. I[=_1:((r) (r:)) c cllser ri + o(Ilez r:) OQI (gg T3), where 1l g())) Ie (0) * 0. This means that there is a continuum of saddle points: {si]ier x II Tyl-
We have g(r)
For a boundary point x E 0x : x' inn X. we partition the index set {1...., d} into L, R. and I such that v; I, for all i E L, x; H, for all i E R. and H; E (z,, F.) for all E I. Define v = (0,)se) for any subset of the indices. Then x is a strong local minimum of f if and only if the following conditions hold: (a) x is strong local minimum in (zyber X I1 2, (b)
Vzf(x) 0; and (c) Vxf(x) < O.
Condition (a) holds if any only if vx((x) = 0 and Vzf(x) > 0. Based on the previous discussion
on interior local minima, it is equivalent to: (i) fxx) 0 and f((x;)))) (x.) 0, for all i E 1; or
(ii) fx) < 0 and fz())),, (z)) < 0, for all i E I.
From eq.
6 condition
(b) is equivalent to: (i) fxx) 0 and f(()).((x)) > 0, for all L; or (x)
(ii) f(x)
< 0 and (2, < 0, for all i E L.
(c) is equivalent to: () f(x) 0 and -{(x;))((()) 0. for all i E R: or (ii) fxx) < 0 and (2., < 0, for all E R.
Similarly, condition
Summarizing the above discussions, we that there is unified way to identify the set X of all strong local minima of f, which is stated in Proposition | The discussion for the set X of local maxima is the exactly the same, except that the signs are flipped. This also means that X and X form partition of the union ( L ()) of the two tensor grids.
If f, is not twice differentiable at some interior critical point I,, we may replace f. (zx) > 0 with the statement that n; is strong local minimum of fú. and replace (ma) 0 with the statement
16

---

Published as a conference paper at ICLR 2025
that 2; is strong local maximum of fi. 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 | each set of candidate coordinates E; is partitioned into mixed type and mono type:
= { E E fi(St,)),, < 0}, = = {Si, E E fi(5.,)) > 0}.
Another partition of =; is by the sign of the corresponding component function value:
E ( E E f((5,,) < 0},
- = (E., E f((E).) > 0}.
These two partitions create finer partition of E into four subsets:
E, (1) = =, n=
_. --(0) = =, n = (0)
(1) = _ n =
E =+10 =
Denote the sizes of mixed and mono type candidate coordinates as ni [ (1) and n (0) then the sizes of the two tensor grids E(1) and ((0) can be written as:
= =1=!0/,
NO == | II
n (1) N == [
II
(0)
n
Define signed sums as the sums signs function values on the two tensor grids:
5(1)
sign(f(E)),
((0)
sign(f(E)).
(EEO)
(E
We now derive efficient formulas to calculate these signed sums, using 5(1) as an example. Denote each coordinate in E as el Denote each point in E() as ( (c! ) where multi-index (.)
= (Ji) (T()) = 1I:-1111...... , The signed sum S(!) can be written as:
s(1)
Σ Jeno
sign
(
E
sign
II f.(l))))
JEmo)
Σ Jenoy
II
sign sign((((e))))
I
sign (i (s((s
i=1
II
Σ ((f()))
> 0)
Σ (()))) <0)
II
11100 -
A formula for (()) can be derived analogously. Denote set sizes:
" (1) | (2), n E,- (0)
" +(1) = 1=:
+() 1
n
then the signed sums can be calculated as:
5O)
I ( +1) n (1)),
( (1)),
II
(m: n; (0)).
+())
The sizes of negative and positive strong local minima of separable function can be written as:
1 ( (0(1))~(0)) -5(1)),
- = 1 Σ 1(f(E) < <)
(
ル - = Σ ((f(E) > 0) 2 (Nv1)) + 5(0))
1
(
Therefore, the size of the strong local minima of a separable function can be written as:
Ñ - 2 (N(1) + (= S(1) + 5(0)).
(9)
17

---

Published as a conference paper at ICLR 2025
B ORDERING THE LOCAL MINIMA OF A SEPARABLE FUNCTION
B.I FILTERING A TENSOR GRID FOR HIGH ABSOLUTE VALUES OF SEPARABLE FUNCTION
The step in Section 3.4 is equivalent to the following problem: given coordinates Z, (Si,," ... -ii,, , and components values F; = {fi,,'' fi,, E {1, ,d}, of a separable function f(x I find points such that [f(() are the & largest in the tensor grid =1 Z;
Because log [(x)] = log I fx(x:)] = Σe log [x())),, we can solve this problem as follows: define two-dimensionall arrays F [F..... Fal and A log [F), solve S = maxk-suml um k), and return { (L La, from A that gives the largest sums, which described next.
E S). Here the maxk_sum algorithms finds the combinations
B.2 TOP COMBINATIONS WITH THE LARGEST SUMS
Consider this problem: given two-dimensional array A -a ;ad], = [ai,1, di.t, with 0; 2 ... 2 di, E (1. ......d),,.... d). find multi-indices the form [/1 ... Io] such that the sums 81 := > i=1 0i,1; are the & largest among all combinations II (1, ... 4.).
An exhaustive search intractable because the number of all possible combinations grows expo- nentially as [[ 44. Instead, we use a min-heap to efficiently keep track of the top combinations. A min-heap is complete binary tree, where each node is greater than its children. The oper- ations inserting an element and removing the smallest element from a min-heap can be done in logarithmic time. Algorithm gives procedure to solve the above problem using min-heaps.
Algorithm 2 maxk_summ Combinations with the A: largest sums
Input:
two-dimensional array A; number of top combinations k.
1: Make the array nonpositive by replacing a, with a, a,,1 for i = 1,... .d.
2: Create min-heap by adding the elements of aj, each considered a combination of length
index J, key (1,,1-
3: At stage 2, .... d: create
min-heap consisting of length-i combinations by adding each element in a, to each combination in the min-heap at the previous stage: index (', T,].
key Σj=1 "j The size of the min-heap at each stage is capped at k by popping the smallest sum from the min-heap when necessary.
Output: combinations the min-heap stage d.
This algorithm has time complexity Olk log k), where 6 = Σi= 4; < II 4, and space com- plexity (dk). In TS-roots, the cost of maxk_sum small compared with the gradient-based multistart optimization of the posterior sample.
C ALGORITHM FOR TS-ROOTS
C.I SPECTRAL SAMPLING OF SEPARABLE GAUSSIAN PROCESS PRIORS
Per Mercer's theorem probability spaces (see e.g.. Rasmussen & Williams (206)), Sec 4.3), any positive definite covariance function that is essentially bounded with respect to some probability measure μ on compact domain X' has spectral representation x(x, x') Lk=n λRCE )x)r(x) (x)),, where (Ax, (x((x) (x)) is a pair of eigenvalue and eigenfunction of the kemne integral operator. The cor- iid (O 1)
responding GP prior sample can written as fx) Exe Wi (x where wy. ~ independent standard Gaussian random variables. Similar spectral representations exist per Bochner's 's theorem, which may have efficient discretizations (Solin Särkkä 2020 Mutny &
Krause 2018).
18

---

Published as a conference paper at ICLR 2025
Given spectral representations of the univariate covariance functions of separable Gaussian Process prior, we can accurately approximate the prior sample as:
(x) 2 II [inaww) films;w))
Wii (I.)).
(10)
k=0
Here N₁ is selected for each variable such that Ns /r,1 S 7i, where 1/s sufficiently small (see Appendix for the value used in this study). The product of the univariate sample functions forms third-order approximation of the multivariate sample function; their moment functions differ at order four and higher even orders. Using spectral representations of the univariate components as eq. (TO) is much more efficient than directly using a spectral representation of the separable GP prior, because the former uses Σ N₁ univariate terms to approximate I[ N; multivariate terms in the latter. Another benefit of this representation is that the sample functions are separable, which can be exploited in downstream tasks.
Spectrum of the Squared Exponential Covariance Function. The univariate squared exponential (SE) covariance function can be written as ro(z, =²; = exp(-}s²) where the relative distance I //l/ and length scale I (0, oe). The spectral representation of such covariance function per Mercer's theorem is k(x, x²) = Ex= ( (a'). With a Gaussian measure fs -N(0, o²) over the domain X' I. we can write the eigenvalues A and eigenfunctions @& (x) of the kernel integral operator as follows. (See e.g., ZZuu et al. (1998) Sec. 4 and Gradshteyn & Ryzhik (2014) 7.374 eq. 8.)
Va² + 4ab, and A 2aa b f.c. For k N, the kth eigenvalue is A V ( and the corresponding eigenfunction is dx(x))
Define constants
(202)-1, b (21)-1,
Ha exp (- 2²² and (z) the kth-order Hermite polynomial defined by Hx(x) (-1)* exp(x²)
(2e) 1/4 (x(vCr) exp (faz²), where (x) = (1' - 1/2
dxs exp(-2"
d
Figure shows approximations to the SE covariance function by truncated spectral representations with the first N eigenpairs and by random Fourier features (Rahimi & Recht (007) with N basis functions. The spectral representation per Mercer's s theorem converges quickly to the true covariance function, while the random Fourier features representation requires a large number basis functions and is inaccurate for N < 1000.
... True N
True
1000
0.75
N
100
2000
k(2, zz)
N 100
0.5
500
5000
0.25
-0.5
0
1
2
4
5
0
1
2
3
4
(a)
27/
(b)
37
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
C.2 UNIVARIATE GLOBAL ROOTFINDING
Algorithm 3 outlines a method for univariate global rootfinding on an interval by solving an eigen- value problem. When the orthogonal polynomial basis is the Chebyshev polynomials, the corre- sponding comrade matrix is called colleague matrix, and we have the following theorem:
Theorem Let p(z)= En og Ta (x), am + O, be a polynomial of degree m, where T the kth Chebyshev polynomial and dg is the corresponding weight. The roots of p(x:) the eigenvalues of
19

---

Published as a conference paper at ICLR 2025
Algorithm 3 roots: Univariate global rootfinding an interval
Input:
polynomial p(z) of degree (or any real function f(æ))
1: transform p(z) into an orthogonal polynomial basis p(ıx) = Ek (x(x)
(or approximate (ax on the interval using such basis)
2: solve all the eigenvalues of the comrade matrix C associated with the polynomial basis
Output: all the real eigenvalues { in the interval, which are the roots of p(z) (or f()))
the following X colleague matrix:
0
1
1/2
0
1/2
1/2
0
1/2
C
r (1I)
20m
1/2
1/2
(10
01
a2
Om
where the elements not displayed are zero.
proof of Theorem provided in Trefethen (2019). Chapter 18. A classical formula to compute the weights } requires O O(m²) floating point operations, which can be reduced to Om log m) using fast Fourier transform. Since the colleague matrix is tridiagonal except in the final row, the complexity of computing its eigenvalues can be improved from ((m²) to O(m²) operations, which can further improved to ((m) via recursive subdivision of intervals (see Trefethen 2019) )
Specifically, if m >> 100, the interval is divided recursively so that each subinterval the function can accurately approximated by polynomial
degree no greater than 100. The roots s gorithm is implemented in the Chebfun package in MATLAB (Battles & Trefethen (0044) and the chebpy package in Python (Richardson) (216)); both packages also implement other related programs such chebfun for Chebyshev polynomial approximation and di f f for differentiation.
C.3 BEST LOCAL MINIMA OF SEPARABLE FUNCTION
Given the univariate component functions of a separable function, Algorithm finds the subset S. of the local minima the function with the n, smallest function values. This procedure requires the maxk sum algorithm in Algorithm 2 the root algorithm in Algorithm and the related programs chebfun and di f, see also Appendix
In Algorithm E, f. h, J, P are two-dimensional arrays, while (T) (T) are matrices. Func- tion evaluations Lines 10 14
20 and 24 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(g (I) < 0 at Line 9 can be evaluated as rowXor(P(L)(() (T)), where P()) is two-dimensional array with p(1) P-(-), P())(I) is ma- trix with d columns, and rowXor is row-wise exclusive or operation. Similarly, the statement f(g() (T) at Line can be evaluated as rewProd(f!)) (T)), where rowProd is row products.
19
C.4 DECOUPLED SAMPLING FROM GAUSSIAN PROCESS POSTERIORS
The decoupled sampling method for GP posteriors Wilson et al. 2020). together with the spectral sampling of separable GP priors, is outlined in Algorithm
C.5 COMPUTATIONAL COMPLEXITY OF TS-
Per Algorithm M the computational cost of the TS roots method dominated by a few tasks: (1) one call of mi nsort (Algorithm 4 (2) n. + n evaluations of the posterior sample path f(-); and (3) ne + ny calls of the gradient-basee based optimizer minimize.
First, consider task (2). Evaluating f(-) involves evaluating: (i) the prior sample path f(-), which involves evaluating its d univariate component functions, each with cost that depends its spectral
20

---

Published as a conference paper at ICLR 2025
Algorithm 4 minsort: Best local minima of a separable function
Input: separable function f(x [ fi(x;): set size no; buffer coefficient (defaults to 3).
f((z,) < C 1....,......... Construct chebfuns for univariate components
(x;; di.f(f((z))):
(x;) +- (z)f( > Compute first and second derivatives Univariate global rootfinding
2: [Eia)) 1/ f roots (T))),
(E.o) +2 (Etn+1) + Tz
>> Include interval lower and upper bounds
Sur.iintt 3: f + (.(E,), 10, ,,
&, (...
Candidate coordinate values (E.)
D Univariate function values Univariate second derivatives critical points Univariate inward derivatives at interval ends Boolean vectors sign parity and positivity Mono and mixed type candidate coordinates: = =
h & 1. ((n,))
== =1, -,,,,,
14.0 + (((E,));
him+1 t +- -f((E
4: J: (ooh; 0);
P. 2f,
5: 4 + {((Ji); 6
E-(-ii)
> Values mono and mixed type candidate coordinates
(0) f(J:)::
(((-l))
6: " (0) sJm(J
ni + suum(-J:)
(1)
" +(0) + (um(li n + sum(( (a)((-P)))J.))
-())
n +(1) sum (Pia(-J)))
-(1) sum( (PP:)& &-J()))))
n,
N(0) - I[ 10) N + II (1)
Sizes of tensor grids Signed sums
5(O) + I
(n.
(n((0)
+(0)
- n, ) SO)) I (mt))
+(1) - n4 (1) [
7: if n. I Ñ = (NC) s()) then
8:
[s. + maxk_suu (1og((1 (( ())),
> The an. largest |f] in =())
9:
I+I (I) 0,:]
Multi-indices of best negative local minima The n, smallest ff X
10:
- -.]+miinkff(e (L)), ng)
11:
S = (II,:)
else
13:
1T (1. (1)
Matrix of index combinations Multi-indices of negative local minima D Sort values in ascending order Negative local minima
+ 11 (f(
0,:]
15:
- (f (T))
16:
X + (i
if no - = - S(1) (()) 5(0) then
17:
18:
s. + maxk_sum (log( (f. 1) 1) a(n ConnnnNi.. -N-)
D Largest [f] in (
19:
(- 0.;
Multi-indices of best positive local minima The no N smallest in ++
20:
b, < minkk (T), -
21:
S. X-+ S g [I[,:))
22:
else
(0) - (1, U
Matrix of index combinations
24:
i + 17 nd [S (II > 0,:]
Multi-indices of positive local minima
25:
b, < sort(f(g()) (I+))
D Sort values in ascending order All local minima
26:
S. 8 U X+ (0) ([+
27:
end
end if
Output: S.
> Candidate exploration set: smallest local minima in ascending order
representation (Appendix C.1 1); and (ii) the canonical basis (.),,(),, which costs O(dn) flops. When the data size n is large, we can pre-filter the observed locations X by the observations y, which is good estimate of f(X) depending on the observation noise. We assume that the number of observed locations after filtering
most comparable to ne. The cost of task (2) thus ((n,dn) flops.
Now consider task (3). Evaluating the gradient of F(-) involves evaluating the gradients of f()) and K. (.). Since both f(-) and k.,,(.) separable functions, their gradients be computed at
21

---

Published as a conference paper at ICLR 2025
Algorithm 5 Decoupled sampling of Gaussian process posterior
Input: eigenpairs {(Ar,,, ds (2) (z)) ( ,N-1 data D = {(x), 1زز
covariance matrix C
K n,n + Σ, canonical basis (x,,(()
1: Wi,k
iid (N.01)) 1)
D Random coefficients for the prior sample D Approximate prior sample Values of prior sample at observed locations D Random noise for the posterior sample b Linear solve via factorization (c.g., Cholesky or SVD) Approximate posterior sample
2: (x) = -II
Mi
kk= waye A
3: f 6 (x*)
[x(
4: E N,,
5: V C-1 (y (y- c))
Output: fx(x) = (x) + (x)
cost comparable to evaluating their function values. Let Ngra be the number of gradient evaluations required by the gradient-based optimizer. The cost of task (3) is thus O(Ei + Thx + graadd)) flops.
For task (I), the cost of the minsort algorithm is dominated by: (i) d calls to chebfun, which evaluates the univariate components ft(.) at a number of points depending on their complexity; (ii) d calls to roots (Algorithm 3. which scales linearly with the polynomial degree m of the chebfun object, itself dependent on the complexity of f():;; and (iii) at most one call to maxk.sum (
2 at a cost of (tn. log no)) where t = Ei= t, and i; is two plus the number of critical points of fi(-) which depends on the complexity of fi()). The complexity of fi).), for the SE kernel for example (see Appendix [.I).. can be quantified as the inverse length scale 0, = 1/4;- We may define average complexity as 0 2 0.. The cost of task (1) is thus OdOnn, log no) flops.
As explained in Appendix D we can set ne and ni to small values and n, to a moderate value, independent of f(.) and thus independent of d, n and 0. The overall cost of TS-roots thus scales as Odn + d0), which is linear in the input dimension d.
D MINIMUM SIZE OF EXPLORATION AND EXPLOITATION SETS
We conduct an empirical experiment to determine minimal values for ne and n, of TS-roots algo- rithm, which are the sizes of S. and Sx, respectively. From this experiment, we also recommend a value for na, which is the size of Sg- Recall that points in S, are sorted in ascending order of prior sample values, while those in S. and S, are sorted in ascending order of posterior sample values.
Let Ie and I2 be the sets of indices of points in S. and S. that converge to the best local minimum of the posterior sample in each optimization iteration, respectively. Let I, be the set of indices of points in S. associated with S.. Our hypothesis is that we have high chance of finding small index value in either I. or Ix If this hypothesis confirmed, then we can set both nee and ni at very small values, which significantly accelerate the inner-loop optimization. To confirm the following two steps. First, we set ne and n. at large values to mimic the effect of removing the set size limits, ensuring accurate solution of the global optimization problem. Then, we show that we have a high chance of finding a small index value from Te and/or Ig in each optimization iteration.
hypothesis, we employ
We test our hypothesis on the 2D Schwefel, 4D Rosenbrock, 10D Levy, 16D Ackley, 16D Powell functions. We n% 5000, The The 1000, and 3 (buffer coefficient). The left and middle columns of Figure show the smallest index values and the variation of index values from le and/or Is of starting points that converge to the best local minimum x" of the posterior sample path in each optimization iteration. The left column also plots the index values from I, corresponding to the smallest index values from Te, exists. The right column shows the histograms of the smallest index values from Te and 7, for all iterations considered. These results show that we have a high chance of finding a small index value from Le and/or L in each iteration. This confirms our hypothesis. In fact, using the first point in S. and the first point in S.-only two points-ww we can discover the global optimum most of the time. Interestingly, the smallest index values appear largely independent of both the optimization iteration and the input dimension. Furthermore, the results suggest that it is safe to set n, = 500; and for almost exact global optimization, we suggest setting n. = 25 and n, = 50.
22

---

Published as a conference paper at ICLR 2025
From
0.5
From F
Index
150
Inder L and
10
0.6
50
0.2
50
100
150
200
100
150
200
5
10
(a)
Iteration
Iseration
Index
50
0.8
Prom
200
Froo
40
0.6
Prom f.
Index
150
Todos Lange L. and
30
0.1
0.2
50
10
0
50
100
150
200
50
100
150
200
1
10
20
30
40
50
(b)
Tteration
Iteration
Index
500
175
0.8
From
From
Fou
Index
300
Inlec mage of and
0.4
75
50
25
0
1000
2001
300
400
100
200
300
400
50
100
150
(c)
Iteration
Index
1000
200
175
From
Eroe.
Froen
о trom
Indox
600
Exfex Waoge of and
Indox
125
0.6
100
400
0.4
1
100
200
300
400
200
300
400
100
130
200
(d)
Iseratiin
Iteration
Iodex
1000
150
x Fram
From
800
a From
125
0.8
From
o From
Index
600
Taage und
100
75
400
0.4
25
1
200
200
400
150
(c)
Iteration
Iteration
Index
Figure 7: Left column: Minimum index values and index range of I, and/or I, for starting points that converge to the best local minimum x" of posterior sample in each optimization iteration, and index values of Io associated with minimum index values from Je. Middle column: Zoom-in plots of index values. Right column: Historam of the minimum index values. (a) 2D Schwefel, (b) Rosenbrock, (c) 10D Levy, (d) 16D Ackley, (c) 16D Powell functions.
E BAYESIAN OPTIMIZATION VIA THOMPSON SAMPLING
general procedure for sequential optimization is given in Algorithm The initial dataset 7 can either be empty contain some observations. In the latter case can write D = [(x', y 3 =m=, where ne E N Three components of this algorithm can be customized: the observation model Obscrve(x), the optimization policy Policy(D), and the termination condition.
BO can be seen as an optimization policy for sequential optimization. A formal procedure is given in Algorithm Three components of this algorithm be customized: the prior probabilistic model
23

---

Published as a conference paper at ICLR 2025
Algorithm 6 Sequential optimization (Garnett 2023)
Input: initial dataset D
1: & <-1
2: repeat
3:
x <- Policy(Di-1)
4: y f- Observe(x*)
5: D f- Dé-1 U {(x* y
6: until termination condition reached
Output: D
Algorithm 7 Bayesian optimization policy
Input: prior stochastic process f for the objective function fures current dataset D
1: determine the posterior f* : ==]DD-
2: derive an acquisition function a (x) from f*
3: global optimization x* +- arg min a* (x)
Output: x*
f. the acquisition function a, and the global optimization algorithm. Any probabilistic model of the objective function fune can be seen probability distribution on a function space, and the prior is usually specified as a stochastic process such GP. The acquisition function a derived from the posterior f| can be either deterministic such as EI and LCB or stochastic, such GP-TS. To simplify notation, we state the global optimization problem (x)) as minimization rather than maximization. The two problems are the same with change of sign to the objective.
When applied to BO, GP-TS generates random acquisition function simply by sampling the pos- terior model. That is, given the posterior f at the kth BO iteration, the GP-TS acquisition function a random function: a axxx) ~ ft.
F BENCHMARK FUNCTIONS
The analytical expressions for the benchmark functions used in Section are given below. The global solutions of these functions detailed in (Surjanovie Bingham 2013)
Schwefel Function:
f(x) = 418.98294
2 sin
(12)
This function is evaluated X' [-500, 500] and has global minimum f= f(xx) 0 x = (420.9687,.... 420.9687| This function is C at x 0.
Rosenbrock Function:
f(x) Σ 100 - x2)²
(x: 1)2].
(13)
This function evaluated on = [-5, 10)4 and has a global minimum f = 0 at x = ,1,
Levy Function:
f(x) = sin? (Tw)) +
(w, 1)² [1+ s0 (TW, + 1)] + (wa - 1)² [1 + sin (2Twa)) (14)
where wl = 1 + 11 1, 1,... d. This function is evaluated on = [--10, 10]4 and has a global minimum 0 at x* [1, [1,... T]T.
24

---

Published as a conference paper at ICLR 2025
Ackley Function:
f(x) -a exp
b
x
exp
d
Σ cos(cr:)) + a exp(1),
(15)
where a - 20, 0.2, and C 2T. This function is evaluated on = [-110, 10j4 and has a global
minimum f 0t [0.... 0 This function not differentiable at x*.
Powell Function:
f(x) Σ [ 10x44-2)² 2 + (24i-1
-xa)²² + (x4-2 20 10 (I41-3 x4:)
(16)
This function evaluated on X [-4, 5j4 and has a global minimum fr= O 10,
,0,
6d Hartmann Function:
4 Σ
f(x)
a; exp
- A
(17)
1=1
j=1
where
[1,
3,3,3,
(18a)
10
17
3.5
1.7
8
0.05
10
17
0.1
8
14
A
(18b)
3.5
1.7
10
17
8
0.05
10
0.1
14
1312
1696
5569
124
8283
58867
2329
4135
8307
3736
1004
9991
P = 10-1
(18c)
2348
1451
3522
2883
3047
6650
14047
8828
8732
5743
1091
381
This function is evaluated on X' = [0, 1] and has a global minimum f =
3..2237 at
x (0.20169. 0. 150011 0. 0.476874, 0.275332, 0.31 0.6573] The rescaled version f(x)
[(x) (x))-2.5 2.58
Picheny et 2013 is used in the experiments.
1.94
G TEN-BAR TRUSS
Consider ten-bar truss shown in Figure g The truss has ten members and is subjected to vertical load P1 60 kN node 2, vertical load P3
kN node 3, and horizontal load P3 10 kN at node 3. The Young's modulus of the truss material E 200 GPa. The length parameter == 1 Let A(x) Ei I₁ and 6x(x) denote the total area of the cross-sectional areas of the truss members and the vertical displacement at node 3, respectively. where x = [z]..... _1o]T is the vector of cross-sectional areas of the truss members. The optimization problem formulated for the truss is to minimize both A(x) and d3 (x)). Since A(x) and (x))) are competing, we define the objective function weight-sum of A(x) and (x)x), such that
A(x) W1 8x(x)
f(x) W1
(19)
Amax 8max
where x E [1,20] cm?, w1 0.6, w2 = 0.4, Amax 200 cm3 and 6max = 3 cm.
H EXPERIMENTAL DETAILS
Data Generation. We generate 20 initial datasets for each problem. The input observations are randomly generated using the Latin hypercube sampling (Owen 1992) 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.
25

---

Published as a conference paper at ICLR 2025
L
5
X4
6
X9
E
X5
X7 X10
X1
3
Pi
P2
Figure 8: Ten-bar truss. Cross-sectional areas of ten truss members are the input variables Zis i E (1..., 10}. Known parameters include length L, Young's s modulus of truss material E, and external loads Pj.j E (1, 3.3. The vertical displacement node is denoted as 63.
Key Parameters for TS-roots and other BO Methods. We use squared exponential (SE) covari- ance functions for our experiments. The spectra of univariate SE covariance functions for all prob- lems (see Appendix C.I) are determined using the Gaussian measure μ = _N 1). The num- ber of terms N;, E {L. [I..... d]. of each truncated univariate spectrum determined such that A1, T7, where = 10 If N; > 1000, we set N, 1000 to trade off between the accuracy of truncated spectra and computational cost. We also set n, 500. The maximum size of 250. The maximum size of the exploitation set is
the exploration set is
200.
The number of initial observations is 10d for all problems. The standard deviation of observation noise 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) characterized by total of 2000 random Fourier features.
To ensure a fair comparison of outer-optimization results, we first implement ro-ro 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 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 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 designated cluster at our host institution. This cluster hosts 9984 Intel CPU cores and 327680 GPU nodes. The inner-loop optimization implemented on PC with an Intel Core 17-1165G7 @ 2.80 GHz and 16 GB memory.
Nvidia GPU cores integrated within 188 compute and 20
For the univariate global rootfinding via Chebyshev polynomials, we use MATLAB's Chebfun pack- age (Battles Trefethen 2004) and its corresponding implementation Python, called chebpy [Richardson] 2016).
I ADDITIONAL RESULTS
Comparison of TS-roots and LogEI Outer-loop Optimization Results. Figure 9 compares the
performance of the outer-loop optimization of TS-roots and LogEl (Ament al. 2023) the five
benchmark functions and a real-world ten-bar truss structure. Across all examples, save 16D Powell, TS-roots outperforms LogEl, further demonstrating the robustness of approach. The compar-
26

---

Published as a conference paper at ICLR 2025
1--roots Logil
18-oots LogEi
LogE
logiVers -f'
You
logly
& 4
200
600
800
600
800
Iteration
Iteration
Iteration
(a)
(b)
(c)
2.50
roots
TS-not Log
TS-roots LogEi
LogEl
2.25
7.5
2.00
1.75
68.5
loglymn-
$1.50
2600
1:25
5.5
1.00
0.75
4.5
200
800
Iteration
Iteration
Iteration
(d)
(f)
Figure 9: Outer-loop optimization results for the (a) 2D Schwefel function, (b) 4D Rosenbrock
function, (c) 10D Levy function, (d) 16D Ackley function, (e) 16D Powell function, and (f) ten-bar truss problem. The plots are histories of medians and interquartile ranges of solution values from 20 runs of TS -ootot and LogEl.
T9 roots EI
288
1S-ooss Ei
TS-oooss 81
TS-DSRF
LCB
TS-NSRF LCB
LCB
loz (Iee ()
TS-RP
7S-RF
TS-RP
Jog (Dxa
2.1
1.7
1.3
150
100
150
200
200
400
600
800
(u)
Iteration
(b)
Eteration
(c)
Iberation
T8
EI
TS-mota TS-DSRE
TS-DRR
1CB
2.5
I5 JF
log Ken
Neen
7.5
0.5
200
400
600
800
0
200
4011
600
800
(d)
Iteration
(c)
Itesation
Figure 10:
Outer-loop optimization results for (a) the 2D Schwefel, (b) 4D Rosenbrock, (c) 10D Levy, (d) 16D Ackley, (c) 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.
atively weaker performance on the 16D Powell function can be attributed its convex landscape, where BO generally underperforms relative to gradient-based methods.
Distance to Global Minimum. Figure 10 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.
27

---

Published as a conference paper at ICLR 2025
Roatlinding Uniform grid
Rootfinding
Rootšmding
Uniifor grid. -.H design
-Uniform grid -.HS design
-10
LHS design
-10
W
-80
A1-5
10
-20
-120
0
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
(a)
Iteration
Iteration
Iteration
75
Routdinding Uniform grid
Rootfinding
Rootfinding
Unniform grid design
-2
8
LHS design
-5
LHS design
Unilfrmm grid
W1.50
-10
25
-15
-75
-20
0
20
40
50
10
30
50
10
20
30
40
50
(b)
Iteration
Iteration
Iteration
Figure 11: Inner-loop optimization results by three different initialization schemes, i.e., rootfinding, uniform grid, and Latin hypercube sampling, for (a) the 2D Schwefel and (b) 4D Rosenbrock func- tions. The plots are cumulative values of optimized GP-TS acquisition functions o cumulative distances between new solution points x and the true global minima x,, of the acquisition func- tions, and cumulative CPU times & for optimizing the acquisition functions.
xx00
-10
(x11))
-1
-20
Rootflading
Roo-fiuding
-2
RoetEnding
Uniform grid
Uniloom grid
-Uniform grid
LLIS design
....H design
LHS design
-40
10
30
50
Iteration
Ieration
Ieration
RootEnning
Rovelinding
Rootfinding
fx
6
Uniforr grid
(x1's))
3
Uniform grid
-Uniform grid
-.HS design
-HS design
2
-HSS design
2
444
2
0
10
20
30
40
50
(a)
Iteration
(b)
Iteration
(c)
Iteration
Figure 12: Inner-loop optimization results by three different initialization schemes, i.e., rootfinding, uniform grid, and Latin hypercube sampling, for (a) the 10D Levy, (b) 16D Ackley, and (c) 16D Powell functions. The plots are cumulative values of optimized GP-TS acquisition functions of and cumulative CPU times tę for optimizing the acquisition functions.
Comparison of Inner-loop Optimization Results. Figures 11 and 12 compare the performance of the inner-loop optimization by three different initialization schemes, 1.e., rootfinding, uniform grid, and Latin hypercube sampling, for oow-dimessional cases of the 2D Schwefel and 4D Rosenbrock functions, and for higher-dimennional cases of the 10D Levy, 16D Ackley, and I6D Powell func- tions. Rootfinding performs better than the uniform grid and Latin hypercube sampling initialization schemes, especially in nigh-dimensional settings.
Sample-average Posterior Function. Figure B shows how we can improve the exploitation of GP-TS when increasing the exploration-exxloiiaao control parameter Ne.
28

---

Published as a conference paper at ICLR 2025
observations
OATS, No
30
GP mean
OSIS, No
195% Pls
No 50
function value
20
10
0
-10
-20
0
10
15
20
0.1
N
0.05
5
10
15
20
p(3" (D)
40
No
20
0
5
10
15
20
200
N
50
100
0
5
10
15
20
E
Figure 13: Sample-average posterior function for different values of No. The posterior function
approaches the GP mean and the conditional distribution of the solution location p(z* is more concentrated when increase N.
10
2.5
N
logs - f*)
N. 60
50
2
N, 30
2
100
log Tugin
Tog
1.5
100
-2
50
100
150
200
100 Iteration
150
200
100
200
300
Iteration
Iteration
N.
6
N:
Xaiia x
N. 50
log Xain
Xman
log
log
10
100
50
100
150
200
50
100
150
200
100
200
300
(8)
Iteration
(b)
Iteration
(c)
Iteration
Figure 14: Performance sample-average TS-roots with different control values N, 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-roo "S-roots for each N. value.
Performance of Sample-average TS-roots. Figure 14 shows the performance of sample-average TS-roots with different exploration-exploitttio control parameters N. for the 2D Schwefel, 4D Rosenbrock, and 6D Ackley functions.
29


Reviews:
[{"review_id": "UL4GbgPGGw", "reviewer": "Reviewer_MHUs", "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. \nSpecifically, 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.\n\nAnother 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.", "questions": "-\tThere 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?\n\n[1] Ament, Sebastian, et al. \"Unexpected improvements to expected improvement for bayesian optimization.\" Advances in Neural Information Processing Systems 36 (2023): 20577-20612.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "-\tThe paper addresses the problems of finding the global optimum of GP posterior samples efficiently, which is an important topic for many BO acquisition functions.\n-\tThere are many empirical analyses for the proposed methods.\n-\tThe method significantly reduces the CPU times to optimize GP-TS compared against genetic algorithm and random multi-start.\n-\tThe 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.", "weaknesses": "-\tThe 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.\n-\tThe idea seems to not very effective in improving MES acquisition function. Fig. 5 shows that TS-Random Fourier and TS-root are statistically similar.\n-\tNo real-world benchmark problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Zlp7MMaY5b", "reviewer": "Reviewer_D2as", "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.", "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?\n\n2. 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?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 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.", "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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "s838Ro4DyH", "reviewer": "Reviewer_kcrQ", "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.", "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?\n\n2. 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?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This paper selects an uncommon topic of BO, namely global optimization of acquisition functions, and conduct a systematic research on it.\n2. 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.", "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.   \n\n2. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "UPzss31Uyn", "reviewer": "Reviewer_AkKy", "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.", "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. \n\n[1] Entropy Search for Information-Efficient Global Optimization (Hennig and Schuler, 2012)", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The method is motivated well (Section 3.2)\n* The method performs better than (or at least on-par with ) the baselines empirically. \n* 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.\n* Figure 1 is very useful", "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. \n* 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)