Triangulation candidates for Bayesian optimization
Authors: Robert B. Gramacy
Annie Sauer
Nathan Wycoff
Venue: NeurIPS 2022 Accept (2022)
Keywords: surrogate modeling
Gaussian process
active learning
sequential design
space-filling design
Delaunay triangulation
convex hull

Abstract: Bayesian optimization involves "inner optimization" over a new-data acquisition criterion which is non-convex/highly multi-modal, may be non-differentiable, or may otherwise thwart local numerical optimizers.  In such cases it is common to replace continuous search with a discrete one over random candidates.  Here we propose using candidates based on a Delaunay triangulation of the existing input design.  We detail the construction of these "tricands" and demonstrate empirically how they outperform both numerically optimized acquisitions and random candidate-based alternatives, and are well-suited for hybrid schemes, on benchmark synthetic and real simulation experiments.

Triangulation candidates for Bayesian optimization
Robert B. Gramacy Department of Statistics Virginia Tech Blacksburg, VA 24060, USA rbgēvt. rdg
Annie Sauer Department Statistics
Virginia Tech
Blacksburg, VA 24060, USA
anniees@vt.co edu
Nathan Wycoff
The McCourt School's Massive Data Institute
Georgetown University
Washington DC, 20057, USA
nathan trobamnwyyoff@geegggetooommmmmm edu
Abstract
Bayesian optimization involves "inner optimization" over new-data acquisition criterion which is con-convex/highly multi-modal, may be non-differentiable, or may otherwise thwart local numerical optimizers. In such cases it is common to replace continuous search with a discrete one over random candidates. Here we propose using candidates based Delaunay triangulation of the existing input design. We detail the construction of these "tricands" and demonstrate empirically how they outperform both numerically optimized acquisitions and random candidate-basee based alternatives, and are well-suited for hybrid schemes, on benchmark synthetic and real simulation experiments.
Introduction
We address the continuous unconstrained optimization problem
argmin f(x) .EB
(I)
where the bounding box B a hyperrectangle, often taken as [0, 1] in coded inputs. The objective f B + R is a blackbox function, meaning that we can only learn about its behavior through expensive, often simulation-base based, evaluation. Such problems are most challenging when f highly non-convex, and thus contains multiple local minima. A tacit goal of solver is to minimize the number times that f is evaluated in search a global solution.
The earliest papers on Bayesian optimization (BO) adapted statistical modeling and design principles tackle this optimization problem [Močkus] 1975| Jones al. 1998 Applications on physics- based simulators f are provided by Pourmohaaad and Lee [2021]: scenarios in machine learning are reviewed by Garnett [2022]. BO common studies of engagement and user experiences in online platforms le.g.,LLthhan and Bakshy 2019]. hyperparameter estimation for deep learning [e.g. Turner Feurer et 2018 and materials design [e.g. Zhang al. 20206), name a few.
al. 2021
To illustrate BO and introduce our contributions, consider f(æ) = sin(r) and = 1, 2T + 1]
evaluated at six equally-spaced inputs x. Next fit surrogate fn. to data (X,, Yn), where y: f(x;) for 1,..., n We privilege Gaussian process (GP) surrogatess [e.g.. Gramacy 2020
36th Conference on Neural Information Processing Systems (NeurIPS 2022).

---

Expected Improvement
0.20
0.15
00
0.05
0.0
Predictive surface (left) via means (solid black), 90% intervals (dashed-red), and sample paths (gray); HI (right) and resulting acquisition (green diamond).
Figure 1:
but our methodology is agnostic that choice so long the predictive equations from fn have similar features c.g., non-linear predictive mean μr, (x) and higher predictive uncertainty/standard deviation o, (x)) away from the training data sites.
The left panel of Figure shows a fitted fe via ne as solid black line with error bars 146 = 2ag in dashed-red. Each of the one-hundred gray lines is a sample from the predictive distribution. This distribution interpolates the training data (black dots), GP hallmark, yet our contributions are not limited to noise free settings. (We entertain a noisy f in Section 3.2) The current best run is near x but there is considerable predictive uncertainty about other parts of the input space B, relative fmiin Locations in [-1, 0J U [4, [2T, 2n + are promising. and that BO acquisition criteria serve operationalize this notion. We shall focus [Thomsoon 1933) several others are variations on similar themes. One is Thompson sampling (TS; (TS:
fa min (yi. ..1...., Un .,
and involves drawing from the predictive distribution Y(x) fn X,. Yn and choosing Tn+11 rrgminges Yxx). TS amounts to randomly selecting one those gray lines and optimizing it lieu of the expensive f, so the criterion is inherently stochastic. The other is expected improvement [EI; 1998 Define improvement as ((x) (O, fruin Y (x)} and take its expectation with respect to Y()) If (x)) ~ In is Gaussian, as with GPs, then this has closed form:
Jones et al.
[(x) (x ( - pn (x)) ((E))))()))) + (x) ((n(()) BI(x) = B([(X))
where 2, (z) = (fmmin (z and //0 are the Gaussian cdf/pdf. For non Gan--Gassinn fa. one can always resort Monte Carlo (MC) integration instead. The right panel of Figure S provides EI for fo. Contrary to TS, the El acquisition Øn+1 = argmax El(x) can be deterministic, properly solved. The maximum El acquisition x7 is shown green diamond on the left panel.
Observe that TS and EI involve "inner-optimizations over criterion, task that may be even more challenging than the original problem O Each gray line is highly multi-modal, as is the El surface. Both have about as many local optima as there are training data points n, whereas f has only two local minima in B. Speedy evaluation of μr, (z) and a, (z) relative to f(x) saves us, but only partly. We still desire good F1+1 without exhaustive search or cumbersome subroutines.
Two strategies are common, sometimes separately, sometimes in tandem as a hybrid. The simplest option to distribute candidate points X'x throughout the input space, evaluate the criterion on XN, and thereby replace a continuous search with a discrete one. In low input dimension a dense grid of candidates be used effectively. In higher dimension one can populate X' with a space-fillingg
design like random Latin hypercube sample [LHS; |LHS: et al. 1979 to manage the computational expense of evaluating the criterion exhaustively. A higher-powered approach to locally apply smooth, convex optimization library such as L-BFGS-B [Byrd et al. 2003). Derivatives may be approximated by miite-differencing or autograd (Paszke al. 2017]. or may have a simple closed form depending on the criterion and nature f.. The former requires more computer work and some consideration of numerical stability; the latter more researcher/programnee effort when possible.
(MC-based f. or El challenges both approaches.) Pure maddidate-based search is most common with

---

TS because of its stochastic nature. Gradient- continuous search is popular in simple GP/EI setups, but multi-start scheme is essential avoid inferior local solutions. This is where the hybrids come in: candidates seeding local solvers.
In this paper contend that both strategies, random candidates and multi-start local optimization, can be replaced by (or hybridized with) more thoughtfully chosen Xx. In the 1d setting of Figure S we could place candidates the midway points between each of the existing n inputs and the boundary B, implementing a kind of multi-pronged bisection search [Brrden and Faires 1985, Section 2. I] and resulting in
+ 1 7 candidates. The best of those Xy by either criterion may not give precise solution to the inner optimization, but it would be an effective one because EI and many of the random gray lines indicate a solution close to one of those midway points. Such locations might be much better than ones identified by a limited candidate or numerical local search.
This 1d example is overly simplistic. Going forward, we shall explicitly target and higher dimensions. In Section 2 we scale-up the midway candidate idea, also suggested by |Soott et al. [20IT] what we call "tricands", based on Delaunay triangulation and the convex hull of X,, We explore tricands' features and limitations and suggest remedies with the BO application in mind. In Section 3 we demonstrate that tricands outperform both random LHS candidates and multi-start gradient-based numerical inner optimization in the conventional setting where surrogate GP predictive equations are available closed form. In Section we consider two nonstationary surrogates requiring Markov chain Monte Carlo (MCMC), for which closed-fon forem acquisition criteria are readily available. Candidates are essential this setting, and our tricands are better than random space-filling ones. Our discussion in Section 5 emphasizes tricands' "plug-n-play" nature- they can inserted into any mandddte-bassed scheme and suggests potential for extension.
2 Delaunay triangulation candidates
Many criteria for BO resemble Eq. 2 balancing exploitation (4s, (x) below fmii) with exploration (large o, (x)). Most surrogates inflate predictive uncertainty (a,, (I). away from training data locations Xn- In the case of GPs, this is what produces the "sausage shaped" predictive intervals shown as red-dashed lines in Figure Our main insight is that careful allocation of candidates between existing training data locations, where on (z) is high. allows for BO acquisitions that do not necessitate cumbersome numerical optimization of posterior predictive quantities but still balance exploitation and exploration. A hard statistical optimization can be replaced with an easier geometric one.
1.0
0 design
triangles fringe
cands
0.8
0.8
66
0.6
x2
0.4
0.4
0.2
0.2
0.0
0.0
random selections
0.0
0.4
0.6
0.8
0.2
0.4
1.0
Figure 2:
Interior and fringe candidates (both), and randomly sub sampled candidates (right).
The idea is sketched graphically Figure whose details will emerge over the next three subsections. The existing design X,, is shown open circles, and the selected candidates are closed red dots and the tips of dashed arrows. call these interior (Section 2.1 and fringe (Section 2.2 candidates, respectively. Both are calculated based on Delaunay triangles and their convex hull, outlined by solid black lines the figure, and explained next. We provide this illustration 2d to ease visualization: however our Section benchmarks include higher dimensional input spaces.

---

2.1 Interior candidates
Delaunay triangulation of X, is angle-maximizinn set of d line-segments connecting geographi- cally nearby points, such that no point lies inside the circumcirele of those points. Such triangulations subdivide the interior the conver hull of Xn, which is the the smallest (convex) set that contains all points. In 2d, lines depicting those subdivisions form triangles. In higher dimension they form tetrahedra, however one often abuses the nomenclature and still refers to triangles. For more details, see, e.g., Lee and Schachter [1980]. The solid lines of Figure Lindicate a triangulation for random X, (left) and an X, derived after BO iterations (right), described in Section B
There are fast algorithms for calculating Delaunay triangulations. In 2d, these have On log n) runtimes. Higher dimensional analysis complicated by the number of triangles, which depends on the geometry of X,. a topic we shall return to shortly. For R we use the geometry package on CRAN [Habel al. 2019; ; for Python we Delaunay in scipy. spatial [Virtanen al. 2020) Both are wrappers around the C library Qhull, implementing "quickhull" [Barber et al.) 1996. Bates and Pronzato 2001 first suggested Delaunay triangulation for BO. That early work only explored two input dimensions, possibly because they did not have convenient access to Qhull. They also did entertain the BO-specffic extensions we provide here, particularly in Sections 2.2
Let the triangles be denoted by 7j. for j = 1,..., nr. Each 7j is (d + 1) x matrix when X,, has columns. Create Mrr new candidates Xx, where the jth candidate located at the barycenter of T
3 = T, d Σ TE.
1 8yk
Tili, A., for k i,...:d.
d
The left expression is vectorized over the second, column dimension of Ty. The second is explicit about coordinates 2 (Z.1....., Zjd). Red dots in Figure provide visual. These locations will almost certainly not be the maximal points of o, (z)) in the vicinity of T,, but they will be close because ž is within T but far from its edges. We refer to these Xer as "interior" candidates.
In 2d, Euler's formula gives that
2n 2 - ((X,,) where h (X,.) is the number of elements of X, on its convex hull. In Figure 2 n = 10 and ((X,,,) = 6 so nt 12. When 2 3, the number of faces of the tetrahedra can grow as n [d/2] depending on the nature of X,. Figure provides
kog. 10 elappeed Becondss
1000
Internal
input dimension
Input
Figure 3: Tricands N (left) over input dimension d and design size n, and compute time (right) using single core of an Intel 17-6900K CPU at 3.20GHz.
an empirical view of the number of candidates (left) and triangulation compute time (right), for varying d and where X,, is distributed uniformly at random in B = [0, 1]d. The dashed lines in the left panel are log nr. Notice that when n is small, ny decreases steadily increases with n. In our BO context this is good thing, because so too does the modality of acquisition criterion. With nonparametric surrogates like GPs, the complexity of the response surface can increase with training data size n. The search effort for the next acquisition should be commensurate with this complexity, an innate characteristic of (interior) tricands. The right panel of Figure 2 described later in Section 2.3 shows how this works after many iterations of EI acquisition.
increases. For fixed d, ny

---

2.2 Fringe candidates
Having candidates Xxp only the interior of the convex hull could limit exploration if the comple- ment of the volume of the hull and B is large. It could limit exploitation if the solution is on/near the boundary of B. One remedy input space. This is not bad approach in low dimension (e.g.. in 2d there are only four corner), but could be prohibitive in higher dimension. Later Section 3 we consider a = 8 example, having 256 corners in B = [0, 1J8, which would nearly blow our entire budget of runs.
to force X,, to contain boundary points, such as the corners of the
We instead prefer the "fringe" candidates pointed by the dashed arrows in Figure 2 There one of these for each facet (edge in 2d) of the convex hull, extending perpendicularly from middle of the facet half way to the boundary of B = [0. 1] The Qhull library furnishes both ingredients: d x d facets Fj and normal d-vectors ūj for j = I..... Rg. Using these quantities, let F, = à/ Ei F2 Fy]
denote the coordinates of the middle of facet j. Then, the jth fringe candidate in Xng is
- Fg + +jayū where aj = min {I - -F
Above, "min" is picking out the nearest boundary to Fj. Division by two mirrors the bisection search analogy for interior candidates, but this could be tuning parameter. Non-unit rectangular B is doable, but requires a more convoluted formula. More general B may present challenges.
Fringe and interior candidates may be combined: Xy = i I. stacked row-wise to form an N x d matrix with N ny + n.p. The number of fringe candidates n is generally small compared to ar. This is shown empirically in the left panel of Figure 3 where N is indicated by the solid line, and nr by the dashed one. As n increases, the gap between N and ne, indicating nr, narrows.
2.3 Targeted sub-sampling
Having N grow exponentially in n may not suit all applications. Entertaining N R 10,000 candidates when 100 and d = 6, referring to Figure 3 is cumbersome and potentially overkill. One way to limit N is random sub-sampling: simply calculate the full Xy based on X,, and downsample, uniformly at random, a subset of size NNub < N. Since our triangulation strategy was designed to focus on exploration (finding locally high o (x)) via midway candidates, we find it advantageous to guarantee retaining some of those candidates which are promising for exploitation (low pn (z)), still without any explicit numerical optimization - only geometry.
One of the rows of Xn corresponds to fmin, the best input found so far: xminn = s.t. i argmin,=11.... Vi (or pa. (z:) in the noisy case). Let T.min {T, x E ...., denote the set of triangles containing a-min, and similarly let A = {Ø, E Xx : E Triin) denote the candidates associated with those triangles. Those are, in a geometric sense, adjacent to amin. A fringe candidate may also be considered adjacent in an analogous way, however we place them in the complement XN \ X Now, rather than sub-sample uniformly from the full set Xx. we partition sampling from points adjacent to a-i i.e., from X'min and from points farther afield in Xy X In so doing, we guarantee that our Neub candidates cover potential for exploitation and exploration, respectively. We prefer a 10:90 split, with up to 10% of Neub coming from candidates adjacent to x, fewer if [2, < Naub. and likewise 90% from its complement.
The right panel of Figure 2 providee an illustration after n = 30 runs optimizing the Goldstein-Pric function with El under random initialization (details in Section 3.1 When n 30 we have N ~ 60, with the precise value depending on hxX,,), Here we consider Nuub = 30. Observe how randomly sub-sampled candidates X' Nan.. circling the original red candidates Nx in green, concentrate near the global minimum (0.5, 0.25) because --تاار is in the vicinity. Other sub-sampled candidates are spread out more widely. This figure illustrates how acquisitions, and thus candidates, gravitate toward promising regions for exploitation without neglecting areas of potential exploration. A ridge of local minima may be found in the sanana-shaped region traced out by a concentration of X,, and XN-.
2.4 Implementation and software
Our implementation, provided for Python and R in our git repository is relatively tidy. For example, trids R therein contains just 71 lines of code, eleven of which are to support optional
5

---

visualizations in 2d such as those in Figure _ Delaunay triangulations, we provide option "Q12" work around numerical instabilities that sometimes arise. When calculating convex hulls, option "n" " returns normal vectors used in calculating fringe candidates (Section 2.2
The heavy lifting is done by Qhull. When performing
Defaults yield both fringe and interior candidates and fix a maximum candidate size of max NNuub 100d, but these are user-adjustable. In experiments coming shortly, we deliberately limit Naub even further so that a fairer comparison can be made to other continuous search and candidate-based methods. When N S Nuub, all interior and fringe candidates are returned. The user may supplement these with additional Nsub - N random candidates if desired. If Neub < NN, the execution flow looks for a variable best, providing the index of x-min in X,, making sure about 10% of tricands come from adjacent triangles. When best= tricands are sub-sampled at random.
3 Classical GP benchmarking
In this first of two sections on benchmarking, we focus on BO via traditional GP surrogates. We follow a homework problem in Section 7.4 of Grramayy (2020) which piggy-backs off of GP, EI, and TS demonstrations earlier in the chapter. Software and other particulars are relegated to Appendix
A.I All of our examples are fully reproducible using the code provided in our git repository. We consider three methods for solving EI acquisitions: a continuous search of the criterion via L-BFGS-B with 5-multi-starts, LHS candidates, and tricands. For TS acquisitions, we similarly employ both LHS candidates and tricands. As non-BO (not surrogate-based) comparators. serving primarily as benchmarks, we entertain "raw" " L-BFGS-B B and Nelder-Mead [Nelder and Mead 1965..
Two examples are showcased here with a third example relegated to Appendix A.2 Our synthetic fs, including those in Section 4 are described in more detail on the pages of the Virtual Library for Simulation Experiments [VLSE; Surjanovic and Bingham 2013 In all of our experiments we code inputs to [o, 1j and summarize results for 100 random restarts where each surrogate is initialized with the same (unique to each random restart) starting design of size no 12, except in Section 3.2 where we use no = 60. This design taken uniformly at random following the advice of Zhang et al. [021]] who caution that small space-filling initial designs can spark pathological behavior in BO, We track "eet observed value" (BOV) fmin as a measure of progress which is summarized by median over = 1, ....., Rend and by boxplots for particular along the way.
Rend
3.1
Goldstein-Price
The -Plistein-Price function is a popular low-dimensional (2d) benchmark for BO [e.g.. Picheny al. 2012 A total of Rend 50 acquisitions are entertained, and all candidate-based methods (tricands and LHS) are limited to fifty candidates. This means max = Nsub = 50 for tricands, with fewer candidates when N Noub. The top row of Figure g summarizes results in three views. In the top-left panel, median progress in BOV is shown over increasing budgets of evaluations (n) as if each subsequent acquisition were the last. In the top-middle and right panels, boxplots capture the distribution of BOV at n = 30 and = 50, respectively. In the top-left panel, methods based on multi-start numerical local search use solid lines; those based on tricands are dashed; those based on LHS candidates are dotted. (Nelder-Men Mead is an exception, being red-dashed.) In the boxplots, tricands use slightly heavier ink so that they stand out. Text printed at the top of the top-right panel indicates the average number of times the criterion, EI or TS, was evaluated in solving the inner optimization sub-problem(s), cumulative over all acquisitions. In the case of numerically optimized El (labeled "El"), these happen within iterations of L-BFGS-B search.
Beginning with medians over in the top-left panel, observe that the dashed lines (those based on tricands) are uniformly superior to all other comparators. EI with tricands is slightly better than TS with tricands. Only at the very end of the run does the raw Meader-Mead alternative look competitive. At the 30th evaluation (top-middle) the boxplots indicate that there are some MC repetitions where BOV based on tricands are underperforming. However, more than half of the time these are the best two (i.e., both El and TS with tricands) of all. TS with tricands is the winner in terms of best-case performance, whereas the El version has slightly better worst-case behavior. Neither local optimizer (red) is competitive. Finally, in the top-right panel we see that the story is similar, except that now EI with tricands edges out its TS analog. Although median performance for Nelder-Mead competitive, more than half of the time its BOV at nend 50 is among the two worst in the experiment.

---

Goldstein-Prica noise
evaluation
evallation
Neldor-Mead BFGS
Pafoordd sureccitei mvr 1000
0.0
Ei-t1l
-10
2021
Ei-hs
TS-Mi
-03
TS-ihe
medder beet oheerned vabue
.13
-10
-1.
-20
00.00
-2.0
-25
-3.
-30
inttal
B-hhs
TS-Ii
TS-hi
NM
BRGS
TS-4n
TS-hrs
BrQS
avaluations
Assemble Order, input-dependent
evaluation
evaluation
El-hyb
(estimated) uoserved
modlan (ostimato))
El-to
EL-26
biacibox evaluations
Figure 4: Goldstein- (top) and ATO (bottom) BOV over 100 MC trials: median (left panel) and distribution intermediate (middle) and final (right) acquisitions. Suffixes "-tri" or "-lhs" denote candidate search. Absence of a suffix (e.g., "El" alone) indicates multi-start L-BFGS-B.
Perhaps the most striking result from the figure is that tricands- EI ("EI-trii" outperforms its multi-start numerical analog ("El") despite factor of five fewer evaluations of the EI criterion (about 9000 compared to 1700). A more aggressive, gradient-based search misutilizes computational resources. Finding a precise solution ultimately furnishing a maximal local value an inferior domain of attraction may come at the expense of finding an accurate approximation near global solution. The same is true, but lesser extent, when comparing TS variations (5% reduction).
See Appendix A.2 for nearly identical results with the higher dimensional Hartmann 6 function.
3.2 Assemble to order (ATO)
The assemble-to-order (ATO) simulator (Hong and Nelson 2006) deploys a queuing system vir- tualizing manufacturer receiving orders for a variety of widgets, each with different demand for component parts. Target inventory levels for eight parts comprise the controllable inputs = 8). The output is profit, scalar distilling inventory and order fulfillment costs under random orders for parts, fulfillment and inventory replenishment times. ATO is coded Matlab and exhibits noise. Although it is reasonably fast (seconds per evaluation), relatively large n is required to separate signal from noise (no 80 and Aend = 300). Cubic computational bottlnecks for GP inference are in play here, despite using thrifty heteroskedastic GP [Binois et al. 2018) More details, including software, are provided Appendix A. Since our Matlab licenses prevented us from running parallel instances, we had to limit our experiment somewhat to obtain results in reasonable amount

---

of time. Therefore we only entertained three acquisition alternatives: random sulti-start numerically optimized EI, tricands with the default Nuub 100dd 800, and a hybrid using the best tricands candidate to initialize a local numerically optimized EI ("EI-hyb").
The bottom panels of Figure 4 track BOV for these three alternatives in view similar the panels above. Note that the original ATO problem involves maximization, so we have negated the output and also scaled it so that its units were more similar to our other benchmark problems. Since the output is random, the BOVs are estimates. We used the surrogate fit at Rennd to assign in-sample predictions via Mn (at). Otherwise. the story is similar to our earlier synthetic results. Tricands' progress is initially slower than numerically optimized EI, but eventually gives lower BOV values. The difference at Rend is slight, but statistically significant. A paired Wilcox test with alternative hypothesis that BOV for "El-tri" is below "El" has a p-value of 0.0037. The hybrid is even better, beating pure tricands and pure numerical optimization 85 and 88 out of 100 times, respectively.
4 Sampling-based surrogates
Tricands are most valuable when the inner-optimization problem cannot be solved by library- numerical methods, even locally. This happens when the surrogate predictive surface is discontinuous and/or when inference requires MCMC. Our examples here involve response surfaces exhibiting nonstationarity, meaning that the input-outpe dynamics evolve over the input space. This demands a more elaborate surrogate. We consider two. A treed Gaussian process [TGP: /Graaayy and Lee 2008) uses uxis-aligned partitioning with GPs. Such divide-and and-conquer excels when dynamics change abruptly across individual inputs creating distinct regimes. MCMC can average over the location of probable partition boundaries. We use the R package tgp on CRAN [Gramacy 2007]. following
Section 4 of Gramacy and Taddy 2010 for candidate- based EI acquisition.
Deep Gaussian processes [DGPs; Damianou and Lawrence, 2013 warp inputs to accommodate a more subtle evolution of dynamics in the input space compared the abrupt regime changes of TGP. Although variational inference is popular for DGPs (Salimbeni and Deisenroth 2017], Sauer 2020 argue that in active learning contexts, such as BO, full posterior integration leads to better uncertainty quantification, and thus better acquisitions. Here we use the R package deepgp on CRAN (Suuer, 202T] which supports EI evaluation on candidates.
al.
Unfortunately, neither tgp nor deepgp facilitate TS acquisition for BO. However, the tgp package furnishes maximum a posteriori sample which can be used as the basis of local search [Gramacy and Taddy 2010 Section 4]. Here we show that tricands provides better candidates for both pure candidate El search and this hybrid scheme.
4.1 Abrupt changes
The Gramacy & Lee (G&L) function benefits from a model supporting hard breaks even though its dynamics evolve smoothly. The 2d input domain (coded to [0, 1J²) is mostly flat except in one quadrant where a local maximum is twinned with a (global) minimum. The domain of attraction of that global minimum covers only about 10% of the input space. It is easily missed by random initial designs and candidate sets. TGP is able to isolate the interesting quadrant with just two axis aligned partitions, thus recognizing that sampling effort should be concentrated there.
The top row of Figure shows results the same three views as earlier. The number of candidates is limited to twenty, otherwise the setup is unchanged. First ignore the red lines and boxplots (those labeled "hyb-") and focus on the black (TGP) and green (DGP). Notice that in both cases the tricands- based comparators, dashed and/or bolded, outperform their LHS-based analogues. In terms of medians over n (top-left panel), that dominance is uniform. In terms of boxplots, the disparity is stark with the exception of few outlying DGP- BOV values with tricands. TGP seems to edge out DGP (with tricands) which we attribute to the abrupt change between the interesting quadrant and the rest of the input space. Now focus on the red "hyb-" results. These the ones where TGP candidate-basee
search is finished with a gradient- based EI search on the maximum a posteriori model. Again, tricands wins. To supplement the visuals, we report that the tricands hybrid BOV value was lowest in 97/100 reps (top-right).
A higher dimensional example using the Michaclwicz function is provided in Appendix A.3

---

Gramacy no
evaluation
00
TGP-Ihs
-3
TGP-tr
Nyb-lha
-032
Ryb-Jni
-0.1
DGP-ths
meed bret desern vet
DGP-ti
-833
beet obseroed
.005
-0.3
-040
+0
-04
-142
+
1
20
TGP-IN
TGP
hyc-tns
hyb-tn
DGP-hhs
DGP
TGP-Ins
TGP-tn
hyc-hrs
hyb-in
DQP-46
Levy evaluations
evalustion
evaluation
TGP-hs
T0P-31
hy0-ths
ty0-tri
DGP-ihs
modon best Cho0ree vollo
DGP-
best otserver
0.10
8002
TGP-9
ny-
hyp-ti
DBP ho
DOP-ai
TGP-Bs
TGP-3
1/o
DGP- hi
biscibox evaluations
Figure 5: Gramacy & Lee (top) and 5d Levy (bottom) with MCMC-based comparators.
4.2 Smooth changes
The Levy function is defined in arbitrary dimension. is ruffled with many peaks and valleys. Here we consider d = 5 for variety; the bottom row of Figure 5 Commensurate with our other examples, we consider two hundred candidates and nend 75.
Both candidate DGP variations outperform their TGP counterparts. Slowly/smoothly varying dy- namics favor warping inputs as opposed to hard partitioning. This is evident in both median and full distribution (boxplot) views. In the median view (bottom-left), notice that the dashed lines are uniformly under the solid ones of the same color. Tricands are providing better median performance over all n. Although the preeminence of tricands is apparent visually, report that it provided better BOV in 81/100 MC restarts among the best (DGP) comparators. When tricands are involved, the hybrid candidate/numerica option not discernibly better than the pure candidate alternative.
5 Discussion
We offer a novel take on space- space-filling candidates for acquisition in BO. The idea is to fill the spaces in-between previous acquisitions. This is motivated by an analogy to bisection search, but also by the nature of GP predictive surfaces which often serve surrogates in BO. GP surrogates have organically inflated uncertainty between training data sites, which makes those spots attractive BO acquisition. This notion is extended to the space between the convex hull of existing training data and the boundary of the study region. In array of benchmark exercises we have demonstrated that these tricands lead to superior performance in BO compared to both higher-powerey tradient-based

---

acquisition schemes and simpler space- space-filling candidates. Tricands' main attraction is that they mimic the behavior of nigher-powered searches with the implementation simplicity of candidates. That simplicity means that tricands may be deployed where the higher- alternatives cannot, such as when the surrogate not continuous or requires MCMC.
Additional discussion on high input dimension and other extensions is provided in Appendix
Acknowledgments and Disclosure of Funding
We are grateful for funding from DOE LAB 17-1697 via subaward from Argonne National Laboratory for SciDAC/DOE Office of Science ASCR and High Energy Physics, and for comments from two teams of anonymous referees which greatly improved the paper.
References
Dario Azzimonti, David Ginshourger, Clément Chevalier, Julien Bect, and Yann Richet. Adaptive design of experiments for conservative estimation of excursion sets. Technometrics, pages 1-14, 2020.
Evan Baker, Pierre Barbillon, Arindam Fadikar, Robert Gramacy, Radu Herbei, David Higdon, Jiangeng Huang, Leah R Johnson, Pulong Ma, Anirban Mondal, et al. Analyzing stochastic computer models: A review with opportunities. arXiv preprint arXiv:2002.0133333 2020.
C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa. The quickhull algorithm for convex hulls, ACM Trons. Math. Softwe 22(4):469 December 1996. ISSN 0098-3500. doi: 10. 10.1145/235815223882 URL.
httppa://do..org//0011113 0151000000001111155555515555 235821
R. Bates and L. Pronzato. Emulator-based global optimisation using lattices and delaunay tesselation. In Sensitivity Analysis of Model Output, pages 189-192, Madrid (Espagne). 2001.
Julien Bect, David Ginsbourger, Ling Li, Victor Picheny, and Emmanuel Vazquez. Sequential design of computer experiments for the estimation of probability of failure. Statistics and Computing, 22(3):773-793, 2012.
Henrik Bengtsson, R.matla:: Read and Write MAT Files and Call MATLAB from Within R, 2018. URL
httpss//OtAi...prrjjecc...m wwtpprooecct.orrppccageeeeeeeeee
https://CRAN www.oRppppe.oooom
matlab R package version 3.6.2.
James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for hyper-pane parameter optimization. Advances neural information processing systems, 24, 2011.
Binois, Gramacy, and M Ludkovski. Practical heteroscedastic Gaussian process modeling for large simulation experiments. Journal of Computational and Gruphical Statisticx, 27(4):808 -2188 2018. doi: 10,1080/1061 101100000018600201110 8600.20118 1458625. URL BtRlphtts:///oo.. org/ 10. 1080 1ww186001160... 2018. 1458625
M Binois, J Huang, RB Gramacy, and M Ludkovski. Replication or exploration' Sequential design for stochastic simulation experiments. Technometrics, 27(4):808 821, 2019. doi: 10.1080000401700.220 14694333...3......11111113333331 URL
2018. 1469433
Mickael Binois Nathan Wycoff. survey on high nigh-dimensional Gaussian process modeling with application
to bayesian optimization, 2021.
Mickaël Binois and Robert B. Gramacy. hetgp: Heteroskedastic Gaussian process modeling and sequential design in r. Journal of Statistical Software, 98(13):1 44, 2021. doi: 10. 1118637/js..0988.113. URL https
www..jstatss.. org/ index. http://sss/rrtttll//////////////
Leo Breiman. Random forests. Machine learning, 45(1):5-32, 2001.
AD Bull. Convergence rates of efficient global optimization algorithms. Journal of Machine Learning Research, 12(0c))2879-2904, 2011.
R.I.. Burden and J.D. Faires. Numerical Analysis. PWS Publishers, 1985.
Richard Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16, 02 2003, doi: 10.1 10.1137/091606,
Rob Carelll. 1hs: Latin Hypercube Samples, 2018 URL https:://RAtt R- package version 0.16.
R
10

---

Clément Chevalier and David Ginsbourger. Fast computation of the multi-points expected improvement with applications in batch selection. In International Conference Learning and Intelligent Optimization, pages 59-69. Springer, 2013.
Clément Chevalier, Julien Bect, David Ginsbourger, Emmanuel Vazquez, Victor Picheny, and Yann Richet. Fast parallel kriging-hased stepwise uncertainty reduction with application to the identification of an excursion set. Technometrics, 56(4))455-465, 2014.
Austin Cole, Robert B Gramacy, James E Warner, Geoffrey Bomarito, Patrick E Leser, and William Leser. Entropy-basse adaptive design for contour finding and estimating reliability, arXiv preprint 2rX.v:2105.113:. 2021.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In Artificial intelligence and statistics. pages 207-215. PMLR, 2013.
Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy, Multi-objective Bayesian optimiza- tion over high-dimensional search spaces, 2021.
David Friksson, Michael Pearce, Jacob Gardner, Ryan Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, pages 5497-5508, 2019.
Matthias Feurer, Benjamin Letham, Frank Hutter, and Eytan Bakshy. Practical transfer learning for bayesian optimization. arXiv preprint arXiv:1802.0223 2018.
Peter Frazier, Warren B Powell, and Savas Dayanik. A knowledge-gradice policy for sequential information collection. SIAM Journal Control and Optimization, 47(5):2410-2439, 2008.
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2022. in preparation.
David Ginshourger, Redolphe Le Riche, and Laurent Carraro. multi-points criterion for deterministic parallel
global optimization based on kriging. NCP07, 2007.
J Gonzalez, M Oshorne, and Lawrence. GLASSES: Relieving myopia Bayesian optimisation. In Proceedings of the 19th International Conference on Artificial Intelligncce and Statistics, pages 790-799, 2016.
Javier González, Zhenwen Dai, Philipp Hennig, and Neil Lawrence. Batch bayesian optimization via local penalization. In Artificial intelligence and stotistics, pages 648 PMLR, 2016.
RB Gramacy. tgp: R package for Bayesian nonstationary, semiparametric nonlinear regression and design by treed Gaussian process models. Journal of Statistical Software, 19(9).6, 2007.
RB Gramacy, 1aGP: large-scale spatial modeling via local approximate Gaussian processes in R. Journal of Statistical Software, 72(1):1-66 2016.
Gramacy and Herbert KH Lee. Optimization under unknown constraints. In Bayesian Statistics, volume Oxford University Press, 2011.
RB Gramacy and HKH Lec. Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, (03(483):1119-11) 2008
RB Gramacy and F Sun. laGP: Local Approximate Gaussian Process Regression, 2018. URL BtlLLttp:://ee http://bobby gramacy. com/ http:////aaaaaaaaaaaaaaaaa R package version 1.5-3.
RB Gramacy and MA Taddy. Categorical inputs. sensitivity analysis, optimization and importance tempering with tgp version 2, an R package for treed Gaussian process models. Journal of Statistical Software, 33(6); 1-48, 2010.
Robert B. Gramacy. Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Chapman Hall/CRC, Boca Raton, Florida, 2020. hstt:://bobbyggrama www.eacccc..oooooeeeeeeeeeeeeeee surrogates/
Habel, Raoul Grasman, Robert B. Gramacy, Pavlo Mozharovskyi. and David C. Sterratl, geometry: Mesh Generation and Surface Tessellation, 2019. URL http://ccA...-ojjccetooooo ottpooieet.rt///pcaaggggeeennttt R package version 0.4.5.
L.J. Hong and B. L. Nelson. Discrete optimization via simulation using compass. Operations Research, 54(1): 115 125,129,20000 2006.

---

Donald R Jones. Mathhiss Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global optimization, 13(4):455-492 1998.
Remi Lam, Karen Willcoo, and David Wolpert. Bayesian optimization with finite budget: approximate dynamic programming approach. Advances in Neural Information Processigg Systems, 29:883-891, 2016.
ER Leatherman, TJ Santner, and AM Dean. Computer experiment designs for accurate prediction. Statitics and Computing, pages 1-13, 2017.
Der-Tsai Lee and Bruee J Schachter. Two algorithms for constructing delaunay triangulation. International Journal of Computer Information Sciences, 93))219-242,, 1980.
Benjamin Letham and Eytan Bakshy. Bayesian optimization policy search via caline-offline experimentation. Mach. Learn. Res., 20:145 2019.
Alexandre Marques, Remi Lam, and Karen Willcox. Contour location via entropy reduction leveraging multiple information sources, In Advances in Neural Information Processing Systems, pages 5217-5227, 2018.
D Mckay, Richard Beckman, and William Conover. A comparison of three methods for selecting vales of input variables in the analysis of output from a computer code. Technometrics, 21:239 245, 1979.
Jonas Močkus. On bayesian methods seeking the extremum. Optimization techniques IFIP technical conference, pages 400-404. Springer, 1975.
JA Nelder and R Mead. A simplex method for function minimization. The Computer Journal, 7(4):308-31 313, 1965.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang. Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation pytorch. NeurIPS Autodiff Workshop, 2017.
Picheny, T Wagner, and D Ginsbourger. A benchmark of kriging- based infill criteria for noisy optimization. Structurul and Multidissiplinuary Optimization, pages 1-20, 2012.
Tony Pourmohamad and Herbert KH Lee. Bayesian Optimization with Application Computer Experiments. Springer, New York, NY, 2021.
Pritam Ranjan, Derek Bingham. and George Michailidis. Sequential experiment design for contour estimation from complex computer codes. Technometrics, 50(4):527-541, 2008.
Hugh Salimbeni and Mare Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. arXiv preprint /tt.:::/7550033333 2017.
Annie Sauer. deepgp: Sequential Design for Deep Gaussian Processes using MCMC, 2021. URL https package version 0.3.0.
Annie Sauer, Robert Gramacy, and David Higdon. Active learning for deep Gaussian process surrogates.
arXiv preprint arXiv:2012.08015, 2020.
Warren Scott, Peter Frazier, and Warren Powell. The correlated knowledge gradient for simulation optimization of continuous parameters using Gaussian process regression. SIAM Journal Optimization, 21(3):996 1026. 2011.
Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995 2009.
Peter Su and Rohert L Scot Drysdale. A comparison of sequential delaunay triangulation algorithms. Computa- tional Geometry, 765-6):311-385, 1997.
S Surjanovic and D Bingham. Virtual library of simulation experiments: test functions and datasets. http: //www. sfu. ca/ jsurjoooo 2013.
MA Taddy, HKH Lee, GA Gray, and JD Griffin. Bayesian guided pattern search for robust local optimization. Technometrics, 51 (14)33899401, 2009.
William R Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3/4):285225 294, 1933.
12

---

Ryan Turner. David Eriksson. Michael McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu. and Isabelle Guyon. Bayesian optimization is superior random search for machine learning hyperparameter tuning: Analysis of 10201, the black-box optimization challenge 2020. arXiv preprint arXiv:2104. 102001, 2021.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau. Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright. Stéfan J. van Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, Ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, lan Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing Python. Nature Methods, 17:261-272, 2020. doi: (038/64192--019-0666
Linnan Wang, Rodrigo Fonseca, and Yuandong Tian. Learning search space partition for black-hox optimization using monte carlo tree search. Advances in Neural Information Processing Systems, 13:1951 1-19522, 2020.
Boya Zhang, Robert Gramacy, Leah Johnson, Kenneth A Rosc, and Erie Smith. Batch-sequential design and heteroskedastic surrogate modeling for delta smelt conservation. arXiv preprint arXiv:2010.065/5, 2020a.
Boya Zhang, D Austin Colc, and Robert Gramacy. Distance-dissributed design for Gaussian process surrogates. Technometrics, 63(1):40-52, 2021.
Yichi Zhang. Daniel Apley, and Wei Chen. Bayesian optimization materials design with mixed quantitative and qualitative variables. Scientific reports, 10(1):1-13, 2020b.
13


Reviews:
[{"review_id": "vr-sE7SztlT", "reviewer": "Reviewer_xvts", "summary": "The paper aims to improve the optimization of an acquisition function in Bayesian optimization (BayesOpt).\nThe motivation stems from the observation that a typical acquisition function trades off between exploration and exploitation and therefore is maximized in a region between or far away from observations, but never close to them.\nThe authors propose a scheme to generate candidates for the search to optimize the acquisition function that relies on Delaunay triangulation, generating two sets of candidates: interior points that lie between observations and fringe points that lie outside the convex hull of the observed set.\nTo improve efficiency, the algorithm also employs a non-uniform sampling strategy from these two sets to assemble a final candidate set to evaluate the acquisition function with.\nExperiments across benchmark objective functions show that the proposed algorithm leads to better optimization performance for two common acquisition functions: expected improvement and Thompson sampling.\nFinally, the paper investigates situations in which the predictive model produces predictions that make the acquisition function unamenable to gradient-based optimization, such as treed or deep Gaussian processes, showing that their algorithm consistently improves performance across the different models.", "questions": "The motivation of the algorithm, looking between the observations as well as regions far away, feels to me similar to a space-filling objective where we'd like to, in addition to the observations already made, cover the rest of the space.\nHave the authors considered a baseline of this kind where it generates points that together with the points in the training set, infills the space and uses generated points as candidates for optimizing the acquisition function?\nI expect this strategy would bias regions far away from observed points and might over-explore, but it may still pose as a spiritual competitor to highlight the benefits of the proposed algorithm and its subsampling scheme.\n\nOn the topic of subsampling, a potential enhancement could be to avoid regions with high predictive mean (in the minimization problem in our case) when sampling the interior points, since those points most likely do not maximize the acquisition function anyway.", "limitations": "Yes, the paper mentions its limitations in the appendix.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This is a good submission.\nThe exposition flows well, and the question of how to optimize the acquisition function in BayesOpt is important yet often ignored.\nThe proposed algorithm is natural and well-motivated.\nIt is clear from their experiment results that their policy consistently performs well across many tasks compared to other baselines.\nI also enjoyed reading the addition discussions in the appendix, which make connections to other problems and areas.\n\nI don't have any major complaints about this submission.\nOne thing that the paper could benefit from is experiments with higher-dimensional ($d > 20$) objective functions.\nI would be interested to see the performance and running time of the algorithm, as well as whether there's any change in behavior, when there are a large number of dimensions and triangulation is presumably more difficult.\nI talk more about a potential baseline that the authors may consider comparing the algorithm against in the Questions section.\n\nOverall, the submission addresses an important problem and the proposed solution seems to work well on the inspected settings, and I lean towards acceptance.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "SF41gU4NCB", "reviewer": "Reviewer_eznT", "summary": "The paper introduces a method for optimizing the acquisition function during Bayesian optimization. The method is based on geometry, rather than numerical optimization or random search, the two most common techniques. It generalizes the well-known \"midway candidate\" idea that has been used as an acquisition optimization initializer by at least 2 software packages, and turns it into a complete strategy that works in full-dimensional Bayesian optimization problems. The results with the new method are very promising and outperform gradient optimization.", "questions": "The questions are those above: How do total wall times compare for the methods, and how does tricands perform without targeted subsampling?", "limitations": "Yes", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is very well written and addresses an important problem. Writers of Bayesian optimization software have put a level of effort into this problem that is not fully reflected in the literature, and so it is easy to underestimate the impact of this paper. But it tackles a very important problem, and one that underlies all of the results and progress in the field. It thus has broad impact in the entire space of Bayesian optimization.\n\nThe method presented is exceptionally innovative. The general idea of the \"midway candidate\" is well-known and has been used as an intializer for gradient optimization, but previous applications of this idea have been limited and not extended naturally to higher dimensions, which is why it has only been successful as an _initializer_ for a subsequent, expensive gradient optimization. The paper uses concepts from geometry to turn this idea into a method that applies to the full range of Bayesian optimization problems, which will be especially useful for extended Bayesian optimization to surrogate models and acquisition functions in which gradients may be difficult to compute. The level of novelty in this paper is very high, for Bayesian optimization particularly.\n\nThe results are very impressive and show that we can get _better_ performance than gradient optimization without having to compute gradients. This is a very strong result that I expect will change the way people think about acquisition optimization.\n\nThere are three areas where I think the paper could improve:\n\n* Evaluation is given in terms of number of candidates evaluated, which is a great thing to directly compare. But it would also be helpful to see a comparison in the overhead (for gradient optimization this means computing gradients / matrix operations, for tricands it means the triangulation). Seeing the total wall-time for each method would allow that.\n\n* How important is the targeted subsampling? Seeing some results without targeted subsampling to see would be helpful, as well as some study of the 10% number.\n\nBoth of these are minor revisions.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "CX7IIwiv-tU", "reviewer": "Reviewer_czZM", "summary": "The authors consider the use of triangulations to generate space-filling candidates in Bayesian optimization. They consider a few low-dimensional synthetic problems and find that their triangulation-based approach as well as a hybrid approach marginally outperforms EI with numerical optimization.", "questions": "- How do you initialize optimization with multi-start L-BFGS-B and how many restarts do you do? People would generally use a large number of quasi-random candidates here, e.g., 2048, as computing the acquisition function values for simple acquisition functions such as EI is relatively cheap. Using an unreasonably small number of initial points may explain why EI flattens out in Figure 4 (bottom) as it could be that all starting points have acquisition function value zero.\n- Packages like Spearmint and BoTorch provide the option to generate candidates as local perturbations around the best point found so far which often works well in practice. This is often important later on in the optimization run when exploitation becomes more important and finding candidates with non-zero acquisition function value becomes more challenging. Can you comment on why you didn't consider such an approach and how the performance compares to using triangulations?\n- That numerical optimization of the acquisition function (which depends on the posterior) is cumbersome is a bit vague and unmotivated. The acquisition function surface is in many cases be easy to optimize (EI, UCB, etc.), especially in cases where the gradient can be computed analytically. Can you comment on what you find to be cumbersome here?\n- The intuition that we want to select candidates where sigma(x) is high sounds reasonable in the low-dimensional setting as long as it is supplemented with candidates that also favor exploitation as those will be important later in the optimization run. As an alternative to using triangulations, do you think it would be possible to screen candidates based on the minimum distance to the training points as that would provide a cheaper but potentially still effective way to select good candidates?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- This is generally a well-written paper that very clearly explains the background material and motivates the approach.\n- The successes of multi-start L-BFGS-B is indeed sensitive to the choice of good starting points. This is particularly important when using acquisition function like EI where a large majority of the domain may have zero acquisition function value, particularly later on in the optimization run. The authors observe that the use of triangulations lead to a small performance improvement.\n\nWeaknesses:\n- The approach is limited to low-dimensional problems (8 is the largest considered in the paper) as the use of triangulations is well-known to scale poorly with the input space dimensionality.\n- I don't think this paper brings much value to the BO community. In particular, the authors don't compare to alternative approaches for generating candidates, which makes it hard to judge the value of the approach.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "0xjysuA3gN0", "reviewer": "Reviewer_NN2f", "summary": "This paper proposes a new space-filling acquisition function for Bayesian optimization based on Delaunay triangulation of the initial data.", "questions": "-", "limitations": "-", "rating": 6, "confidence": 1, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The idea of using Delaunay triangulation to get new observation locations is intuitively appealing, and the paper clearly explains their method. However the paper, as it stands, is barely adequate. No theoretical analysis justifies their choices, the experimental analysis is unconvincing, and the literature survey is meager.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept