Unexpected Improvements to Expected Improvement for Bayesian Optimization
Authors: Sebastian Ament
Sam Daulton
David Eriksson
Maximilian Balandat
Eytan Bakshy
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Bayesian Optimization
Gaussian Process
Multi-Objective Optimization
TL;DR: We analyze and fix pathologies of Expected Improvement and its variants for single and multiple-objective problems, leading to significant improvements in empirical optimization performance.

Abstract: Expected Improvement (EI) is arguably the most popular acquisition function in Bayesian optimization and has found countless successful applications, but its performance is often exceeded by that of more recent methods. Notably, EI and its variants, including for the parallel and multi-objective settings, are challenging to optimize because their acquisition values vanish numerically in many regions. This difficulty generally increases as the number of observations, dimensionality of the search space, or the number of constraints grow, resulting in performance that is inconsistent across the literature and most often sub-optimal. Herein, we propose LogEI, a new family of acquisition functions whose members either have identical or approximately equal optima as their canonical counterparts, but are substantially easier to optimize numerically. We demonstrate that numerical pathologies manifest themselves in “classic” analytic EI, Expected Hypervolume Improvement (EHVI), as well as their constrained, noisy, and parallel variants, and propose corresponding reformulations that remedy these pathologies. Our empirical results show that members of the LogEI family of acquisition functions substantially improve on the optimization performance of their canonical counterparts and surprisingly, are on par with or exceed the performance of recent state-of-the-art acquisition functions, highlighting the understated role of numerical optimization in the literature.

Unexpected Improvements to Expected Improvement for Bayesian Optimization
Sebastian Ament
Samuel Daulton Meta sdaulton@meta. com
Meta
amentêmeta.com
David Eriksson
Maximilian Balandat Meta bal @mdat com
Eytan Bakshy
Meta derikssonSmetta.com com
Meta ebakshy@meta.com com
Abstract
Expected Improvement (EI) is arguably the most popular acquisition function Bayesian optimization and has found countless successful applications, but its performance is often exceeded by that more recent methods. Notably, El and its variants, including for the parallel and multi-objective settings, challenging optimize because their acquisition values vanish numerically in many regions. This difficulty generally increases as the number of observations, dimensionality of the search space, or the number of constraints grow, resulting in performance that is inconsistent across the literature and most often sub-optiin sub-optimal. Herein, we propose LogEI, new family of acquisition functions whose members either have identical or approximately equal optima as their canonical counterparts, but are substantially easier to optimize numerically. We demonstrate that numerical pathologies manifest themselves in "classic" analytic EI, Expected Hypervolume Improvement (EHVI), as well as their constrained, noisy, and parallel variants, and propose corresponding reformulations that remedy these pathologies. Our empirical results show that members of the LogEl family of acquisition functions substantially improve the optimization performance of their canonical counterparts and surprisingly, are on par with or exceed the performance of recent state-of-the-ard acquisition functions, highlighting the understated role of numerical optimization in the literature.
Introduction
Bayesian Optimization (BO) a widely used and effective approach for sample-efficient optimization of prrusivv-tt-evaluall black-box functions [25] 28), with applications ranging widely between
aerospace engineering (48) biology and medicine 49 materials science [3], civil engineering 44,
and machine learning hyperparameter optimization
666 72 BO leverages probabilistic surrogate
model in conjunction with an acquisition functtin determine where to query the underlying objective function. Improvement-bent based acquisition functions, such as Expected Improvement (EI) and Probability Improvement (PI), are among the earliest and most widely used acquisition functions for efficient global optimization of non-convex functions 42 (88) El has been extended to the constrained [27 [29]. noisy [52), and multi-objective 20 setting, as well as their respective batch variants 16
13 77). and is standard bascline in the BO literature (25) [66.. While much of the literature has focused on developing new sophisticated acquisition functions, subtle yet critical implementation details of foundational BO methods are often overlooked. Importantly, the performance El and its variants is inconsistent even for mathematically identical formulations and, show in this work, most often sub-optimal.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Although the problem of optimizing EI effectively has been discussed in various works, e.g. 25 B1 77,, prior focus has been on optimization algorithms and initialization strategies, rather than the fundamental issue of computing El.
In this work, identify pathologies in the computation of improvement- -aased acquisition functions that give rise to numerically vanishing values and gradients, which - to our knowledge are present in all existing implementations of EI, and propose reformulations that lead to increases in the associated optimization performance which often match or exceed that of recent methods.
Contributions
1. We introduce LogEI, new family of acquisition functions whose members either have identical or approximately equal optima as their canonical counterparts, but are substantially easier to optimize numerically. Notably, the analytic variant of LogEI, which mathematically results in the same BO policy as EI, empirically shows significantly improved optimization performance.
2. We extend the ideas behind analytical LogEl to other members of the EI family, including constrained EI (CEI), Expected Hypervolume Improvement (EHVT), well their respective batch variants for parallel BO, qEI and qEHVI, using smooth approximations the acquisition utilities obtain non-vanishing gradients. All of our methods are available as part of BoToch (6).
3. We demonstrate that newly proposed acquisition functions substantially outperform their respective analogues on a broad range of benchmarks without incurring meaningful additional computational cost, and often match exceed the performance of recent methods.
Vanishing gradients of EI
Fraction of samples with gradient
1.0
0.8
-100
Mijaaive
-15
Observtion
Toombeen
0.6
2200
Pottrier
Acqulstion values for El add Leg!!
0.4
-10
0.2
dm
d=
-10
Laget
lepito
DI Octrobde prccidem
0
20
40
60
80
100
120
Number of training data points
Left: Fraction of points sampled from the domain for which the magnitude of the gradient EI vanishes to 10 10010 as a function of the number of randomly generated data points for different dimensions d on the Ackley function. As n increases. El and its gradients become numerically zero across most of the domain, see App. D.2 for details. Right: Values of and LogEl on quadratic objective. EI takes on extremely small values on points for which the likelihood of improving over the incumbent small and numerically exactly zero double precision for large part of the domain (z [5, 13.5)). The left plot shows that this tends to worsen as the dimensionality of the problem and the number of data points grow, rendering gradient- optimization of EI futile.
Figure 1:
Motivation
Maximizing acquisition functions for BO is challenging problem, which is generally non-convex and often contains numerous local maxima, see the lower right panel of Figure I While zeroth-order methods sometimes used, gradient-basee based methods tend to be far more effective at optimizing acquisition functions on continuous domains, especially in higher dimensions.
In addition to the challenges stemming from non-convexity that are shared across acquisition functions, the values and gradients improvement-baste basse acquisition functions are frequently minuscule in under Gaussian posterior
large swaths of the domain. Although El is never mathematically distribution/ often vanishes, even becoming exactly zero in floating point precision. The same
eexcept at any point that is perfectly correlated with previous noiseless observation

---

applies to its gradient, making EI (and PI, see Appendix A exceptionally difficult to optimize via gradient-based methods. The right panels of Figure T Iilustrate this behavior on a simple one- dimensional quadratic function.
To increase the chance of finding the global optimum of non-convex functions, gradient-hasee based optimization is typically performed from multiple starting points, which can help avoid getting stuck in local optima (70). For improvement-beed -aased acquisition functions however, optimization becomes increasingly challenging as more data is collected and the likelihood of improving over the incumbent diminishes, see our theoretical results in Section and the empirical illustration in Figure I and Appendix D.2 As a result, gradient- based optimization with multiple random starting points will eventually degenerate into random search when the gradients at the starting points are numerically zero. This problem is particularly acute in high dimensions and for objectives with a large range.
Various initialization heuristics have been proposed to address this behavior by modifying the random- restart strategy. Rather than starting from random candidates, an alternative naïve approach would be use initial conditions close the best previously observed inputs. However, doing that alone inherently limits the acquisition optimization to a type of local search, which cannot have global guarantees. To attain such guarantees, it is necessary use an asymptotically space-filling heuristic; even if not random, this will entail evaluating the acquisition function in regions where no prior observation lies. Ideally, these regions should permit gradient- -ased optimization of the objective for efficient acquisition function optimization, which necessitates the gradients to be non-zero, In this work, we show that this can be achieved for large number improvement-basee acquisition functions, and demonstrate empirically how this leads to substantially improved BO performance.
2 Background
We consider the problem of maximizing expensive-to-evalus black-box function firee X + RM over some feasible set X C 1Rd Suppose we have collected data D, = and y; furee (x:) + (x,) and v; is noise corrupting the true function value (x,)),,). The response furue may be multi-output as is the case for multiple objectives or black-box constraints, in which case yi, V, E R We use Bayesian optimization (BO), which relies surrogate model f that for any batch X == [xi..... X,) of candidate points provides a probability distribution over the outputs f(X)
where x; E X
(f(x)),.... (x.,,.,,,(x))) (x,)). The acquisition function then utilizes this posterior prediction to assign an acquisition value to x that quantifies the value of evaluating the points in x, trading off exploration and exploitation.
2.1 Gaussian Processes
Gaussian Processes (GP) 65 are the most widely used surrogates in BO, due to their high data efficiency and good uncertainty quantification. For our purposes, it suffices to consider a GP as a mapping that provides multivariate Normal distribution over the outputs f(x) for any
f(x) ~ (N((x))) E(x)),
P X
IRaM
Σ X2 + STM
(I)
In the single-outcome (M = 1) setting, f(x) ~ (N((xx),, E(x)) with μ X + R and Σ X - S². In the sequential (q = 1) case, this further reduces to univariate Normal distribution: f(x) ~ o(xx( with μ X + R2 and o X R..
2.2 Improvement-basee Acquisition Functions
Expected Improvement For the fully-sequential (q = 1), single-outcome (M 1) setting, "classic" defined as
El 59
El.. (x) [I[(x) - yn v'I+] ox)) (u(x)= (1(x)-))) 3y
(x))
where [:]+ denotes the max(0, operation, y" = max; yi is the best function value observed so far, also referred to the incumbent, h22) dbz) 2T)2), and Φ t are the standard Normal density and distribution functions, respectively. This formulation is arguably the most widely used acquisition function in BO, and the default in many popular software packages.
3

---

Constrained Expected Improvement Constrained BO involves one or more black-box constraints
and is typically formulated as finding m fime,1 (x) such that f S for i f2..... M.
reasibility-weighting the improvement 27
29
is a natural approach for this class of problems:
M
CEIy- (x) E((x)
[f((x) v"I+ II
1=2
where is the indicator function. If the constraints {ff
modeled conditionally inde- pendent of the objective /i this can be simplified as the product of El and the probability of feasibility.
Parallel Expected Improvement In many settings, one may evaluate fime on q > 1 candidates in
parallel to increase throughput. The associated parallel or batch analogue of EI 30 75 is given by
q
(X)
((x)
max {15(x)) (((xx))>")++}
j=l....g
Unlike EI, qEI does not admit a closed-form expression and is thus typically computed via Monte Carlo sampling, which also extends to non-Gaussian posterior distributions [6 75):
q El " (X)
max {15((x)) - y"}+},
(5)
where (x() [xx are random samples drawn from the joint model posterior at x.
Expected Hypervolume Improvement In multi-objective optimization (MOO), there generally is no single best solution; instead the goal is to explore the Pareto Frontier between multiple competing objectives, the set of mutually-optimal objective vectors. A common measure of the quality of finitely approximated Pareto Frontier P between M objectives with respect to a specified reference point E IR is its hypervolume HV(P,r) rectangle bounded by vertices and yı, and À is the Lebesgue measure. An apt acquisition function for multi-bbjective optimization problems is therefore the expected hypervolume improvement
A(Uy.ep (r. yvl). where [r. y:] denotes the hyper-
EHVI(x) = Ex(x) (HV(PU((X]
HV(P, +)
due to observing a batch f(X) := [((x)). ... f(x,)) of q new observations. EHVT can be expressed in closed form if 1 and the objectives are modeled with independent GPs [80), but Monte Carlo approximations are required for the general case (qEHVI) [13.
2.3 Optimizing Acquisition Functions
Optimizing an acquisition function (AF) is a challenging task that amounts to solving a non-convex
optimization problem, which multiple approaches and heuristics have been applied. These include gradient-free methods such as divided rectangles [41. evolutionary methods such as CMA-ES (32.
first-order methods such as stochastic gradient ascent, see c.g., Daulton al. ([5). Wang et al. [75]
and (quasi-)secone order methods (25) such as L-BFGS-B [10]. Multi-start optimization is commonly employed with gradient- baaed methods mitigate the risk of getting stuck in local minima. Initial points for optimization are selected via various heuristics with different levels of complexity, ranging from simple uniform random selection to BoTorch's initialization heuristic, which selects initial points by performing Boltzmann sampling on set of random points according to their acquisition function value 66) See Appendix B for more complete account of initialization strategies and optimization procedures used by popular implementations. We focus gradient-based optimization as often leveraging gradients results in faster and more performant optimization [13
Optimizing AFs for parallel BO that quantify the value of batch of g > points is more challenging than optimizing their sequential counterparts due to the higher dimensionality of the optimization problem - gd instead of d - and the more challenging optimization surface. A common approach to use sequential greedy strategy that greedily solves a sequence of single
simplify the problem point selection problems. For i = 1,..... candidate x; is selected by optimizing the AF for 1, conditional on the previously selected designs (x.....,X; x;-1} } and their unknown observations, e.g. by fantasizing the values at those designs [77]. For submodular AFs, including EI, PI, and EHVI, sequential greedy strategy will attain regret within a factor of 1/e compared to the joint optimum, and previous works have found that sequential greedy optimization yields improved BO performance compared to joint optimization [3] (7)) Herein, we find that our reformulations enable joint batch optimization to be competitive with the sequential greedy strategy, especially for larger batches.

---

2.4 Related Work
While there is a substantial body of work introducing a large variety different AFs, much less focus has been on the question of how to effectively implement and optimize these AFs. Zhan and Xing [8]] provide a comprehensive review of a large number of different variants of the El family, but do not discuss any numerical optimization challenges. Zhao et al. [82 propose combining a variety of different initialization strategies to select initial conditions for optimization of acquisition functions and show empirically that this improves optimization performance. However, they do not address any potential issues or degeneracies with the acquisition functions themselves. Recent works have considered effective gradient-based approaches for acquisition optimization. Wilson et al. [77] demonstrates how stochastic first-order methods can be leveraged for optimizing Monte Carlo acquisition functions. Balandat et al. [6] build on this work and put forth sample average approximations for MC acquisition functions that admit gradient-based optimization using deterministic higher-order optimizers such as L-BFGS-B.
Another line of work proposes to switch from BO to local optimization based on some stopping criterion to achieve faster local convergence, using either zeroth order [60] gradient-base [57] optimization. While McLeod et al. [57] are also concerned with numerical issues, we emphasize that those issues arise due ill-conditioned covariance matrices and orthogonal to the numerical pathologies of improvement-bart acquisition functions.
3 Theoretical Analysis of Expected Improvement's Vanishing Gradients
In this section, shed light the conditions on the objective function and surrogate model that give rise to the numerically vanishing gradients in EI, as seen in Figure [[] In particular, we show that BO algorithm closes the optimality gap f - y², where f* is the global maximum of the function five. and the associated GP surrogate's uncertainty decreases, El is exceedingly likely to exhibit numerically vanishing gradients.
Let Px be a distribution over the inputs x, and f ~ Py be an objective drawn from a Gaussian process. Then with high probability over the particular instantiation of the objective, the probability that an input x ~ Px gives rise to an argument (μ(x) y²) /ar()xx to h in Eq. 2 that is smaller than a threshold B exceeds (x) < I tx)). where e, depends on the optimality gap '" y² and the maximum posterior uncertainty maxx on (x). This pertains to El's numerically vanishing values and gradients, since the numerical support S,, (h) = fx: : [/(x]] nu} of naïve implementation of in 2 is limited by a lower bound В(π) that depends on the floating point precision 77. Formally, (h) [B)), oc) even though (h) IR mathematically. As consequence, the following result can be seen as a bound on the probability of encountering numerically vanishing values and gradients in El using samples from the distribution Px to initialize the optimization of the acquisition function.
Theorem 1. Suppose f is drawn from Gaussian process prior Pj, y² S f., pns On are the mean and standard deviation of the posterior Pf(f(Dm) and B E R. Then with probability 1 --8, 8,
- y² Mhn (x)
Px
<B 2 <
(x)
where = (f" -f* + (V- logg)))))) B) maxx on (x).
For any given - and especially carly- iteration, en does not have to be small, as both the optimality gap and the maximal posterior standard deviation can be large initially. Note that under certain technical conditions on the kernel function and the asymptotic distribution of the training data D,,, the maximum posterior variance vanishes guaranteeably as n increases, see [50] Corollary 3.2]. On its own, Theorem [ gives insight into the asympttiic behavior by exposing dependence to the distribution of objective values f. In particular, if the set of inputs that give rise to high objective values (< f*) is concentrated, PxJ(x) < f e) will decay very slowly as e increases, thereby maintaining a lower bound on the probability of close to As an example, this the case for the Ackley function, especially the dimensionality increases, which explains the behavior in Figure
5

---

4
Unexpected Improvements
In this section, propose re-formulations of analytic and MC-based improvement-basee based acquisition functions that render them significantly casier to optimize. We will use differing fonts, e.g. log and log, to differentiate between the mathematical functions and their numerical implementations.
4.1 Analytic LogEl
Mathematically, El's values and gradients are nonzero the entire real line, except in the noiseless case for points that are perfectly correlated with previous observations. However, naive implementa- tions of h are numerically zero when 2 (u(x) - y²) (x) small, which happens when the model has high confidence that little improvement can be achieved at x. We propose an implementation of
log o that can be accurately computed for much larger range of inputs. Specifically, we compute
LogEI,- (x) = =1g__(((((x) - y²)/o(x)) + (x)(a(x))),
(8)
where log_h is mathematically equivalent to log h and can be stably and accurately computed by 1og[d(z)+ + +4)
ラ> -1 - -1/VE< <ive<==-1 (9) -1/V
h(z) - - hwaaaaooooeeeeeloooooooooo/////-
+ 02)
222/2-
1og(
where e (2r) and C₂ 12g(z/2)/2, /,, is the numerical precision, and log1mexp, erfcx are numerically stable implementations of log (1 - exp(2)) and exp(2" (rff(()))) respectively, see App. A Progenitors of Eq. 9 are found in SMAC 1.0 B5 and RoBO [46] which contain a log- transformed analytic EI implementation that is much improved, but can still exhibit instabilities as a grows negative, see App. Fig. I0 To remedy similar instabilities, we put forth the third, asymptotic case in Eq. (9). ensuring numerical stability throughout, see App. AA2 for details. The asymptotically quadratic behavior of log_h h becomes apparent in the last two cases, making the function particularly amenable to gradient-based optimization with significant practical implications for El-based BO.
4.2 Monte Carlo Parallel LogEl
Beyond analytic EI, Monte Carlo formulations of parallel El that perform differentiation on the level of MC samples, don't just exhibit numerically, but mathematically zero gradients for a significant proportion of practically relevant inputs. For qEI, the primary issue is the discrete maximum over the q outcomes for each MC sample in Gp. particular, the acquisition utility of expected improvement in Eq. 4 on single sample & of f is max,(5.(xx)) the acquisition utility in two stages: 1) May = soft softplus,, (E:(x)) y²) Rs (xx( - v²]+ and 2) ||u,lla /u/-lla/7mme z max, Mij- Notably, while we use canonical softplus and p-norm approximations here, specialized fat-tailed non-line non-linearities are required to scale to large batches, see Appendix [A.4 Since the resulting quantities are strictly positive, they can be transformed to log-space permitting an implementation of qLogEI that is numerically stable and can be optimized effectively. In particular,
y²]+ Mathematically, we smoothly approximate
qLogEI, y- (X) = =log| / (Σ)= softplus, (xxx) - 5²) 1/TmOXX
Tina
df
(10)
2i logsumexp, (Tmxc (omaxegg (logsoftplus, To (6'(x)) - -y0)) Tmax))
. 9 is the index for the candidate in the batch, and logsoftplus is a numerically stable implementation of log(log(1 exp(:))). See Appendix A.3 for details, including the novel fat-tailed non-lincarities like fatplus.
where i is the index of the Monte Carlo draws from the GP posterior, j L, 1.....
While the smoothing in (10) approximates the canonical qEI formulation, the following result shows that the associated relative approximation error can be quantified and bounded tightly as a function of the temperature parameters To, Tmax and the batch size q. See Appendix C for the proof.
Lemma 2. [Relative Approximation Guarantee] Given To, Trax 0, the approximation error of gLogEl to qEI is bounded by
(Logg))( qEI(X)] S Tm_ 1) QEI(X) + log(2) Toun.
(11)
we show the importance of setting the temperatures sufficiently low for qLogEl to achieve good optimization characteristics, something that only becomes possible by transforming all involved computations to log-space. Otherwise, the smooth approximation to the acquisition utility would exhibit vanishing gradients numerically, as the discrete max operator does mathematically.
In Appendix D.10

---

4.3 Constrained ET
Both analytic and Monte Carlo variants LogEl can be extended for optimization problems with black-box constraints. For analytic CEI with independent constraints of the form [i( S 0. the constrained formulation in Eq. (3) simplifies to LogCEI(x) (x) LogEl x) + Σ, log( (x) S 0)), which can be readily and stably computed using LogEl in Eq. (8) and, if Si is modelled by a GP, stable implementation of the Gaussian log cumulative distribution function. For the Monte Carlo variant, apply similar strategy as for Eq. (1O) to the constraint indicators in Eq. (3) 1) smooth approximation and 2) an accurate and stable implementation its log value, see Appendix
4.4 Monte Carlo Parallel LogEHVI
The numerical difficulties of qEHVI in 66 are similar those of qBI, and the basic ingredients of smoothing and log-transformation still apply, but the details are significantly more complex since qEHVI uses many operations that have mathematically zero gradients with respect to some of the inputs. Our implementation is based on the differentiable inclusion-exclusion formulation of the hypervolume improvement [3]. As a by-product, the implementation also readily allows for the differentiable computation of the expected log hypervolume, instead of the log expected hypervolume, note the order, which can be preferable in certain applications of multi-objective optimization (26).
5 Empirical Results
We compare standard versions of analytic EI (EI) and constrained EI (CEI), Monte Carlo parallel EI (qEI), as well as Monte Carlo EHVI (qEHVI) in addition to other state-of-the of-thc-art baselines like lower-bound Max- Value Entropy Search (GIBBON) 61 and single- and multi-objective Joint Entropy Search (JES) [36 71. All experiments are implemented using BoTorch [6] and utilize multi-start optimization of the AF with scipy's s L-BFGS-B optimizer. In order to avoid conflating the effect of BoTTorch's s default initialization strategy with those of our contributions, we use 16 initial points chosen uniformly at random from which to start the L-BFGS-B optimization. For a comparison with other initialization strategies, see Appendix D We run multiple replicates and report mean and error bars of +22 standard errors of the mean. Appendix D.] contains additional details.
Single-objective sequential BO
We compare EI and LogEl the 10-dimensional convex Sum- of-Squares (SoS) function fx) E (z; 0.5) using 20 restarts seeded from 1024 pseudo- random samples through Bo'Trrch's 's default initialization heuristic. Figure shows that due to vanishing gradients, El is unable make progress even on this trivial problem.
Regret
Aqquistion Eunction value
Acquiition dunction valle
10-4
EI
10
LoglEl
exp(LagED
Regret
10"
10-2
10"
10-4
50
100
125
150
175
200
50
100
125
175
200
Nomber of evallutions
Number of evaluations
Figure 2: Regret and El acquisition value for the candidates selected by maximizing El and LogEl on the convex Sum-of-Square problem. Optimization stalls out for El after about 75 observations due to vanishing gradients (indicated by the jagged behavior of the acquisition value), while LogEl continues to make steady progress.
In Figure 3] we compare performance on the Ackley and Michalewicz test functions 667. Notably, LogEl substantially outperforms El on Ackley as the dimensionality increases. Ackley is challeng- ing multimodal function for which it is critical to trade off local exploitation with global exploration, a task made exceedingly difficult by the numerically vanishing gradients of El in a large fraction of the search space. We see similar albeit less pronounced behavior on Michalewicz, which reflects the fact that Michalewicz is somewhat less challenging problem than Ackley.
7

---

2=
4=
16
d=
nodom
14
Boxt dbicirll value
Logut CIRBON
Jallyy
IES
10
10
100
150
0.8
dbsteed value
Michllnnat
-1.0-
-1.2
-5-
-1.6
-7
-1.8
100 1.50
150
200
150
200
250
Function ovaliation
Fuaction evaliations
Function ovaliaaions
Figure 3: Best objective value as a function of iterations on the moderately and severely non-convex Michalewicz and Ackley problems for varying numbers of input dimensions. LogEl substantially outperforms both and GIBBON, and this gap widens as the problem dimensionality increases. JES performs slightly better than LogEl on Ackley, but for some reason fails on Michalewicz. Notably, JES is almost two orders of magnitude slower than the other acquisition functions (see Appendix D)
Temninn-Compessaiie Compecsionnnnnn Swing
Vecual Design
Welded Design Nlack-hox
Redueer Design black-hus
with blcc--hi
0.25
33000
ezi
Iogoti
17000
Ingent
10
logget
5000
logoni
(outuu
0.15
kg
14000
kg
sibo
kg
4500
0.10
11000
6
0.05
3500
3000
5000
50
100
Berations
Itenutions
aterations
Figure 4: Best feasible objective value as a function of number of function evaluations (iterations) on four engineering design problems with black-box constraints after an initial 2d pseudo-random evaluations.
BO with Black Box Constraints 4 shows results on four engineering design problems with black box constraints that were also considered in (22). We apply the same bilog transform the trust region-based SCBO method [22) to all constraints to make them easier to model with a GP. We that LogCEI outperforms the naive CEI implementation and converges faster than SCBO. Similar to the unconstrained problems, the performance gains of LogCEI over CEI grow with increasing problem dimensionality and the number of constraints. Notably, we found that for some problems, LogCEI in fact improved upon some of the best results quoted in the original literature, while using three orders of magnitude fewer function evaluations, see Appendix D.7 for details.
Parallel Expected Improvement with gLogEI Figure reports the optimization performance of parallel BO on the 16-dimensional Ackley function for both sequential greedy and joint batch optimization using the fat-tailed non-linearitier of App. [A.4] In addition to the apparent advantages of ologE over qEI, a key finding is that jointly optimizing the candidates of batch acquisition functions can yield highly competitive optimization performance, see App. D.3 for extended results.
High- Hiigh-dimensional BO with gLogEI Figure 6 shows the performance of LogEl on three high- dimensional problems: the 6-dimensional Hartmann function embedded a 100-dimensional space, a 100-dimeniional rover trajectory planning problem, and 103-dimensional SVM hyperparameter tuning problem. We use a 103-dimensional version of the 388-dimensional SVM problem considered by Eriksson and Jankowiak (21) where the 100 most important features were selected using Xgboost.
8

---

16
150
desteberrva valus
12.5
quijet
qei soq
125
Ackryy ISD
10.0
glogei soq
glogri
10.0
7.5
gibbon
7,5
75
5.0
5.0
2.5
0
100
200
100
200
100
200
Tunction ovaluations
Function evaliations
Function evalustions
Figure 5: Best objective value for parallel BO function of the number evaluations for single- objective optimization on the 16- (6-dimensiona Ackley function with varying batch sizes Notably, joint optimization of the batch outperforms sequential greedy optimization.
-1.0
100D Fishedded Hartmannó (q
0.19
103D
(g=
100D Rover (g=
Best obsorved value
-1.5
0.18
3.5
-2.0
Tox RMSE
0.17
Reward
0.16
-2.5
0.15
-3.0
0.14
-3.5
25
75
0.13
25
75
100
2.0
50 Berations
75
100
Berations
glogei saas
qci
glogei
gei
random
Figure 6: Best objective value function of number of function evaluations (iterations) three high-dimensional problems, including Eriksson and Jankowiak [2]] 's SAAS prior.
Figure 6 shows that the optimization exhibits varying degrees of improvement from the inclusion of gLogEl, both when combined with SAASBO (2] and a standard GP. In particular, qLogEI leads to significant improvements on the embedded Hartmann problem, even leading BO with the canonical GP ultimately catch up with the SAAS MaA---ror--quuppee model. On the other hand, the differences on the SVM and Rover problems are not significant, see Section 6 for discussion.
Multi-Objecfio optimization with gLogEHVI Figure compares gLogEHVI and qEHVI two multi-objective test problems with varying batch sizes, including the real-wordd-inspired cell network design for optimizing coverage and capacity [9 The results are consistent with our findings in the single-objective and constrained cases: qLogEHVI consistently outperforms qEHVI and even JES [7T) for all batch sizes. Curiously, for the largest batch size and DTLZ2, gLogNEHVI's improvement over the reference point (HV > 0) occurs around three batches after the other methods, but dominates their performance in later batches. See Appendix for results on additional synthetic and real- world-inspinn multi-objective problems such as the laser plasma acceleration optimization
38 and vehicle design optimization [54 68
Discussion
To recap, El exhibits vanishing gradients 1) when high objective values are highly concentrated the search space, and 2) as the optimization progresses. In this section, we highlight that these conditions are not met for BO applications, and that LogEl's performance depends on the surrogate's 's quality.
On problem dimensionality While our experimental results show that advantages of LogEI generally grow larger as the dimensionality of the problem grows, we stress that this fundamentally due to the concentration of high objective values in the search space, not the dimensionality itself. low intrinsic dimensionality,
Indeed, we have observed problems with high ambient dimensionality where LogEl does not lead to significant improvements over EI, e.g. the SVM problem in Figure

---

0.06
0.06
Call Call Network 30D
Best observed HV
0.05
0.05
0.04
0.04
0.04
0.03
0.03
0.03
0
100
200
100
200
100
200
0.4
DTL 16D
observed
0.3
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0.1
0.0
0.0
0
100
200
100
200
100
200
Function evaluations
Function evalustions
Function evaluations
Rand
SNERVI
qLogNERVI
JES
Figure 7: Batch optimization performance on two multi-objective problems, as measured by the hypervolume of the Pareto frontier across observed points. This plot includes JES [71). Similar the single-objective case, the LogEl variant qLogEHVI significantly outperforms the baselines.
On asymptotic improvements While members of the LogEI family can generally be optimized better, leading to higher acquisition values, improvements in optimization performance might be small in magnitude, e.g. the log-objective results on the convex 10D sum of squares Fig. 2 or only begin to materialize in later iterations, like for y = 16 on DTLZ2 in Figure
On model quality Even if good objective values are concentrated small volume of the search space and many iterations are run, LogEl might still not outperform EI if the surrogate's ' predictions are poor, or its uncertainties are not indicative the surrogate's mismatch the objective, see Rover in Fig. 6 In these cases, better acquisition values do not necessarily lead to better BO performance.
Replacing EI Despite these limitation, we strongly suggest replacing variants of EI with their LogEI counterparts. If LogEI were dominated by El on some problem, it would be an indication that the EI family itself is sub-optimal, and improvements performance can be attributed to the exploratory quality of randomly distributed candidates, which could be incorporated explicitly.
7 Conclusion
Our results demonstrate that the problem vanishing gradients is major source of the difficulty of optimizing improvement-basee careful reformulations and implementations. As a result, we see substantially improved optimization performance across a variety of modified EI variants across broad range of problems. In particular, we demonstrate that joint batch optimization for parallel BO can be competitive with, and times exceed the sequential greedy approach typically used in practice, which also benefits from our modifications. Besides the convincing performance improvements, one of the key advantages of our modified acquisition functions is that they are much less dependent on heuristic and potentially brittle initialization strategics. Moreover, our proposed modifications do not meaningfully increase the computational complexity of the respective original acquisition function.
acquisition functions and that we can mitigate this issue through
While our contributions may not apply verbatim to other classes of acquisition functions, our key insights and strategies do translate and could help with e.g. improving information-basee 34 76), cost-aware [5]] 66), and other types of acquisition functions that are prone to similar numerical challenges. Further, combining the proposed methods with gradient-awar first-order BO methods (5 16 23 could lead to particularly effective high-dimensional applications BO, since the advantages of both methods tend to increase with the dimensionality of the search space. Overall, we hope that our findings will increase awareness in the community for the importance of optimizing acquisition functions well, and in particular, for the care that the involved numerics demand.
10

---

Acknowledgments and Disclosure of Funding
The authors thank Frank Hutter for valuable references about prior work on numerically stable compu- tations of analytie EI, David Bindel for insightful conversations about the difficulty of optimizing EI, as well as the anonymous reviewers for their knowledgeable feedback.
References
[1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Goodfellow, Andrew Harp. Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens. Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. TensorFlow: Large- scale machine learning on heterogeneous systems, 2015. URL http:::t///// wttp:///wwwtttw..tllleeee tensorflos org/
[2] Sebastian Ament and Michael O'Neil. Accurate and efficient numerical calculation of stable densities via optimized quadrature and asymptotics. Statistics and Computing. 28:171-185. 2018.
[3] Sebastian Ament, Maximilian Amsler, Duncan R. Sutherland, Ming-Chhang Chang, Dan Gue- varra, Aine B. Connolly, John M. Gregoire, Michael O. Thompson, Carla P. Gomes, and R. Bruce van Dover. Autonomous materials synthesis via hierarchical active learning nonequilibrium phase diagrams. Science Advances, 7(51) (ebbg4930, 2021. doi: 10. 1266 URL wwt.:1112/ssccaa http:/wwwwwwcceee////////.e ort/ abgp abg4930.
[4] Sebastian Ament, Andrew Witte, Nishant Garg, and Julius Kusuma. Sustainable concrete via bayesian optimization, 2023. URL httts::tp/rsrst .ww.vv.or///bbb. arg/abs/23111 22100.231...00 18288 NeurIPS 2023 Workshop on Adaptive Experimentation in the Real World.
[5] Sebastian E Ament and Carla P Gomes. Scalable first-order Bayesian optimization via structured automatic differentiation. In Kamalika Chaudhuri, Stefanie Jegelka. Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference Machine Learning, volume 162 Proceedings of Machine Learning Research, pages 500-516. URL https: http://///oceeeeessseeee /
PMLR, 17-23 Jul 2022.
mlr www.srrpreee/eee v16s////vv1/// ame/t22a.......m. html
[6] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton. Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch: Framework for Efficient Monte-Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems 33, 2020.
[7] Ricardo Baptista and Matthias Poloczek. Bayesian optimization of combinatorial structures, 2018.
[8] Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Max-value entropy search for multi-objective bayesian optimization with constraints, 2020.
[9] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Miln and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http://giie htttp://ggthtbb.ccom/ogg/////// cow.google//ax.
[10] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16(5): 190-1208, 1995.
[11] Carlos A Coello Coello and Efrén Mezura Montes. Constraint-handling in genetic algorithms through the use dominance-basee tournament selection. Advanced Engineering Informatics, 16(3):193-28) 193-20,,,, 2002.

---

[12] Alexander I. Cowen-Rivers, Wenlong Lyu. Rasul Tutunov, Zhi Wang, Antoine Grosmit, Ryan Rhys Griffiths, Alexandre Max Maraval, Hao Jianye, Jun Wang, Jan Peters, and Haitham Bou Ammar. Hebo pushing the limits of sample-efficient hyperparameter optimisation, 2022.
[13] Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hypervolume improvement for parallel multi-objective bayesian optimization. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 9851-9864. Curran As-
sociates, Inc., 2020. URL httops://proceeddngg: heg.reerrsnnns...u.pp.s.s.. neurips. wtt.riiss.cppppppeeeee httpprrpp...ccpppp///////iiaeee. pappp//00//ttt/le// .pf
[14] Samuel Daulton, Sait Cakmak, Maximilian Balandat, Michael A. Osborn,, Enlu Zhou, and Eytan Bakshy. Robust multi-objective Bayesian optimization under input noise. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, edi- tors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 4831-4866. PMIR, 17-23 Jul 2022. URL https:://proceediiin
mlr press/ ppress/vl8e 2rress//122/daultoo22aa html
[15] Samuel Daulton, Xingchen Wan, David Eriksson, Maximilian Balandat, Michael A. Osborne, and Eytan Bakshy. Bayesian optimization over discrete and mixed spaces via probabilistic reparameterization. In Advances in Neural Information Processing Systems 35, 2022.
[16] Filip De Roos, Alexandra Gessner, and Philipp Hennig. High-dimensional gaussian process inference with derivatives. In International Conference on Machine Learning, pages 2535-2545. PMLR, 2021.
[17] Kalyan Deb, L. Thiele, Marco Laumanns, and Eckart Zitzler. Scalable multi-objective op- timization test problems. volume 1, pages 825-830, 2002. ISBN 0-7803-7282-4. doi: 10.1 10.1109/CEEC.20022007777
[18] Aryan Deshwal, Sebastian Ament, Maximilian Balandat, Eytan Bakshy, Janardhan Rao Doppa, and David Eriksson. Bayesian optimization over high-dimensional combinatorial spaces via dictionary-basee embeddings. In Francisco Ruiz, Jennifer Dy, and Jan- On-Willem van de Meent, editors, Proceedings of The 26th International Conference Artificial Intelligence and Statis- tics, volume 206 of Proceedings of Machine Learning Research, pages 7021-7039, PMLR, 25-27 Apr 2023.
[19] Ryan M. Dreifuerst, Samuel Daulton, Yuchen Qian, Paul Varkey, Maximilian Balandat, Sanjay Kasturia, Anoop Tomar, Ali Yazdan, Vish Ponnampalam, and Robert W. Heath. Optimizing coverage and capacity in cellular networks using machine learning, 2021.
[20] M. T. M. Emmerich, K. C. Giannakoglou, and B. Naujoks. Single- and multiobjective evo- lutionary optimization assisted by gaussian random field metamodels. IEEE Transactions Evolutionary Computation, 10(4):421-43 2006.
[21] David Eriksson and Martin Jankowiak. High-dimensional Bayesian optimization with sparse axis-aligned subspaces. In Uncertainty in Artificial Intelligence. PMLR, 2021.
[22] David Eriksson and Matthias Poloczek. Scalable constrained Bayesian optimization. In International Conference Artificial Intelligence and Statistics. PMLR, 2021.
[23] David Eriksson, Kun Dong, Eric Lee, David Bindel, and Andrew Wilson. Scaling gaussian process regression with derivatives. Advances in neural information processing systems, 31, 2018.
[24] David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems 32, NeurIPS, 2019.
[25] Peter Frazier. A tutorial bayesian optimization. arXiv preprint arXiv:1807.022 100.::0107.02111.02000 2018.
12

---

[26] Tobias Friedrich, Karl Bringmann, Thomas VoB, and Christian Igel. The logarithmic hypervol- umee indicator. In Proceedings of the 11th workshop proceedings on Foundations of genetic algorithms, pages 81-92, 2011.
[27] Jacob Gardner, Matt Kusner, Zhixiang, Kilian Weinberger, and John Cunningham. Bayesian optimization with inequality constraints. In Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pages 937-945, Beijing, China, 22-24 Jun 2014. PMLR.
[28] Roman Garnett. Bayesian Oplimization. Cambridge University Press, 2023, to appear.
[29] Michael A. Gelbart, Jasper Snoek, and Ryan P. Adams. Bayesian optimization with unknown constraints. In Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence, UAI, 2014.
[30] David Ginsbourger, Rodolphe Le Riche, and Laurent Carraro. Multi-points Criterion for Deterministic Parallel Global Optimization based Gaussian Processes. Technical report, March 2008. Uttpttts:://aal htt.atppp//aall.seemmeeaaaaaa-ap
[31] Robert B Gramacy, Annie Sauer, and Nathan Wycoff. Triangulation candidates for bayesian optimization. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 35933-35945. Curran Associates, Inc., 2022.
[32] 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):1-18, 2003. doi: 10.1162/163366632111887777
[33] Tim Head, Manoj Kumar, Holger Nahrstaedt, Gilles Louppe, and Iaroslav Shcherbatyi. aciktt-opimire/sckktt--ppiiiie October 2021. URL http::://o/:oooorr//////// wtt.://..oorg00..2//eeeoood 5565057
[34] José Miguel Henrández-Lobato, Matthew W. Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. In Proceedings of the 27th International Conference on Neural Information Processing Systems Volume 1, NIPS' 14, pages 918 926, Cambridge, MA, USA, 2014. MIT Press.
[35] Frank Hutter, Holger H Hoos, and Kevin Leyton- Brown. Sequential model-based optimization for general algorithm configuration. In Learning and Intelligent Optimization: 5th International Conference, LION 5, Rome, Italy, January 17-21, 2011. Selected Papers 5, pages 507-523. Springer, 2011
[36] Carl Hvarfner, Frank Hutter, and Luigi Nardi. Joint entropy search for maximally-informed bayesian optimization. In Advances in Neural Information Processing Systems 35, 2022.
[37] Wolfram Research, Inc. Wolfram alpha, 2023. httpss/tt:///ww http::////wwwwoofrraaaaaalaa..a com/
[38] F. Irshad, S. Karsch, and A. Döpp. Multi-objective and multi-fidelity bayesian optimization of laser-plasma acceleration. Phys. Rev. Res., 5:01 Jan 2023. doi: 10. R0.1100//hyyReeeRsseachhh 5.013063. URI https:////ik aps
51300663
[39] Faran Irshad, Stefan Karsch, and Andreas Doepp. Reference dataset of multi-objective and multi- fidelity optimization in laser-plasma acceleration, January 2023. URI.https::
org/10. 10000521//ze0 2281/zenodo.. 7565882
[40] Shali Jiang, Daniel Jiang, Maximilian Balandat, Brian Karrer, Jacob Gardner, and Roman Gar- nett. Efficient nonmyopic bayesian optimization via one-shot multi-step trees. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balean, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 8039-188449 Curran Associates, Inc., 2020.
[41] Donald Jones, C. Perttunen, and B. Stuckman. Lipschitzian optimisation without the lipschitz constant. Journal of Optimization Theory and Applications, 79:157-181, 1993. doi: 10.1007/BFO0941199.
13

---

[42] Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13:455 492, 1998.
[43] Kirthevasan Kandasamy, Karun Raju Vysyaraju, Willie Neiswanger, Biswajit Paria, Christo- pher R. Collins, Jeff Schneider, Barnabás Póczos, and Eric P. Xing. Tuning hyperparameters without grad students: Scalable and robust bayesian optimisation with dragonfly. Mach. Learn. Res., 21(1), jan 2020.
[44] Jungtaek Kim, Seungjin Choi, and Minsu Cho. Combinatorial bayesian optimization with random mapping functions to convex polytopes. In James Cussens and Kun Zhang, editors, Proceedings of the Thirty-Lighth Conference on Uncertainty in Artificial Intelligence, volume 180 of Proceedings of Machine Learning Research, pages 1001-1011. I. PMLR, 01-05 Aug 2022.
[45] Diederik P Kingma and Max Welling. Auto-Encoding Variational Bayes. arXiv e-prints, page arXiv:1312.6114 Dec 2013.
[46] A. Klein, S. Falkner, N. Mansur, and F. Hutter. Robo: A flexible and robust bayesian opti- mization framework in python. In NIPS 2017 Bayesian Optimization Workshop, December 2017.
[47] Dieter Kraft. A software package for sequential quadratic programming. Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt, 1988.
[48] Rémi Lam, Matthias Poloczek, Peter Frazier, and Karen Willcox. Advances bayesian optimization with applications in aerospace engineering. In 2018 AIAA Non-Deterministic Approaches Conference, page 1656, 2018.
[49] Robert Langer and David Tirrll. Designing materials for biology and medicine. Nature, 428, 04 2004.
[50] Armin Lederer, Jonas Umlauft, and Sandra Hirche. Posterior variance analysis of gaussian processes with application to average learning curves. arXiv preprint arXiv: arXiv:1966.00 1906.014064 2019.
[51] Eric Hans Lee, Valerio Perrone, Cedric Archambeau, and Matthiss Seeger. Cost-aware Bayesian Optimization. arXiv e-prints, page arXiv:2003.10870 March 2020.
[52] Benjamin Letham, Brian Karrer, Guilherme Ottoni, and Eytan Bakshy. Constrained bayesian optimization with noisy experiments. Bayesian Analysis, 14(2):495 -199, 06 2019. doi: 10. 1214/18-BA110..
[53] Qiaohao Liang, Aldair F. Gongora, Zekun Ren, Armi Tihonen, Zhe Liu, Shijing Sun, James R. Deneault, Daniil Bash, Flore Mekki-Berrada, Saif A. Khan, Kedar Hippalgaonkar, Benji Maruyama, Keith A. Brown, John Fisher III, and Tonio Buonassisi. Benchmarking the perfor- mance of bayesian optimization across multiple experimental materials science domains. npj Computational Materials, 7(1):188, 2021.
[54] Xingtao Liao, Qing Li, Xujing Yang. Weigang Zhang, and Wei Li. Multiobjective optimization for crash safety design of vehicles using stepwise regression model. Structural and Multidisci- plinary Optimization, 35:561-569, 06 2008. doi:
[55] Wenlong Lyu, Fan Yang. Changhao Yan, Dian Zhou, and Xuan Zeng. Batch Bayesian optimiza- tion via oblti-objective acquisition ensemble for automated analog circuit design. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 3306-3314 PMLR, 10-15 Jul 2018. URL MRRL: mlr. wta.r.press//eee v80/1yu18a. html
[56] Martin Mächler. Accurately computing log (1- exp (-1 al)) assessed by the rmpfr package. Technical report, Technical report, 2012.
[57] Mark McLeod, Stephen Roberts, and Michael A. Osbore. Optimization, fast and slow: optimally switching between local and Bayesian optimization. In Jennifer Dy and Andreas Krause, editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 3443-3452, PMLR, 10-15 Jul 2018.
14

---

[58] Jonas Močkus. On bayesian methods for seeking the extremum. In Optimization Techniques IFIP Technical Conference: Novosibirsk, July 1-7, 1974, pages 400-404. Springer, 1975.
[59] Jonas Mockus. The application of bayesian methods for seeking the extremum. Towards global optimization, 129 1978.
[60] Hossein Mohammadi, Rodolphe Le Riche, and Eric Touboul. Making ego and ema-es comple- mentary for global optimization, In Clarisse Dhaenens, Laetitia Jourdan, and Marmion, editors, Learning and Intelligent Optimization, pages 287-292, Cham, 2015. Springer International Publishing.
[61] Henry B. Moss, David S. Leslie, Javier González, and Paul Rayson. Gibbon: General-purpose information-basee bayesian optimisation. J. Mach. Learn. Res., 22(1), jan 2021.
[62] Changyong Oh, Jakub Tomezak, Efstratios Gavves, and Max Welling. Combinatorial bayesian optimization using the graph cartesian product. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 2914-2924. Curran Associates, Inc., 2019.
[63] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf. Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high- performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[64] Victor Picheny, Joel Berkeley, Henry B. Moss, Hrvoje Stojic, Uri Granta, Sebastian W. Ober, Artem Artemev, Khurram Ghani, Alexander Goodall, Andrei Paleyes, Sattar Vakili, Sergio Pascual-Diaz, Stratis Markou, Jixiang Qing, Nasrulloh R. B. S Loka, and Ivo Couckuyt. Trieste: Efficiently exploring the depths of black-box functions with tensorflow, 2023. URL https
//arxiv. abs /202.. 08436
[65] Carl Edward Rasmussen. Gaussian Processes in Machine Learning, pages 63-71. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.
[66] Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. In Advances in neural information processing systems, pages 2951-2959, 2012.
[67] Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved May 14, 2023, from http::/www..ttt sfu. ca/"ssur jano
[68] Ryoji Tanabe and Hisao Ishibuchi. An easy-to-use real-world multi-objective optimization problem suite. Applied Soft Computing, 89:106078, 2020. ISSN 1568-4946.
[69] The GPyOpt authors. github. com
GPyOpt: A bayesian optimization framework in python. http://
2016.
[70] Aimo Töm and Antanas Zilinskas. Global optimization, volume 350. Springer, 1989.
[71] Ben Tu, Axel Gandy, Nikolas Kantas, and Behrang Shafei. Joint entropy search for multi- objective bayesian optimization. Advances Neural Information Processing Systems 35, 2022.
[72] Ryan Turner, David Eriksson, Michael McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu, and Isabelle Guyon. Bayesian optimization is superior to random search for machine learning hyperparameter tuning: Analysis of the black-box optimization challenge 2020. In NeurIPS 2020 Competition and Demonstration Track, 2021.
[73] Andreas Wächter and I Lorenz T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106 (1):25-57, 2006.
15

---

[74] Xingchen Wan, Vu Nguyen, Huong Ha, Binxin Ru, Cong Lu, and Michael A. Osborn.. Think global and local: Bayesian optimisation over high- -ighesimnssooall categorical and mixed search spaces. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 10663-10674. PMLR, 18-24 Jul 2021.
[75] Jialei Wang. Scott Clark, Eric Liu, and Peter I. Frazier. Parallel bayesian global optimization of expensive functions, 2016.
[76] Zi Wang and Stefanie Jegelka. Max-value Entropy Search for Efficient Bayesian Optimization. ArXiv e-prints, page arXiv: 770.001968, March 2017.
[77] James Wilson, Frank Hutter, and Mare Deisenroth. Maximizing acquisition functions for bayesian optimization. In Advances in Neural Information Processing Systems 31, pages 9906-9916.. 2018.
[78] Jian Wu and Peter I. Frazier. optimization. In Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS' 16, page 3134-3142. Curran Associates Inc., 2016.
The parallel knowledge gradient method for batch bayesian
[79] Jian Wu, Matthias Poloczek, Andrew Gordon Wilson, and Peter I Frazier. Bayesian optimization with gradients. Advances in Neural Information Processing Systems, pages 5267-5278, 2017.
[80] Kaifeng Yang. Michael Emmerich, André Deutz, and Thomas Bäck. Multi-objective bayesian global optimization using expected hypervolume improvement gradient. Swarm and Evolution- ary Computation, 44:945 956, 2019.
[81] Dawei Zhan and Huanlai Xing. Expected improvement for expensive optimization: a review. Journal of Global Optimization, 78(3):507-544, 2020.
[82] Jiayu Zhao, Renyu Yang, Shenghao Qiu, and Zheng Wang. Enhancing high-dimensional bayesian optimization by optimizing the acquisition function maximizer initialization. arXiv preprint arXiv:2302.08298 2023.
[83] Eckart Zitzler, Kalyanmoy Deb, and Lothar Thiele. Comparison multiobjective evolutionary algorithms: Empirical results. Evol. Comput., 8(2):133-1955 jun 2000. ISSN 1063-6560. doi: 10.1 10.11622106666600556222. URL http:://.....oooooo......6666660
16

---

A Acquisition Function Details
A.1 Analytic Expected Improvement
Recall that the main challenge with computing analytic LogEl to accurately compute log h, where
h(2) = (== + (2). with ((>) exp(-2²/2)/ and ((2) J Φ(u)du. To express log
numerically stable form becomes increasingly negative, we first take the log and multiply
out of the argument to the logarithm:
((z) (>z)
log l(z) - log log 1
(12)
Fortunately, this form exposes the quadratic factor, ((=)/0(z) can be computed via standard imple- mentations of the scaled complementary error function erfcx, and log (1 =#<=)/((2)) the last term of Eq. (12), can be computed stably with the loglmexp implementation proposed in [56]:
[xo(-=xpp1(x)) - log2 <
log1mexp(r) 1oglp(--eppl - log 2 r
(13)
where expm1, logip are canonical stable implementations of exp(x)) and log(1 + x)) respectively. In particular, we compute
Log(o(z) + ==))))
2> =>-1
log_h(2) - /2 e + loggmxxppllgggrrtttt==/00/0000
(14)
- 22 - 210g(12))
2 2 -1/vE.
where C1 (0rg(2r)//3)))) 0z = 1og(#/2)/2. and e is the floating point precision. We detail the reasoning behind the third asymptotic case of Equation (14) in Lemma of Section A.2 below.
Numerical Study of Acquisition Values Figure& shows both the numerical failure mode of naïve implementation of EI, which becomes exactly zero numerically for moderately small 2, while the evaluation via log_h in Eq. (14) exhibits quadratic asymptotic behavior that is particularly amenable to numerical optimization routines.
Expected Improvement
LogEl
Naivo
Fiälure
4000
-lg(EE))
3000
2000
1000
-100
-50
Figure 8: Plot of the log computed via log h and log_h in Eq. (14). Crucially, the naïve implementation fails as 2 (j(x) f (x)) becomes increasingly negative, due being exactly numerically zero, while our proposed implementation exhibits quadratic asymptotic behavior.
Equivalence of Optimizers Lemma states that if the maximum of El is greater than 0, LogEl and EI have the same set of maximizers. Furthermore, if max El(x) O, then x= arg max El(x). In this case, LogEI is undefined everywhere, so it has no maximizers, which we would yield the same BO policy EI, for which every point is maximizer. Lemma 3. If maxxxx EI(x) >
0, then arg maxxex EI(x) arg maxxxx,EIxaa (x)).(x).
Proof. Suppose maXxex [I 0. Then arg max EI(x) arg max x [I)
For all x X such that El(x) 0. LogEI(x) log (x)) Since log monontonic, we have that arg max 2 = arg max lzg(2). Hence, arg maxxex,B())>0 El(x)
arg max Lx).
17

---

Difference of Laurent Expansion to logimexp Implementation
10a
10-1
D(1/222)
0(1/2%4)
Order
Absolute Difference
Order
10
Order
Order 10
10
Order 12
Order 14
Order
10-10
Order 20
Order 22
Order 24
Order
10
Order 28
Order 30
Order 32
20.0
17.5
15.0
12.5
10.0
7.5
5.0
2.5
Figure 9: Convergence behavior of the asymptotic Laurent expansion Eq. (66) different orders.
A.2 Analytical LogEl's Asymptotics
As 2 grows negative and large, even the more robust second branch in Eq. as well as the implementation of Hutter et al. [35] and Klein et al. [46] can suffer from numerical instabilities (Fig. [0] left). In our case, the computation negative 2. For this reason, we propose an approximate asymptotic computation based on Laurent expansion at -c. As result of the full analysis in the following, we also attain a particularly simple formula with inverse quadratic convergence in 2, which is the basis of the third branch of Eq. (14)
the last term of Eq. (2) is problematic for large
*0(=)
log
[+
((z) 2 log(l=)) + O(I== 91).
(15)
In full generality, the asymptotic behavior of the last term can be characterized by the following result.
Lemma 4 (Asymptotic Expansion). Let < -1 and K N, then
log (1+ =0(=) = log
K
h-1
Σ
(-1
II 1=0
(2) + 1)
2 -2k
2k
+ 10(1==2(K-
(16)
((2)
Proof. We first derived Laurent expansion of the non-log-trane transformed =0(2)/6(2), key quantity in the last term of Eq. (12), with the help of Wolfram Alpha (37]:
(x)) ((2)
Σ
II ===
(2) + 1)
24 +((==-22K)
(17)
remains to derive the asymptotic error bound through the log log-transformation of the above expansion. Letting L(z, K)
- K
(-1)
1 (2j + +1]
we get
2x)
log
log (L(z, K) (12)122)
o(2)
log L(, K) + log(1 (1 L((x,K)
(18)
log ((((1=) 0 =12_4(:;40-0000012222//1111111 LL2,,
logg [o,
The penultimate equality is due to log(1+x) + O (2)), the last due to L(z, K) O(|2|
18

---

+P
Comparison with RoBO LogEi
Difference HEBO's Approximation LogE!
Lage
RoBO Loge:
02
30"4
-0.4
30
10gF00)
-08
Diffrrmnee Loge
10"
10
-08
10
-12
10:45
(Canecled HERO LogE))-
-6000
-3000 -2000 -1000
-20.0
17.5
-15.0
-125
-30.0
7.5
-5.0
-25
Left: Comparison of LogEl values in single-precision floating point arithmetic, function of = (μ(x) (z) to RoBO's LogEl implementation [46]. Notably, RoBO's LogEl improves greatly the naïve implementation (Fig. 8). but still exhibits failure points (red) well above floating point underflow. The implementation of Ei.. (14) continues to be stable in this regime. Right: Comparison of LogEl values to HEBO's approximate LogEl implementation 12 At the time of writing, HEBO's implementation exhibits discontinuity and error of > 1.02 at the threshold -6, below which the approximation takes effect. The discontinuity could be ameliorated, though not removed, by correcting the log normalization constant (turquoise). The figure also shows that the naive implementation used by HEBO for z -6 starts to become unstable well before --6
Figure 10:
SMAC 1.0 and RoBO's Analytic LogEI To our knowledge, SMAC 1.0's implementation of the logarithm of analytic EI due Hutter et al. [35), later translated RoBO (46), was the first to improve the numerical stability of analytic El through careful numerics. The associated imple- mentation mathematically identical to log EI, and greatly improves the numerical stability of the computation. For large negative however, the implementation still exhibits instabilities that gives rise to floating point infinities through which useful gradients cannot be propagated (Fig. IO left). The implementation proposed herein remedies this problem by switching to the asymptotic approximation of Eq. (15) once it accurate to machine precision e. This is similar to the use of asymptotic expansions for the computation of a-stable densities proposed by Ament and O'Neil (2).
HEBO's Approximate Analytic LogEI HEBO 12 contains an approximation to the logarithm
of analytical EI as part of its implementation of the MACE acquisition function (55), which the time of writing - is missing the log normalization constant of the Gaussian density, leading to a large discontinuity at the chosen cut-off point of -6 below which the approximation takes effect, see here for the HEBO implementation. Notably, HEBO does not implement an exact stable computation of LogEI like the one of Hutter et al. (35) Klein et 46) and the non-asymptotic branches of the current work. Instead, it applies the approximation for all x < -6, where the implementation exhibits maximum error of 1.02, if the implementation' 's normalization constant were corrected, maximum error of > 0.1. By comparison, the implementation put forth herein is mathematically exact for the con-asymptotic regime 2 > 11/ /C and accurate numerical precision in the asymptotic regime due to the design of the threshold value.
A.3 Monte-Carlo Expected Improvement
For Monte-Carke cannot directly apply similar numerical improvements for the analytical version, because the utility values, the integrand of Eq. (4), on the sample level are likely to be mathematically zero. For this reason, we first smoothly approximate the acquisition utility and subsequently apply log transformations to the approximate acquisition function.
To this end, natural choice is softplus (z) To log(1 + exp(2/To) for smoothing the max(0, 2). where To a temperature parameter governing the approximation error. Further, we approximate the max; over the y candidates by the norm || | and note that the approximation error introduced by both smooth approximations can be bound tightly as a function of two "temperature" parameters To and Tmaxs see Lemma
19

---

Importantly, the smoothing alone only solves the problem of having mathematically zero gradients, that of having numerically vanishing gradients, we have shown for the analytical case above. For this reason, we transform all smoothed computations to log space and thus need the following special implementation of log o softplus that can be evaluated stably for a very large range of inputs:
[1og softplus,](a
z/T > l x E/T
logsoftplus, (1)
x//T log(r)
where T is a temperature parameter and / depends on the floating point precision used, around -35 for double precision in our implementation.
Note that the lower branch of logsoftplus approximate. Using a Taylor expansion (og(1+2) O,
see that the approximation error is O(=²), and therefore, Log(log((1 + exp(x))) = æ + (exp(xx)² which converges to a exponentially quickly as x + -00 In implementation, / is chosen that significant digit is lost in dropping the second order term from the lower branch.
/ 2/2+ + (=²²) around
Having defined logsoftplus, we further note that
log | 1/Tmae = log Σ
= Tmax log
Σ
exp((ogg(x))
Tmax Tonaxoogsmexpp [lo((:))/ mmux))
Therefore, we express the logarithm of the smoothed acquisition utility for q candidates as
Tmax 1 logsumexpppppp (logsoftplus (mans//gruurrrryoogototttop ((5' (I5) - y²) Tmax)-
Applying another logsumexp to compute the logarithm of the mean of acquisition utilities over set of Monte Carlo samples {Erhe gives rise to the expression Eq. (1O).
In particular for large batches (large q). this expression can still give rise to vanishing gradients for some candidates, which is due the large dynamic range of the outputs of the logsof optlus when << 0. To solve this problem, we propose a new class of smooth approximations to the "hard" non-linearities that decay as (11/z²) as I + -0 in the next section.
A.4 Class of Smooth Approximations with Fat Tails for Larger Batches
regular softplus(z)) = log(1 + exp(z)) function smoothly approximates the ReLU non-line con-lincaritty and - in conjunction with the log transformations is sufficient to achieve good numerical behavior for small batches of the Monte Carlo acquisition functions. However, more candidates are added, log softplus(x) = log(log(1 + exp(x))) is increasingly likely to have a high dynamic range as for > is chosen to be small, (-z/r) can vary orders of magnitude within a single batch. This becomes problematic when we approximate the maximum utility over the batch of candidates, since logsumexp only propagates numerically non-zero gradients to inputs that are no smaller than approximately (max, xj 700) in double precision, another source of vanishing gradients.
< log softplus, (x) -z/T.
To solve this problem, we propose a new smooth approximation to the ReI.U, maximum, and indicator x - -o,, instead of exponentially, like the canonical softplus. The high level idea is to use (1 + - which is proportional to the Cauchy density function (and is also known Lorentzian), ways that maintain key properties of existing smooth approximations convexity, positivity, etc while changing the asymptotic behavior of the functions from exponential to O(1/z²) as + -c, also known "fat tail". Further, we will show that the proposed smooth approximations satisfy similar maximum error bounds as their exponentially decaying counterparts, thereby permitting a similar approximation guarantee as Lemma 2 with minor adjusments to the involved constants. While the derivations herein are based on the Cauchy density with inverse quadratic decay, possible to generalize the derivations to e.g. a-stable distribution whose symmetric variants permit accurate and efficient numerical computation [2]
functions that decay only polynomially
20

---

Fat Softplus
We define
8 (z) = a11 + + log(1 + exp(x)),
(19)
for positive scalar a. The following result shows that we can ensure the monotonicity and convexity both important properties of the ReLU that we would like to maintain in our approximation of by carefully choosing
1.0
farckus
fanikss
softs us
0,8
-20
0.6
0.4
-60
0.2
-100
0.0
Figure 11: The fat softplus approximates Figure 12: The fat softplus has an (1/=²² asymp- max(2, 0) similarly tightly as the regular softplus totic decay, versus the (exx(x)) decay as x + and is also monotonic, convex, and positive. The -oc, moderating the dynamic range of the quanti- ties involved in parallel LogEl.
plot used temperature of To 0.01.
Lemma (Monotonicity and Convexity). 40+ (x) is positive, monotonically increasing, and strictly convex for satisfying
et/v3
0 0<a <
2 (1+
Proof. Positivity follows due to a 2 0, and both sumands being positive. Monotonicity and convexity be shown via canonical differential calculus and bounding relevant quantities.
particular, regarding monotonicity, we want to select a so that the first derivative bounded below by zero:
ex
2r
(x)+(x) = 1 C²
0
(1+ (1+x
First, note that 0g4p+(2) is positive for and any a, since both terms positive in this
regime. For 2 Ite= (1 -)- 1/2, and -1/(1 + 22)² 2 -1/(1 x²), that
2x
0, (x) 2 2
a (1 - x²)
Forcing 2 at "1++27 O, and multiplying by (1 + x gives rise to quadratic equation whose roots
1
2x
2a + / 2 Thus. there real roots for 1/2. Since the derivative is certainly positive for the negative reals and the guaranteed non-existence of roots implies that the derivative cannot eross zero elsewhere, 0 < 1/2 is sufficient condition for monotonicity of ++
Regarding convexity, our goal to prove a similar condition on a that guarantees the positivity of the second derivative:
6x + (1+e-03
(x (1 + x233
Note that 6-2-2 (I+r") is symmetric around 0, is negative in (- 1/3, v1/3) and has a minimum of -2 at 0. [1+e Zyz is symmetric around zero and decreasing away from zero. Since the rational
11/8
polynomial is only negative in (- (-31/3, /1/3), can lower bound u+
in
(1+e 1/73
21

---

(- (-v1/3, /1/3). Therefore,
(x)+
(1+ 2a
V173 (1+e -14e-///aaa> /2 0.115135,
Forcing (1+e- ///y
2a > 0 and rearranging yields the result. Since
the convexity condition is stronger than the monotonicity condition and therefore subsumes it.
Importantly yp decays only polynomially for increasingly negative inputs, and therefore log yp only log- arithmically, which keeps the range of oo constrained to values that are more manageable numerically, Similar to Lemma one can show that
(x/T(x/)) ReLU(>)] & (a + log(2)) T.
(20)
There are large number of approximations variants of the ReLU that have been proposed as activation functions of artificial neural networks, but to our knowledge, none satisfy the properties that we seek here: (I) smoothness, (2) positivity, (3) monotonicity, (4) convexity, and (5) polynomial decay. For example, the leaky ReLU does not satisfy (1) and (2), and the ELU does not satisfy (5).
Fat Maximum The canonical logsumexp approximation max; I; suffers from numerically vanishing gradients max; 22 - min, z is larger a moderate threshold, around 760 in double precision, depending the floating point implementation. In particular, while elements close the maximum receive numerically non-zero gradients, elements far away are increasingly likely to have a numerically zero gradient. To fix this behavior for the smooth maximum approximation,
propose
I1+
(*t
T4 - maxy Ti ( )
>
Pmax (x) = max zj + T log Σ 1+
(21)
This approximation to the maximum has the same error bound to the true maximum as the logsumexp approximation:
Lemma 6. Given T > 0
max F. ST T Tmmax((/// (x/T) & max v; + + log(d).
(22)
Proof. Regarding the lower bound, let = arg max, zy. For this index, the associated summand in 21 is 1. Since all sumands are positive, the entire sum is lower bounded by 1, hence
T log Σ 1 +
)
T; max, Iz
27
> t log 1) =0
Adding max, Tÿ to the inequality finishes the proof for the lower bound.
Regarding the upper bound, (21) can be maximized when x; max, æj for all i, in which case each (z: - max, 22,² is minimized, and hence each summand is maximized. In this case,
>T
27
T log Σ
1 +
2 max, 23
S T log
Σ1 = 7 log(d).
Adding max, 2ÿ to the inequality finishes the proof for the upper bound.
Fat Sigmoid Notably, we encountered similar problem using regular (log)-sigmoids smooth the constraint indicators for El with black-box constraints. In principle the Cauchy cummulative distribution function would satisfy these conditions, but requires the computation of arctan, a special function that requires more floating point operations to compute numerically than the following function. Here, we want the smooth approximation
satisfy 1) positivity, 2) monotonicity. 3) polynomial decay, and (xx) 1/2 - (--x). yet = /1/3, then we define
1 + (x - 7)))
x < O,
(xx)
1 1- 3 1+ (xr (cr 7)²)
x 2 0.
22

---

1.0
0.8
-20
-40
0.4
-60
0.2
-80
cig
0.0
skgma
-10
-5
10
-10 -5
0
10
Figure 13: We construct the fat sigmoid ap- Figure 14: The fat sigmoid approximation 0.01) decays as O(1/x²) instead of ((exp(=)) as + -00 minimizing the dynamic range of the numerical quantities in constrained qLogEI.
proximation (purple) by splicing together two Lorentzians (blue and teal) one their inflec- tion points.
c is monotonically increasing, satisfies u(z) 1 as oc, (0) 1/2. and dxx) O(1/x²) as æ + -oo. Further, we note that the asymptotics are primarily important here, but that we can also make the approximation tighter by introducing temperature parameter т, and letting 1- (=(x) = (x/r). The approximation error of i- (x() to the Heaviside step function becomes tighter point-wise as 0+, except for the origin where (x(x) 1/2, similar to the canonical sigmoid.
A.5
Constrained Expected Improvement
For the analytical case, many computational frameworks already provide numerically stable imple- mentation of the logarithm of the Gaussian cummulative distribution function, in the case of Py PyTorhh torch. srrch.speiia... llg_llognnneeel ndtr, which can be readily used in conjunction with our implementation of LogEI, as described Sec. [4.
For the case of Monte- Carlo parallel EI, we implemented the fat-tailed function from Sec. A.4 to approximate the constraint indicator and compute the per-candidate, per-sample acquisition utility using
(logsoftplus, (Sx(x)) v²) Σ log a (-5" K
5: (x;)
Teoms
where ((ah) is the ith sample of the kth constraint model, and Teoms the temperature parameter control- ling the approximation to the constraint indicator. While this functionality in our implementation, our benchmark results use the analytical version.
A.6 Parallel Expected Hypervolume Improvement
The hypervolume improvement can be computed via the inclusion-exclusion principle, see [3 for details, focus on the numerical issues concerning qEHVI here. To this end, we define
(m) R..i.....
min [uz f(x;;, ...,, f(xi,))
where f the vector-valued objective function, and uy is the vector of upper bounds of one of K ryper-rectangles that partition the wwre-Paaotiinaaaaee space, see [13) for details the partitioning. Letting 1 be the corresponding lower bounds of the hyper-rectangles, the hypervolume improvement then be computed
M (m)
HVI((
Σ
(-1)+
II
FE,X
(m0)
(m))
(23)
k=1j=1 XyEX3
where X, (X, C Xemi : |Xy] = j} the superset of all subsets of Xami of size and (m) (m) for X, {xi,...., }.
2R,X, 2
23

---

To find numerically stable formulation of the logarithm of this expression, we first re-purpose the Pmax function to compute the minimum in the expression of (m) " like so (P (4)
-P-max(-z) Further, we use the 4+ function of Sec. A.4 as for the single objective case to approximate R k 14 We then have
log II P
(m)) 8+ P-
Σ log + &+ 1m
(m)
(24)
Since we can only transform positive quantities to log space, we split the sum in Eq. 23) into positive and negative components, depending on the sign of (-1)++1 and compute the result using numerically stable implementation of log(cxp((1)) (cxp(log of positive terms) - exp(log of negative terms). The remaining sums over k and q can be carried out by applying logsumexp to the resulting quantity. Finally, applying logsumexp to reduce over an additional Monte-Carde sample dimension yields the formulation of qLogEHVI that we use in our multi-objective benchmarks.
A.7 Probability of Improvement
Numerical improvements for the probability of improvement acquisition which is defined as a(x) = [(x)=p" where the standard Normal CDF, can be obtained simply by taking the logarithm
((x)
using a numerically stable implementation of log( (()) = logerfo logerfe is computed as
+ - loog()), where
logerfc(x) = [1og(erf((z)) x x 0
1oog (eeffeff()))) x x> 0.
A.8 q-Noisy Expected Improvement
The same numerical improvements used by gLogEI to improve Monte-Carlo expected improvement (qEI) in Appendix can be applied improve the fully Monte Carlo Noisy Expected Improvement 6,.52 acquisition function. As in gLogEI, we can (i) approximate the max(0, x) using softplus to smooth the sample-level improvements to ensure that they are mathematically positive, (ii) approximate the maximum over the q candidate designs by norm || I and (iii) take the logarithm of the resulting smoothed value mitigate vanishing gradients. To further mitigate vanishing gradients, we can again leverage the Fat Softplus and Fat Maximum approximations. The only notable difference in the qEI and &NEI acquisition functions is the choice of incumbent, and similarly only a change of incumbent is required to obtain qLogNEI from qLogEI. Specifically, when the scalar y" in Equation (Tpp is replaced with vector containing the new incumbnent under each sample wer obtain the LLogNEI acquisition value. The hhh element of the incumbent vector nty for NEogNEI is max max E' (x where Xq+1. Xx+t are the previously evaluated designs and (((xx) ) is the value of the 7thh point under the :th sample from the joint posterior over XI, X We note that we use a hard maximum to compute the incumbent for each sample because we do not need compute gradients with respect to the previously evaluated designs X.+1, .... Xn+
We note that we obtain computational speed ups by (i) pruning the set of previously evaluated points that considered for being the best incumbent to include only those designs with non-zero probability of being the best design and (ii) caching the Cholesky decomposition of the posterior covariance over the resulting pruned set of previously evaluated designs and using low-rank updates for efficient sampling [4].
For experimental results of gLogNEI see Section D.4
B Strategies for Optimizing Acquisition Functions
As discussed in Section 2.3 a variety of different approaches and heuristics have been applied to the problem of optimizing acquisition functions. For the purpose of this work, we only consider continuous domains X. While discrete and/or mixed domains are also relevant in practice and have
24

---

received substantial attention in recent years see e.g. Baptista and Poloczek 77]. Daulton et al. [5]. Deshwal et al. [8], Kim al. 44, Oh et al. [62]. Wan et al. 74 our work here on improving acquisition functions is largely orthogonal to this (though the largest gains should be expected when using gradient-based optimizers, as is done in mixed-variable BO when conditioning on discrete variables, or when performing discrete or mixed BO using continuous relaxations, probabilistic reparameterization, or straight-through estimators (15).
Arguably the simplest approach to optimizing acquisition functions is by grid search or random search. While variants of this combined with local descent can make sense in the context of optimizing over discrete or mixed spaces and when acquisition functions can be evaluated efficiently in batch (e.g. on GPUs), this clearly does not scale to higher-dimensione continuous domains due to the exponential growth of space to cover.
Another relatively straightforward approach is zeroth-order methods such DIRECT 41 (used e.g. by Dragonfly (43) or the popular CMA-ES [2). These approaches are easy implement as they avoid the need compute gradients of acquisition functions. However, not relying on gradients is also what renders their optimization performance inferior to gradient based methods, especially for higher-dimensional problems and/or joint batch optimization in parallel Bayesian optimization.
The most common approach to optimizing acquisition functions on continuous domains is using gradient descent-type algorithms. Gradients are either be computed based on analytically derived closed-form expressions, or via suto-differentiation capabilities of modern ML systems such as PyTorch 633 Tensorflow [1, or JAX [9].
For analytic acquisition functions, common choice of optimizer L-BFGS-B 10] a quasi-second order method that uses gradient information to approximate the Hessian and supports box constraints. If other, more general constraints are imposed on the domain, other general purpose nonlinear optimizers such SLSQP 47 or IPOPT 773 are used (e.g. by BoTocch). For Monte Carlo (MC) acquisition functions, Wilson et al. (77 proposes using stochastic gradient ascent (SGA) based on stochastic gradient estimates obtained via the reparameterization trick (45). Stochastic first-order
algorithms are also used by others, including e.g. Wang et al. [75] and Daulton et al. [15]. Balandat et al. [6] build on the work by Wilson et al. [77] and show how sample average approximation (SAA) can be employed to obtain deterministic gradient estimates for MC acquisition functions, which has the advantage of being able to leverage the improved convergence rates of optimization algorithms designed for deterministic functions such as L-BFGS-B. B. This general approach has since been used for variety of other acquisition functions, including e.g. Daulton al. 13 and Jiang al. 40
Very few implementations of Bayesian Optimization actually use higher-order derivative information, as this either requires complex derivations of analytical expressions and their custom implementation, or computation of second-order derivatives via automated differentiation, which is less well supported and computationally much more costly than computing only first-order derivatives. One notable exception is Cornell-MOS 788 79), which supports Newton's method (though this is limited to the acquisition functions implemented in C++ within the library and not easily extensible to other acquisition functions).
B.1 Common initialization heuristics for multi-start gradient-descent
One of the key issues to deal with gradient- optimization in the context of optimizing acquisition functions is the optimizer getting stuck in local optima due to the generally highly non-convex objective. This is typically addressed by means of restarting the optimizer from a number of different initial conditions distributed across the domain.
A variety of different heuristies have been proposed for this. The most basic one is to restart from random points uniformly sampled from the domain (for instance, scikit- tootiiiitttttteeeeeee [33] uses this strategy). However, as we have argued in this paper, acquisition functions can be (numerically) zero in large parts of the domain, and so purely random restarts can become ineffective, especially in higher dimensions and with more data points. A common strategy is therefore to either augment or bias the restart point selection to include initial conditions that are closer to "promising points". CPyOpt [69] augments random restarts with the best points observed so far, or alternatively points generated via Thompson sampling. Spearmint (66) initializes starting points based on Gaussian perturbations of the current best point. BoTorch (6) selects initial points by performing Boltzmann sampling on a set
of random points according to their acquisition function value; the goal of this strategy is to achieve
25

---

biased random sampling across the domain that is likely to generate more points around regions with high acquisition value, but remains asymptotically space-filling. The initialization strategy used by Trieste (64] works similarly to the one in BoTorch, but instead of using soft-randomization via Boltzmann sampling, simply selects the top-k points. Most recently, Gramacy et al. [31] proposed distributing initial conditions using Delaunay triangulation of previously observed data points. This is an interesting approach that generalizes the idea of initializing "in between" observed points from the single-dinensional case. However, this approach does not scale well with the problem dimension and the number of observed data points due to the complexity of computing the triangulation (with wall time empirically found to be exponential in the dimension, see 31 Fig. 3] and worst-case quadratic in the number of observed points).
However, while these initialization strategies can help substantially with better optimizing acquisition functions, they ultimately cannot resolve foundational issues with acquisition functions the. Ensuring that acquisition functions provides enough gradient information (not just mathematically but also numerically) is therefore key to be able to optimize it effectively, especially in higher dimensions and with more observed data points.
C Proofs
Theorem 1. Suppose f is drawn from a Gaussian process prior Pj. 1² S Pns On are the mean and standard deviation of the posterior (f()) and B E R. Then with probability 1 - 6,
Mn (x)
V¹
Px
B Px (f(x) < 5 cn)
on (x)
where (f- + (V 2og 28) B) maxx on (x).
Proof. We begin by expanding the argument to the h function in Eq. (2) as a sum of 1) the standardized error of the posterior mean (4n to the true objective ʃ and 2) the standardized difference of the value of the true objective f at x to the best previously observed value y" max, yr:
(x(x)) - """ Pn (x) s(x) + [(x) xx
(25)
(x)
(x)
We proceed by bounding the first term on the right hand side. Note that by assumption, f(x) (fnn( (x), o. )x)²) and thus (ux (x) - f(x))
(x) ~ ~N(,,) 1). For a positive C >> 0 then, we use
a standard bound on the Gaussian tail probability to attain
Pn [n(x)) f(x)
P
C)
S se 121222
(26)
(x)
Therefore, (p(x) - f(x))/oo < C with probability - 8 if C= - log(25).
Using the bound just derived, and forcing the resulting upper bound to be less than B yields a sufficient condition to imply (x) - V² < (x):
(x - 3"' S CC+ f(x) - 3² < B
(27)
on (x)
(x(x)
Letting fz ;=
;= f(xx), re-arranging and using 32 = fO + (y- ff") we get with probability 1 8,
f(x) < 8 - ( - (v/ -2 =olgg(2)) - (x)x
(28)
Thus, we get
Mn (x)
- 3"
Px
< 2 Px [((x) $ -(f -z") - (V-22og(26)) 3)))(z))
(x)
2 Pa [(x) (f' - v¹) 22V22log(22 B) max on (x)
(29)
Note that the last inequality gives a bound that is not directly dependent on the evaluation of the posterior statistics of the surrogate at any specific x. Rather, it is dependent on the optimality
gap f y² and the maximal posterior standard deviation, or a bound thereof. Letting en (f* - -")) - (o-2log(28) - B) maxx on (x) finishes the proof.
26

---

Lemma 2. [Relative Approximation Guarantee] Given To, Tmax > 0, the approximation evvor of gLogEl qEI is bounded by
|e (L QEI(X) S (q"max 1) qEI(X) + log(2)roq"mmx Tmaxx
(11)
Proof. Let Big (x,,) yr where i {L,..., n), and for brevity of notation, and let 1se, 1sp refer to the logsumexp and logsoftplus functions, respectively, and (eLU(x) = max(2; 0). We then bound n/ - qEI(X) by
exp( (se (2ig//Tmax)) Fmax))) max (21g)
exp(T Cop((maxlse,, (1sp,, (mmx))) max ReLU zin)
Σ
[|softplus, (2). 111/Tmxx max (eLU(
(30)
$ [softplus, (2). 111/Tmox max softplus,, (Ziq)
max softplus, (Zig) max ReLU =iig)
First and second inequalities are due to the triangle inequality, where for the second we used la cl S a b| + cl with b = max, softplus(zg))
To bound the first term in the sum, note that ||x| S ||x||g I ||x thus 0 S (/xxllg ||x| S d=/4 x||x||a 1x.. and therefore
|softplus,, (21))11/Tmx max moxtpllplus (Zig)
S (g - 1) max maxxoftplus, (Zig)
$ (q mmxx 1) ReaxxeeLU(s (21g) + 12g(2)
The second term in the sum be bound due |softplus (x) ReLU(x)] < log(() (see
Lemma below) and therefore.
max softplus, (Zig) max (zin))) S 1og(2)7o-
Dividing Eq. (30) by n to compute the sample mean finishes the proof for the Monte-Carlo approx- imations to the acquisition value. Taking n + oo further proves the result for the mathematical definitions of the parallel acquisition values, i.c. Eq. (4)
Approximating the ReLU using the softplus, (x) T logg(1 + exp(x/r)) function leads to an approximation етrог that is at most T in the infinity norm, i.e. ||ootpluu,, ReLU|| = (2) The following lemma formally proves this.
Lemma 7. Given
0. we have for all R.
ssotplus, (x) - ReLU(x)) & log(2) т.
(31)
Proof.
Taking the Subb)dderivative of softplus, ReLU, we get
0,softplus, (x) ReLU(x) = (1 + - =/r)- 1 x>0 0
0 x
which is positive for x and negative for all x > 0. hence the extremum must be x, at which point softplus, (0) ReLU(0) = log(2)r. Analyzing the asymptotic behavior,
(smp==us (x) (z) 0, and therefore softplus, (x) ReLU(x) for R. O
Approximation guarantees for the fat-tailed non-linearities of App. A.4 can be derived similarly.
27

---

D Additional Empirical Details and Results
D.1 Experimental details
All algorithms are implemented in BoTorch. The analytic EI, qEI, cEI utilize the standard BoTorch implementations. We utilize the original authors' implementations of single objective JES (36). GIBBON [61], and multi-objective JES [7]], which are all available in the main BoTorch repository. All simulations are ran with 32 replicates and error bars represent 12 times the standard error of the mean. We use a Matern-5/2 kernel with automatic relevance determination (ARD), i.e. separate length-scales for each input dimension, and a top-hat prior on the length- magtt--ccales in [0.01, 100). The input spaces are normalized to the unit hyper-cube and the objective values are standardized during each optimization iteration.
D.2
Additional Empirical Results on Vanishing Values and Gradients
The left plot of Figure D in the main text shows that for a large fraction of points across the domain the gradients of EI are numerically essentially zero. In this section we provide additional detail on these simulations as well as intuition for results.
The data generating process (DGP) for the training data used for the left plot of Figure T is the following: 80% of training points are sampled uniformly at random from the domain, while 20% sampled according a multivariate Gaussian centered at the function maximum with a standard deviation of 25% of the length of the domain. The idea behind this DGP is to mimic the kind of data would see during Bayesian Optimization loop (without having to thousands of BO loops to generate the Figure I). Under this DGP with the chosen test problem, the incumbent (best observed point) is typically better than the values the random test locations, and this becomes increasingly the case as the dimensionality of the problem increases and the number of training points grows. This is exactly the situation that is typical when conducting Bayesian Optimization.
orann d-H8, 2-10
Histogram
800
devovili too 1
Bamavki NDI < le
700
Taeavalii kot lell
-17
600
-19
600
-22
100
URIRA smmple innnppe
-33
100
-10
-25
Figure 15: Model fits for typical replicate used in generating the Fig. [let)) While there is am- ple uncertainty in the test point predictions (blue, chosen uniformly at random), the mean predic- tion for the majority of points is many standard deviations away from the incumbent value (red).
Figure 16: Histogram of (x)) the argument to h in 2 corresponding to Figure [5] Vertical lines are the thresholds corresponding to values = below which h((z) is less than the respective threshold. The majority of the test points fall below these threshold values.
For a particular replicate, Figure 15 shows the model fits in-sample (black), out-of-sample (blue), and the best point identified so far (red) with 60 training points and a random subset of 50 (out of 2000) test points. One can see that the model produces decent mean predictions for out-of-sample data, and that the uncertainty estimates appear reasonably wall-calibrated (e.g., the credible intervals typically cover the true value). A practitioner would consider this a good model for the purposes of Bayesian Optimization. However, while there ample uncertainty the predictions of the model away from the training points, for the vast majority of points, the mean prediction is many standard deviations away from the incumbent value (the error bars are standard deviations). This is the key reason for El taking zero (or vanishingly small) values and having vanishing gradients.
28

---

To illustrate this, Figure 16 shows the histogram of 2(x:) values. the argument to the function h in 2 It also contains the thresholds corresponding to the values below which h(2) is less than the respective threshold. Since a (ar) is close to for most test points (mean: 0.87, std: 0.07), this is more or less the same as saying that (x)) is less than the threshold. It is evident from the histogram that the majority of the test points fall below these threshold values (especially for larger thresholds),
showing that the associated acquisition function values (and similarly the gradients) are numerically almost zero and causing issues during acquisition function optimization.
D.3 Parallel Expected Improvement
9.
150
150
130
15.0
Aslo 440
Boi Meriidd vale
12.5
100
100
164
7.5
7.5
50-
5.
20
250
20
100
250
90
80
43
30-
32
Panciin
Pandkan
Plackkn cokatton
Plackan evalation
Figure 17: Parallel optimization performance the Ackley and Levy functions in 16 dimensions. GLogEl outperforms all baselines on Ackley, where joint optimization of the batch also improves on the sequential greedy. On Levy, joint optimization of the batch with LlogEl starts out performing worse in terms of BÔ performance than sequential, but wins out over all sequential methods as the batch size increases.
Suat www.naadaa aad
14
34
12
1-
300
166
1.00
10:
80
80
160
Pannkos
Figure 18: Breakdown of the parallel optimization performance of Figure 17 per method, rather than per batch size. On Levy, qLogEI exhibits a much smaller deterioration in BO performance due to increases parallelism than the methods relying on sequentially optimized batches of candidates.
Figure 17 reports optimization performance of parallel BO on the 16-dimensional Ackley and Levy functions for both sequential greedy and joint batch optimization. Besides the apparent substantial advantages of qLogEI over qEI on Ackley, key observation here is that jointly optimizing the candidates of batch acquisition functions can yield highly competitive optimization performance, especially as the batch size increases. Notably, joint optimization of the batch with gLogEl starts out performing worse in terms of BO performance than sequential on the Levy function, but outperforms all sequential methods as the batch size increases. See also Figure T8 for the scaling of each method with respect to the batch size q.
29

---

Noise: 1.0% Range(f)
Noise: 2.0% Range(f)
Noise: 5.0% Rangel()
Best observed value
rand
Hartmann 6D
qnei
-2-
qei
glognei-
glogei
-33
gibbon
100
200
200
100
200
16-
Best observed value
14
14
Ackley 8D
12.
12
12
10
10
8
8
6
4
2
2
100
200
100
200
16
Best observed value
14-
Ackley 16D
12-
10.
8-
8
6
6
4
2
2
100
200
100 200
Function evaluations
Function evaluations
Function evaluations
Figure 19: Optimization performance with noisy observations on Hartmann 6d (top), Ackley 8d (mid), and Ackley 16 (bottom) for varying noise levels and 9 1. We set the noise level proportion of the maximum range of the function, which R 3.3 for Hartmann and 20 for Ackley. That is, the 1% noise level corresponds standard deviation of 0.2 for Ackley. qLogNEI outperforms both canonical El counterparts and Gibbon significantly in most cases, especially in higher dimensions.
D.4 Noisy Expected Improvement
Figure 19 benchmarks the "noisy" variant, gLogNEI Similar to the noiseless case, the advantage of the LogEl versions over the canonical counterparts grows as with the dimensionality the problem, and the noisy version improves the canonical versions for larger noise levels.
D.5 Multi-Objective optimization with gLogEHVI
commpress qLogEHVI and qEHVI on different test problems with objectives, and ranging from 2-30 dimensions. This includes real world inspired problems: cell network design for optimizing coverage and capacity [19], laser plasma accelcration optimization [38]. and vehicle design optimization [54] 68). The results are consistent with our findings in the singlc-objective and constrained cases: gLogEHVI consistently outperforms qEHVI, and the gap is larger on higher dimensional problems.
Figure compares
D.6 Combining LogEl with TuRBO for Highh-Dimensionll Bayesian Optimization
In the main text, we show how LogEl performs particularly well relative to other basclines in high dimensional spaces. Here, show how LogEl can work synergistically with trust region based methods for high-dimensional BO, such as TuRBO 24]
Fig. 21 compares the performance of LogEI, TuRBO-1 + LogEl, TuRBO-1 + EI, well as the based implementation for the 50d Ackley test problem. Combining TuRBO-1 with LogEI results in substantially better performance than the baselines when using small number of function evaluations. Thompson sampling (TS) ultimately performs better after 10, 000 evaluations, but this experiment shows the promise of combining TuRBO and LogEl in settings where we cannot do thousands of function evaluations. Since we optimize batches of g = 50 candidates jointly, also increase the number of Monte-Caie Carlo samples from the Gaussian process
original Thompson-sampe
30

---

Laser Plasma 4d=4,M=3)
Cell Network (d=30,
DTLZ2 (d=6,M=
0.06
0.4
Hypervolume
2-
0.04
0.2
200
0
100
50
100
ZITT (d=6,M
BraninCurrin (d=
Vehicle Safety (d=5, M =3)
250
120
Hypervolume
50
110
0
50 Function Evaluations
100
0
50
100
50
100
Function Evaluations
Function Evaluations
Rand
qEHVI
gLogEHVI
JES
Figure 20: Sequential (g = 1) optimization performance on multi-objectivee problems, as measured by the hypervolume of the Pareto frontier across observed points. This plot includes JES (71) Similar to the single-objective case, qLogEHVI significantly outperforms all baselines on all test problems.
50D Ackley, q=50
16
Regret
2
1
0.0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
Number evaluations
TURBO-1 1 dEL1 nottrt
TuR8o-1 + glnggE), msturt
el olog matort
-... qEl, 1 restart
TuR8O-1 (TR)
TURBO-1 GEL
TURBO-1 gLogEi, restatss
eLagEl,
restarts
Figure 21: Combining LogEl with TuRBO the high-dimensional the 50d Ackley problem yields significant improvement in objective value for small number of function evaluations. Unlike for qEI, no random restarts are necessary achieve good performance when performing joint optimization of the batch with qLogEl (q = 50).
from 128, the BoTorch default, 512, and use the fat-tailed smooth approximations of Sec. A.4 ensure strong gradient signal to all candidates of the batch.
Regarding the ultimate out-performance of TS over qLogEl, think this is due model mis- specification since the smoothness of a GP with the Matern-5/2 kernel cannot express the non- differentiability of Ackley at the optimum.
D.7 Constrained Problems
While running the benchmarks using CEI in section we found that we in fact improved upon best known result from the literature. We compare with the results Coello and Montes [T], which are generated using 30 runs of 80,000 function evaluations each.
31

---

For the pressure vessel design problem, Coello and Montes [1 quote a best-case feasible objective of 6059.99 9959.1143441. Out of just 16 different runs, LogEI achieves a worst-case feasible objective of 5659. 1 1008 after only 110 evaluations, and best case of 5651.8862, a notable reduction in objective value using almost three orders of magnitude fewer function evaluations.
For the welded beam problem, Coello and Montes [[]] quote 1.728226, whereas LogEI found best case of 1.7496 after 110 evaluations, which is lightly worse, but we stress that this is using three orders of magnitude fewer evaluations.
For the -esssiooncoeeesiionn problem, LogEI found feasible solution with value 0.0129 after 110 evaluations compared to the 0.012681 reported in in ..]
We emphasize that genetic algorithms and BO are generally concerned with distinct problem classes: BO focuses heavily on sample efficiency and the small-data regime, while genetic algorithms often utilize a substantially larger number of function evaluations. The results here show that in this case BO is competitive with and can even outperform genetic algorithm, using only tiny fraction of the sample budget, see App. D.7 details. Sample efficiency is particularly relevant for physical simulators whose evaluation takes significant computational effort, often rendering several tens of thousands of evaluations infeasible.
D.8 Parallel Bayesian Optimization with cross-batch constraints
In some parallel Bayesian optimization settings, batch optimization is subject to non-trivial constraints across the batch elements. A natural example for this are budget constraints. For instance, in the context of experimental material science, consider the case where each manufactured compound requires a certain amount of different materials (as described by its parameters), but there is only fixed total amount of material available (e.g., because the stock is limited due to cost and/or storage capacity). In such a situation, batch generation will be subject to a budget constraint that not separable across the elements of the batch. Importantly, in that case sequential greedy batch generation is not an option since it is not able to incorporate the budget constraint. Therefore, joint batch optimization is required.
Here we give one such example in the context of Bayesian Optimization for sequential experimental design. We consider the -ie-dimeniionll silver nanoparticle flow synthesis problem from Liang et al. [53]. In this problem, to goal is to optimize the absorbance spectrum score of the synthesized nanoparticles over five parameters: four flow rate ratios of different components (silver, silver nitrate, trisodium citrate, polyvinyl alcohol) and a total flow rate Qtot-
The original problem was optimized over a discrete set of parameterizations. For our purposes we created a continuous surrogate model based on the experimental dataset (available from https by fitting an RBF interpolator (smoothing factor of 0.01) in scipy on the (negative) loss. We use the same search space as Liang et al. [53], but addition to the box bounds on the parameters we also impose an additional constraint on the total flow rate Que = 2000 /L/mmin across the batch: Σ Qiot S Q (the maximum flow rate per syringe / batch element is 1000nL/mmnn. This constraint expresses the maximum throughput limit of the microfluidic experimentation setup. The result of this constraint is that we cannot consider the batch elements (in this case automated syringe pumps) have all elements of a batch of experiments operate in the figh-flow regime at the same time.
7/github. wwtp.//--aa//eechhheeeeo
In our experiment, we use batch size of 3 and start the optimization from 5 randomly sampled points from the domain. We run 75 replicates with random initial conditions (shared across the different methods), error bars show + two times the standard error of the mean. Our baseline is uniform random sampling from the domain (we use a hit-and-run sampler to sample uniformly from the constraint polytope Σ Qiot S QmaF). We compare qEI vs. qLogEl, and for each of the two we evaluate (i) the version with the batch constraint imposed explicitly in the optimizer (the optimization in this case uses scipy's 's SLSQP solver), and (ii) heuristic that first samples the total flow rates {Qi.r 1i_ uniformly from the constraint set, and then optimizes the acquisition function with the flow rates fixed to the sampled values.
The results in Figure show that while both the heuristic ("random Qtat") ") and the proper constrained optimization ("batch-constrained Qu.at") substantially outperform the purely random baseline, it requires uisng both LogEl and proper constraints to achieve additional performance gains over the other 3 combinations. Importantly, this approach is only possible by performing joint optimization of
32

---

Best observed point
Best inferred point
-0.3
-0.4
negative loss
-0.5
-0.6
Ranlom
Random
gEL, batch-constrineed Qae
betch-constrained Qui
-077
LogEE batch-constrained Qui
gLogEl, batch-constrained Q
gEL random Qe:
rundom Qai
-0.8
Log radom Q..
gLogEl, random Que
Number evaluations
Number of batches
Figure 22: Optimization results on the nanomaterial synthesis material science problem with cross- batch constraints. While qLogEI outperforms qEI under the proper constrained ("batch-cone ceatch-consttained Qiot") optimization, this is the case for the the heuristic ("random Qiot") ",o"",, demonstrating the value of both joint batch optimization with constraints and LogEl.
the batch, which underlines the importance of qLogEl and its siblings being able achieve superior joint batch optimization in settings like this.
D.9 Details on Multi-Objective Problems
We consider variety of multi-objective benchmark problems. We evaluate performance three synthetic biobjective problems Branin-Currin (d = 2) [8, ZDTI (d = 6) [83] and DTLZ2 (d = 6) three real world inspired problems. For
[17. As described in we also evaluated performance the laser plasma acceleration problem, we used the public data available at Irshad et al. [39] to fit an independent GP surrogate model each objective. We only queried te surrogate at the highest fidelity to create a single fidelity benchmark.
D.10 Effect of Temperature Parameter
In Figure 23
examine the effect of fixed r for the softplus operator optimization performance. find that smaller values typically work better.
D.11 Effect the initialization strategy
Packages and frameworks commonly utilize smart initialization heuristics to improve acquisition function optimization performance. In Figure 24 we compare simple random restart optimization, where initial points are selected uniformly at random, with Bo BoTTorch's default initialization strategy, which evaluates the acquisition function on large number of points selected from a scrambled Sobol sequence, and selects points at random via Boltzman sampling (c.g., sampling using probabilities computed by taking a softmax over the acquisition values [6. Here we consider 1024 initial candidates. We find that the BoTorch initialization strategy improves regret for all cases, and that gLogEI, followed by UCB show less sensitivity to the choice of initializations strategy. Figure25 examines the sensitivity of qEI the number of initial starting points when performing standard random restart optimization and jointly optimizing the points in the batch. We find that, consistent with empirical and theoretical results in the main text, qEI often gets stuck in local minima for the Ackley test function, and additional random restarts often improve results but do not compensate for the fundamental optimality gap. The performance of gLogEl also improves the number of starting points increases.
33

---

d-2
tog 2000
gett 00000
101
tofute lets
10-1
6X10
4x10
10-2
5x.10
10-
Figure 23: Ablation study on the convergence characteristics of LogEI on Ackley and sum of squares (SOS) problems in 2 and 16 dimensions. The study shows that it is important to choose small Tu for the best convergence properties, which results in very tight approximation to the original ReLU non-linearity in the integrand. Critically, setting to as low as 10 is only possible due to the transformation of all computations into log-space. Otherwise, the smoothed acquísition utility would exhibit similarly numerically vanishing gradients as the original ReLU non-linearity.
CELL NETWORK
BRANIN-CURRIN
DTLZ2
LASER PLASMA
ZDTI
VEHICLE SAFETY
JES
21/6 1.1)
(9/-
33/6 1.0)
$7.3 (+/- 0.7)
72.7 1.0)
47.0 (+/- 16)
QEHVI
(4/- 0.0)
(+/- 0.0)
1.0 (+)- 0.0)
3.0 (+/- 0.1)
(+/- 0.0)
(+/- 0.0)
OLOGEHVI
9.2 (4/- 0.8)
10.0 (+/-0.4) 0.4)
5.8 (+/- 0.21
31.6 (+/ 1.7)
7.2 (4/- 0.7)
2+/-0.11
RAND
0.2 (+/- 0.0)
0.2 (+/- 0.01
0.2 (+)- 0.01
0./ 0.0)
0.3 (+7- 0.0)
0.3 (+/ 0.0)
Table 1: Multi-objectve acquisition function optimization wall time in seconds on CPU (2x Intel Xeon E5-2680 2.40GHz). We report the mean and + standard errors.
34

---

2D
8D
16D
16
Best abserved value
10
14
14
12
Ackley
10
10
5
8
2
100
200
100
200
100
200
-2
Best observed value
Michalewiez
-1.0
-1.5
100
200
Function Evaluations
Function Evaluations
Function Evaluations
EI Random Init.
Bl Bolizmann Init.
LogEl Random Init.
LogEl Boltzmann Init.
GIBBON Random Init.
GIBBON Boltzmann Inil.
Random Init.
Boltzmann Init.
UCB Random Init.
UCB Boltzmann Init.
Figure 24: Sensitivity to the initialization strategy. Random selects random restart points from the design space uniformly at random, whereas Boltzmann initialization the default BoTorch initialization strategy which selects points with higher acquisition function values with higher probability via Boltzmann sampling.
35

---

Jery
Steat: Stattoatts
Steptcram
-00
230
200
beratioo
ookley
Stade
-0-
-00
12
-10
wokley
points
-10
eturrig polces
461-
Kaninng park
staring poin
481-18
stattng pount
Sundic Sundicrge
-40
-70
Figure 25: Sensitivity number of starting points with multi-start optimization for the 16D Ackley and Levy test problems. Note: We plot negative regret, so higher is better.
36


Reviews:
[{"review_id": "t7Gv3kNBxp", "reviewer": "Reviewer_QK5L", "summary": "In this paper, the authors identify both through examples and theoretical analysis several numerical pathologies inherent to the computation and optimization of the Expected Improvement (EI), a popular acquisition function at the heart of Bayesian Optimization (BO) algorithms. \nThey subsequently propose a numerical reformulation, LogEI, which achieves substantially better performances on a quite extensive range of benchmarks. This reformulation applies to all the member of the EI family: constrained EI for constrained BO, parallel EI for batch BO, and expected hypervolume improvement for multi-objective BO.", "questions": "I do not have any questions.\n\nEdit: I have read the rebuttal and the discussions between the authors and other reviewers, my score remains unchanged.", "limitations": "None.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written and well-organised.\n- The proposed numerical fix for EI is likely to have a great impact as it may benefit to all public implementations of EI. Furthermore, it does not incur excessively longer computation times except perhaps for multi-objective BO (roughly one order of magnitude larger).\n- LogEI seems to produce more consistent results with respect to the initial optimization starting point compared to its canonical counterpart EI, thus reducing the need for heuristics for that matter.\n- All claims are backed up by an impressive amount of numerical experiments and ablation studies in each setting (vanilla BO/ constrained BO/ batch BO/multi-objective BO).", "weaknesses": "I did not spot any weakness.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your positive review."}, {"review_id": "M5SZNhJs4n", "reviewer": "Reviewer_o9Cr", "summary": "This paper identifies a numerical pathology with the expected improvement (EI) family of acquisition functions: the vanishing gradients of the acquisition function leads to failure in acquisition function optimization. A set of modified EI acquisition functions that fix the numerical pathologies have been proposed. In experiments, the proposed EI acquisition functions are shown to outperform the canonical EI acquisition functions and to perform on par with other state-of-the-art acquisition functions on Bayesian optimization benchmarks.\n", "questions": "- It is great to have an error bound of the qLogEI. However, for acquisition functions, preserving the relative order of values is more important than absolute difference. I wonder to what extent qLogEI preserves the relative order of values compared to qEI.\n\n", "limitations": "The limitation of the proposed method has been discussed in the last section.\n", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- This paper focuses on a previous neglected aspect of Bayesian optimization, acquisition function optimization, and identifies the numerical pathologies associated with EI acquisition functions. This paper provides a theoretical analysis of the vanishing gradient issue.\n- Improvements have been proposed to EI, Monto Carlo Parallel EI, Constrained EI and EHVI. \n- Experiments clearly show the numerical pathology of EI optimization and the superior performance of the improved EI version.\n- The proposed improved EI acquisition functions perform on par with state-of-the-art acquisition functions on high dimensional synthetic functions in both sequential and batch settings.\n", "weaknesses": "- The proposed treatment only works for the acquisition functions with vanishing gradient issues.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your encouraging your review. \n\n> It is great to have an error bound of the qLogEI. However, for acquisition functions, preserving the relative order of values is more important than absolute difference. I wonder to what extent qLogEI preserves the relative order of values compared to qEI.\n\nWhile the smooth approximation to the ReLU in the integrand of qEI is monotonically increasing and upper-bounds the ReLU, the smooth approximation to the max operator can change the relative ordering. Therefore, while the relative ordering is fully maintained for the analytical LogEI version, this cannot generally be guaranteed for the Monte-Carlo (batch) version. \n\nNevertheless, the absolute error bound guarantees that the maximizer of qLogEI attains a similar acquisition value to the maximizer of qEI. Whether or not that is a similar point in the input space depends on the surrogate model during the particular iteration. If the EI acquisition value is a good indicator of optimization performance, such a change in ordering would not lead to a significant decrease in optimization performance, and if EI is not a good indicator of optimization performance on a particular problem, it would not be recommended to use EI in any case. \n\nNotably, our empirical experiments demonstrate that the benefits of a smooth acquisition optimization landscape with strong gradients far outweigh the potential for a small change in relative ordering of points with promising acquisition values."}, {"review_id": "ti5YE6UJBH", "reviewer": "Reviewer_u8yk", "summary": "The paper proposes LogEI, family of acquisition functions with improved numerical stability over EI that makes it more suitable for gradient-based acquisition function optimization, all while retaining similar optima as EI. Pathologies of EI are visualized and analyzed, and the approximation error between qEI and qLogEI is theoretically bounded. Empirically, LogEI clearly outperforms EI on most tasks, suggesting that it can act as a drop-in replacement for EI.", "questions": "- Is there any setting where LogEI can _not_ act as a plug-in replacement for EI, or where performance would be expected to be worse?\n- Have the authors experimented with LogEI non-continous search spaces, and if so, what are the findings?", "limitations": "Some suggestions for addressing limitations have been stated in the Weaknesses section, but are otherwise adequately addressed.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "__Good motivation__: Acquisition function optimization is an often-overlooked aspect of Bayesian optimization, and the paper does a good job of displaying the difficulties of acquisition function optimization (Fig. 1) and how the proposed approach remedies the issue.\n__Simple, effective and extensible solution:__ Simple solutions that work are great, and LogEI (and its extensions) is good example. \n__Very good empirical performance:__ The improvement over EI in the results is striking on most tasks, suggesting it is simply a superior acquisition function to the default.", "weaknesses": "__Anecdotal evidence for similarity with EI:__ Intuitively, It is sensible that LogEI is similar to EI. However, there is little evidence (Fig. 1) and theory (Lemma 2) to support this.  I would greatly appreciate a similar Lemma for the analytic variant, and examples of when the two may not be identical. The performance gain of LogEI compared to EI is rather substantial on Ackley (~4x on Ackley-16!) and Mich, which suggests that the two may in fact not be very similar (but that LogEI may in fact simply be superior). Specifically, I don't believe the statement in Row 9, \"LogEI, whose members either have _identical or approximately equal optima_ as their canonical counterparts\" is well supported.\n\n__Existing LogEI and Lacking references to Related work:__ The idea of a LogEI is not novel (LogEI was an acquisition function in SMAC at one point). Admittedly, that implementation regards log-transformations of the objective [1, 2] and would not help with numerical stability in the same manner. Nevertheless, I think these warrant citation and  _comparison in the experiments_, given the similarities. Moreover, it limits the novelty of the approach.\n\n__Relevance in high dimensions:__ Currently, I am not convinced by the justification for LogEI in high dimensions. To me, LogEI aims to address the pathology where the acquisition function is zero, but I don't see that happening in high-dimensional problems due to the exuberance of high-uncertainty regions (which would make the acquisition function _non-zero, but constant?_). So, why is the proposed approach particularily important in high dimensions - i.e., why does LogEI help when the acquisition function is constant as opposed to (almost) zero? This would, in my opinion, require a separate motivation than Fig. 1, empirical results aside. Moreover, I think the _zero-value_ versus _constant-value_ distinction is very important, and should be emphasized more.\n\n  With this in mind, I find Fig. 2 striking and odd. 60 data points (which is when almost all points have a zero-valued gradient) on an 8-dimensional problem is not a large amount of data (not even a 2x2x...x2 grid), yet the uncertainty is small to the point of \"EI and its gradients become numerically zero across most of the domain\"? With all due respect, are the authors sure that this is not _just_ the gradient (and not the function value), or that the model has oddly long lengthscales?\n\n__Minor:__\n  - __Noisy tasks:__ Adaptation of LogEI to noisy problem settings are missing\n  - __Lack of conventional benchmarks__: As a potential drop-in replacement for EI, seeing its performance of the method on the most conventional low-dimensional tasks (Branin, Hartmanns) would be informative. Moreover, it would be helpful for future benchmarking.\n\nFor all of these bullets, I believe that evidence to clarify (and not necessarily disprove) the remarks would substantially strengthen the paper. \n\n[1]. An experimental investigation of model-based parameter optimisation: SPO and beyond. F Hutter, HH Hoos, K Leyton-Brown, Kevin P. Murphy. _GECCO '09: Proceedings of the 11th Annual conference on Genetic and evolutionary computation_. 2009.\n\n[2] Sequential model-based optimization for general algorithm configuration. F Hutter, HH Hoos, K Leyton-Brown. _Learning and Intelligent Optimization: 5th International Conference_. 2011.\n\nLogEI in SMAC: https://github.com/automl/SMAC3/blob/29355618b35dcf4b3ce3e773d633109f036dba17/smac/optimizer/acquisition.py#L503", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed feedback about areas that deserve additional discussion.  We seek to clarify the points raised in the following. \n\n__Equivalence of optima of analytical LogEI and EI__\n\nIn the CR, we will clarify this statement through a brief Lemma. If the maximum of EI is greater than 0, LogEI and EI have the same set of maximizers. Furthermore, if $\\max_{x \\in \\mathbb X} EI(x) = 0$, then $\\mathbb X = \\arg \\max_{x \\in \\mathbb X} EI(x)$. In this case, LogEI is undefined everywhere, so it has no maximizers, which we note would yield the same BO policy as EI (where every point is a maximizer). The Lemma is as follows:\n\n_Lemma_: If $\\max_{x \\in \\mathbb X} EI(x) > 0$, then $\\arg \\max_{x \\in \\mathbb X} EI(x) = \\arg \\max_{x \\in \\mathbb X, EI(x) > 0} LogEI(x)$.\n\n_Proof_:\nSuppose $\\max_{x \\in \\mathbb X} EI(x) > 0$. Then $\\arg \\max_{x \\in \\mathbb X} EI(x) = \\arg \\max_{x \\in \\mathbb X, EI(x) > 0} EI(x)$.\nFor all $x \\in \\mathbb X$ such that $EI(x) > 0$,  $LogEI(x)  = \\log(EI(x))$. Since $\\log$ is monontonic, we have that $\\arg \\max_{z \\in \\mathbb R_{>0}} z  = \\arg \\max_{z \\in \\mathbb R_{>0}} \\log(z)$. Hence, $\\arg \\max_{x \\in \\mathbb X, EI(x) > 0} EI(x) = \\arg \\max_{x \\in \\mathbb X, EI(x) > 0} LogEI(x)$. \n\nFigure 1 shows that LogEI holds approximately the same value as the logarithm of EI, except in regions where the log of EI is numerically zero.\n\n\n\n__Figure 2, Model quality, and constant vs. zero values__\n\nFirst, we’d like to clarify the data generating process (DGP) for the training data: Training points are not chosen uniformly at random, but rather, 80% are sampled uniformly at random from the domain, and 20% are sampled according to a MVN centered at the function maximum with a standard deviation of 25% of the length of the domain. The idea behind this DGP is to mimic the kind of data one would see during a BO loop (for illustration purposes, without having to run thousands of BO loops to generate the figure). We will clarify this in the CR.\n\nFigure 2 in the MT reflects the fact that under this data generating process with the chosen test problem, the incumbent (best observed point) is much better than the values at the random test locations, and this becomes increasingly the case as the dimensionality of the problem increases and the number of training points grows. For a particular replicate, Figure 1 in the attached PDF shows the model fits in-sample (black), out-of-sample (blue), and the best point identified so far (red) for our DGP with 60 training points and a random subset of 50 (out of 2000) test points. One can see that the model produces decent mean predictions for out-of-sample data, and that the uncertainty estimates appear reasonably well-calibrated (e.g., the credible intervals typically cover the true value). Because of this, we do not see any clear evidence that there is something odd about the model and its length scales.\n\nWhat Figure 1 in the attached PDF does show is that while there is ample uncertainty in the predictions of the model away from the training points, for the vast majority of points, the mean prediction is many standard deviations away from the incumbent value (the error bars are +/- 2 standard deviations). This is the key reason for EI taking on zero (or vanishingly small) values and having vanishing gradients.\n\nTo illustrate this, Figure 2 in the attached PDF shows the histogram of `z(x)` values, the argument to the function `h` in eq (2). It also contains the thresholds corresponding to the values `z` below which `h(z)` is less than the respective threshold. Since `sigma(x)` is close to 1 for most test points (mean: 0.87, std: 0.07), this is more or less the same as saying that `EI(z(x))` is less than the threshold. It is evident from the histogram that the majority of the test points fall below these threshold values (especially for larger thresholds), showing that the associated acquisition function values (and similarly the gradients) are numerically almost zero and causing issues during acquisition function optimization. \n\n__Relevance in high dimensions__\n\nLogEI tends to be particularly effective in high dimensions and we can see how the discussion of Lemma 1 did not adequately address this. \n\nIt is important to note that dimensionality alone is not the key factor at play, but the measure (“volume”) of points that attain significantly sub-optimal objective values in the search space (the right-hand side of Lemma 1). \n\nWhile the average posterior uncertainties tend to be larger in a larger proportion of the search space as  the dimensionality increases, they are usually much smaller than the empirical standard deviation of the objective values, and this means that the predictive mean can still be many predictive standard deviations away from the incumbent value for large swaths of the search space.\n\nFig. 1 and 2 in the attached pdf are an empirical validation of this intuition: while the posterior uncertainties are large and cover the errors, the corresponding EI values are vanishingly small.\n\n__Existing work__\n\nAs you already noted in your review, SMAC’s “LogEI” acquisition function tackles a fundamentally different problem than our LogEI. In the former, the surrogate model is fit to log-transformed outcomes in the hope that a log-transformation improves the model fit and optimization performance. Despite the nominal similarity, transformations of the outcomes before fitting the surrogate model are entirely orthogonal to this work, as they do not seek to solve the fundamental problem of the acquisition function itself being hard to optimize due to vanishing gradients. In fact, the techniques underlying our LogEI acquisition function are complementary to the outcome transformation (e.g. a log transform in SMAC’s LogEI) and would help resolve numerical issues and vanishing gradients. We will clarify these differences in the section on related work. See Fig. 4 in the PDF for a comparison of LogEI with SMAC's LogEI.\n"}, {"review_id": "mLlyzOwMQ4", "reviewer": "Reviewer_UCkN", "summary": "This paper addresses a major weakness of Bayesian expected improvement acquisition functions, which are ubiquitously used for black-box optimization tasks such as computational hyperparameter tuning, materials science, and biomedical research. It is very common to use gradient-based optimizers to find local maxima of the acquisition surface, however expected-improvement acquisition functions have the very unfortunate pathology of a completely flat acquisition surface (where both the acquisition value and acquisition gradient are 0) on large regions of input space, particularly as optimization progresses and the best known solution improves. This pathology makes Bayesian optimization extremely sensitive to implementation decisions, particularly the initialization scheme of the acquisition maximization subproblem, which hinders BayesOpt practitioners in academia and industry. This paper rightfully places numerical precision and stability as one of the primary considerations in acquisition function design, and proposes simple, intuitive modifications to expected improvement acquisition functions that significantly improve performance.", "questions": "Can you comment on whether you expect batch acquisition value optimization to outperform sequential greedy optimization when there are strong locality constraints placed on the inner loop problem (i.e. d(x_0, x_t) < \\varepsilon for all x_t optimization iterates)? Intuitively it seems that the performance of batch acquisition optimization once again comes down to heuristics for choosing the right collection of points as the initial solution (e.g. [2]), since the iterates may not be able to move far enough from the initialization for the improved acquisition landscape to make much difference.", "limitations": "I think the discussion section could be expanded a bit. In particular I think the following rather vague sentence could be made more specific:\n\n\"While our contributions may not apply verbatim to other classes of acquisition functions, our key insights and strategies do translate and could help e.g. with improving information-based [20, 42], cost-aware [26, 36], and other types of acquisition functions that are prone to similar numerical challenges.\"\n\nI take this to mean that this paper has primarily focused on resolving numerical difficulties arising from the use of the max operator, and other acquisition functions may have numerical issues from operators that are not the max. If this is the case I think it could be stated more clearly, as it would give readers a clearer picture of avenues for future work.", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "I am strongly in favor of accepting this paper. The basic problem the authors are addressing is one I have often encountered myself, and I have even tried some similar ideas as those presented in this paper to try to address the problem, however I had to shelve the project due to competing demands for my time. I'm delighted to see the problem addressed so thoroughly here.\n\nThe greatest strength of this paper is the emphasis placed on how acquisition function design interacts with the optimization algorithms used to find their maxima. Generally speaking I feel this aspect of acquisition function design is often neglected in many Bayesian optimization papers, to the great detriment of the field.\n\nGiven the widespread use of BayesOpt across industries, and the use of EI-style acquisition functions in particular, I think this paper could have significant practical impact on multiple fields.", "weaknesses": "This work is ready for publication without any significant revisions. I would encourage reviewers in general to think about the opportunity cost of burdening authors with minor or tangential concerns, slowing the development of follow-up work.\n\nIt's worth noting that the weaknesses of EI-style acquisition functions are fairly well documented in latent-space BayesOpt papers, such as [1] and [2]. Both of those works employed a heuristic I didn't see mentioned in the paper, which is to scale the max_{x_i \\in D} f(x_i) term in the acquisition function with some factor < 1 (e.g. 0.9). Some brief discussion in the related work on this point could help better communicate the potential impact of this paper.\n\n[1] Tripp, A., Daxberger, E., and Hernandez-Lobato, J. M. ´\nSample-efficient optimization in the latent space of deep\ngenerative models via weighted retraining. Advances in\nNeural Information Processing Systems, 33, 2020.\n\n[2] Stanton, S., Maddox, W., Gruver, N., Maffettone, P., Delaney, E., Greenside, P., & Wilson, A. G. (2022, June). Accelerating bayesian optimization for biological sequence design with denoising autoencoders. In International Conference on Machine Learning (pp. 20459-20478). PMLR.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your detailed review and encouragement.\n\n__Regarding the scaling of the incumbent__ \n\nThank you for pointing us to the heuristic of scaling the incumbent by a factor. We have encountered and experimented with this heuristic in the past, and while it can be employed to try to avoid numerical degeneracies, it has drawbacks compared to the solution proposed in the present paper:\n- For a fixed scaling factor, the resulting acquisition function – to the best of our knowledge –  does not have a principled grounding, requires setting a hyper-parameter that the optimization is sensitive to, and it is unclear which scaling factor is sufficient to remedy the numerical problems a-priori.\n- Using a homotopy approach to sequentially increase the scaling factor can be more robust, but this requires solving a sequence of acquisition function optimization problems, which is more computationally expensive than optimizing LogEI once. We experimented with this approach in the early stages of this project.\n\nLogEI successfully circumvents the downsides of the incumbent-scaling approach without adding the computational overhead of solving multiple acquisition function optimization problems. \n\nWe agree that discussing this will further improve the paper and the motivation of the methods. Will do so in the related work section. \n\n__Regarding greedy vs sequential batch optimization__\n\nWe are not sure we fully understand your comment regarding “strong locality constraints”, but believe it may refer to explicitly constraining the distance between new candidates and previously seen points in the context of latent space BO. Such a heuristic has been used in the literature to avoid exploring parts of the latent space that the decoder cannot map back well in the original input space. Assuming this is indeed the setting, it isn’t immediately clear to us whether these imposed locality constraints would have a significant effect on whether joint optimization is beneficial over sequential greedy. Certainly, as $\\varepsilon \\to 0$ we would expect the benefit to disappear (and similarly we would expect to recover it fully as $\\varepsilon -> \\infty$, assuming the decoder performance is good across the entire latent space), but it’s hard to say in general what the behavior would be. It’s an interesting question that deserves further study.\n\n__Regarding the extensibility of the methods__\n\nWe seek to highlight two distinct aspects:\n\n- Analytical LogEI relies primarily on careful and stable implementations of Gaussian log-probabiliby functions, their sums, differences, and fractions. We believe information-based acquisition functions like Gibbon and Joint Entropy Search (JES) could benefit from similar treatments. A particularly striking result of our experiments revealed that JES fails to outperform random search on the Michaelwicz function in greater than 8 dimensions, which could be caused by numerical issues of the acquisition function.\n\n- Monte-Carlo (“q”) LogEI primarily relies on novel smooth approximations to the ReLU non-linearity, as well as max operator, in addition to the associated log-tricks. We believe that especially the fat-tailed smooth operators (Appendix A.3) could find wider application even outside of Bayesian optimization, in the design and optimization of deep architectures, which are known to lead to saturating nonlinearities and vanishing gradients. We aim to explore this in follow up work.\n"}]

Author rebuttal:
We thank the reviewers for their detailed and predominantly positive reviews. We are attaching a one-page pdf with additional figures to help answer questions that arose during the review process, and are responding to each reviewer's questions in detail below and in the comments.

__Generality of the Methodology__

> The proposed treatment only works for the acquisition functions with vanishing gradient issues.

Our work brings to light the fundamental challenges of optimizing the popular EI family of acquisition functions (AFs), and proposes remediations tailored to the analytic, batch, constrained, and multi-objective setting, which cover extensive application domains. Through the lens of this popular family of AFs, we investigate how the formulation of AFs contribute to one’s ability to optimize them, and in general how experimental results may be subject to the quality of the AF optimization procedure. While we consider a broad class of EI-based acquisition functions in this work, we hope that it can inspire further investigation, thought, and potential remediation in the development of other types of acquisition functions, promising examples including entropy-based acquisition functions.

__Extension to Noisy Observations__

We point out that LogEI extends naturally to the Noisy Expected Improvement (NEI) acquisition function. We did not include NEI in our original manuscript since the construction is a trivial extension of qEI. The only difference between qEI and qNEI is that the computation of "best_f" values are based on samples of the GP at previously observed points, rather than taking the empirically observed objective value (cf. Balandat et al. NeurIPS 2020, S5.2). We then directly use the forward pass of qLogEI to compute qLogNEI. 

In the attached PDF (Fig. 3), we show empirical optimization performance of qLogNEI compared to qLogEI, EI, NEI, and Gibbon on Hartmann 6D, Ackley 8D, and Ackley 16D for varying noise levels. We set the noise level as a proportion of the total range of the respective function, which is ~3.2 for Hartmann and ~20 for Ackley. Thus, a noise level of 1% * Range(f) is equivalent to Gaussian noise with a std of 0.2 for Ackley. 

On Hartmann 6D, qNEI, qLogNEI, and Gibbon consistently find the optimum in the allocated number of iterations. qEI and qLogEI exhibit higher variance in their optimization traces especially at the highest noise level, which is expected since the `best_f` value they rely on becomes highly stochastic and might be far from the true (noiseless) best objective value corresponding to the queried inputs. 

Notably, qLogNEI outperforms both canonical EI counterparts and leads to significantly improved optimization on Ackley. qLogNEI also leads to notable improvements over Gibbon for the higher-dimensional functions. 

We will include a description of qLogNEI, include its implementation in the code release, and add a treatment of these results in the SM.

__Non-continuous search spaces__

Our work largely focuses on the problem of vanishing gradients, which is most pronounced for problems involving continuous or mixed spaces.  We expect the proposed methods to have less impact in fully discrete or mixed spaces, especially when some combinations of input parameters can be exhausted entirely. Nevertheless, since our work ensures that the acquisition function values (not just the gradients) do not numerically become zero, it will result in better optimization performance in settings where the feasible choices are all “far” from the incumbent. Promising approaches to gradient-based optimization of fully-discrete or mixed spaces with difficult-to-enumerate combinations such as Probabilistic Reparameterization or straight-through estimators (Daulton et al. NeurIPS 2022) may particularly benefit from LogEI.  We will add a discussion of this to the CR. 

__Conventional Benchmarks__

As a potential drop-in replacement for EI, seeing its performance of the method on the most conventional low-dimensional tasks (Branin, Hartmanns) would be informative. Moreover, it would be helpful for future benchmarking.

We included benchmarks that had flexible numbers of input dimensions, many constraints, or many outcomes to highlight how these factors impact the vanishing gradient problem and benefit from LogEI.  We agree that including such canonical test problems are worth documenting and will include them in the CR. In the interim, please see the attached PDF for the results on Hartmann and Branin (Fig. 4). LogEI outperforms alternatives on Hartmann (6 dimensional) and matches the performance of EI on Branin (2 dimensional).

__LogEI as Replacement of EI__

LogEI should work for any case where EI works, as the two acquisition functions have the same optima (in the analytic case) — LogEI just makes optimizing for those optima much easier. One could imagine cases where the model is incorrect or so degenerate that other model-free policies such as random search or evolutionary strategies work better.  To our knowledge, such concerns are true of most model-based acquisition functions studied in the literature. Otherwise, we confidently recommend LogEI as a drop-in replacement of EI. 
On a practical note, our implementations indeed already share the same interface as BoTorch’s EI variants, making such a replacement particularly convenient in one popular BO framework.


Meta-review: The reviewers and meta reviewer all appreciated the quality of the work with a clear, well-written manuscript. The paper focuses on a simple trick (a drop-in replacement) to improve the widely-used expected-improvement acquisition function. The reviewers and meta reviewer are convinced that this paper will inform future research and the work of practitioners. 

They thank the authors for their response and their efforts during the rebuttal phase, which was helpful to improve the submission and clarify concerns (e.g., relation to previous work). The reviewers and meta reviewer unanimously recommend the paper for acceptance.
Decision: Accept (spotlight)