Batched Energy-Entropy acquisition for Bayesian Optimization

Authors: felix  teufel
carsten  stahlhut
jesper  ferkinghoff-borg
machine  intelligence
novo  nordisk
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 14:08:51

Abstract: Bayesian optimization (BO) is an attractive machine learning framework for performing sample-efficient global optimization of black-box functions. The optimization process is guided by an acquisition function that selects points to acquire in each round of BO. In batched BO, when multiple points are acquired in parallel, commonly used acquisition functions are often high-dimensional and intractable, leading to the use of sampling-based alternatives. We propose a statistical physics inspired acquisition function that can natively handle batches. Batched Energy-Entropy acquisition for BO (BEEBO) enables tight control of the explore-exploit trade-off of the optimization process and generalizes to heteroskedastic black-box problems. We demonstrate the applicability of BEEBO on a range of problems, showing competitive performance to existing acquisition functions.

Main paper: 
 introduction: 
 Figure 1 : q-UCB does not allow for controlling its explore-exploit trade-off with large batches. A GP surrogate (background) was initialized with 100 random points of the Ackley function. q-UCB was run with κ = 0.1 and κ = 100, BEEBO with T ′ = 0.05 and T ′ =50. Batch size Q=100.Bayesian Optimization (BO) has since its inception [1, 2] made a profound contribution to the realm of global optimization of black-box functions through the usage of Bayesian statistics. For global optimization problems pursuing x * = argmax x∈X f true (x), BO has surfaced as a premier strategy for efficiently handling especially complex and costly unknown functions, f true (x). While BO is traditionally formulated in a single-point scenario, where individual points are queried and results are observed sequentially, there are situations where batched acquisition is needed. Such situations arise when f true (x) is expensive to evaluate in either time or cost, but can be effectively evaluated in parallel by dispatching multiple experiments, reducing the overall optimization time. This is often the case in e.g. drug discovery, materials design or hyperparameter tuning for deep models [3, 4, 5, 6, 7] .The realization that BO could be employed for the training of deep neural networks, as suggested by [3] , sparked renewed research interest, with advancements encompassing a variety of areas, including the generalization to accommo-date noisy inputs [8, 9] , heteroskedastic noise [10, 11] , multi-task problems [12] , multi-fidelity [13] , high-dimensional input spaces [14] , and parallel methods with batch queries [15, 16] . Generally, these desired properties are addressed by customizing one of the two key components in BO, either the surrogate model or the acquisition function. The surrogate model f approximates the black-box function f true using the available data. In BO, the surrogate is formulated from a Bayesian perspective, allowing us to quantify the model's uncertainty when evaluating new points. Typically, the model of choice is a Gaussian Process (GP) [17] . The acquisition function is responsible for guiding the selection of new input point(s) to evaluate at each optimization step, utilizing the surrogate model to identify promising regions in the input domain and exploring the unknown function further.Any acquisition process needs to trade off exploration (reducing uncertainty to learn a better surrogate model) against exploitation (selecting points with a high expected f true (x) based on the current surrogate). In this work, we are particularly interested in acquisition processes that make this trade-off controllable using a hyperparameter. Controllability can be a desirable property if e.g. domain knowledge relating to the difficulty of the optimization process and the quality of the surrogate model is available, or if the strategy needs to be adjusted depending on future experimental budgets. Similarly, it can be desirable to acquire multiple x with high f true (x) in a batch (as opposed to just finding the optimum x * , with the remaining x being considered explorative). This is useful when optima identified in BO can be subject to constraints that are unknown at optimization time, but may render x * intractable [18] . Such constraints arise when the f true explored in BO is a necessary simplification of the actual objective. Practical examples include e.g. the synthesizability of a material at larger scale, when BO experiments are performed at lab scale; or the in vivo activity of a molecule with BO experiments performed in vitro.A wide range of batch mode acquisition functions has been proposed, with approaches often leveraging random sampling strategies or Monte Carlo (MC) integration, which can adversely affect controllability for large batches (Figure 1 ). In contrast, we here introduce BEEBO (Batched Energy-Entropy acquisition for BO), a statistical physics inspired acquisition function for BO with GP surrogate models that natively generalizes to batched acquisition. BEEBO enables • Parallel gradient-based optimization of the inputs, without requiring sampling or Monte Carlo integrals. • Tight control of the explore-exploit trade-off in batch mode using a single temperature hyperparameter. • Risk-averse BO under heteroskedastic noise.We demonstrate the application of BEEBO on a wide range of test problems, and investigate its behaviour under heteroskedastic noise. 
related works: 
 Batch variants of traditional strategies Parallel acquisition in BO has seen a variety of approaches, often starting from established single-point acquisition functions like probability of improvement (PI), expected improvement (EI), knowledge gradient (KG) or upper confidence bound (UCB) [2, 19, 20, 21, 22] . Reformulating these to batch mode with Q query points, we obtain q-PI, q-EI, and q-UCB [23, 24] . While the single-point specifications provide an analytical form and enable gradient-based optimization, batch expressions are more challenging and require different optimization strategies, typically involving greedy algorithms [25] or deriving an integral expression over multiple points.For instance, in the popular EI acquisition function, a single point is selected by maximizing the expression a EI (x) = E[max (0, f (x) -f * t )] = max (0, f (x) -f * t ) P (f |x) df . Here f * t represents the best observed evaluation of f true so far. With a surrogate model in the form of a GP, the acquisition function depends only on the predictive mean and variance functions, µ (x) and C (x). Effectively, we need to evaluate the cumulative normal distribution, which quickly becomes intractable for large batch sizes and approximating the gradient of the q-EI acquisition function typically requires MC estimation [26, 27] . However, proper MC integration can be laborious and is sensitive to both the dimension of the problem and the choice of batch size Q. Specifically, MC methods face the curse of dimensionality problem when applied to high-dimensional integrals, as they require an exponentially increasing number of sample points to maintain accuracy, making them computationally impractical for such tasks [28, 29] . Of particular interest is Wilson et al. [24] , in which they adopt the reparameterization trick [30, 31] on acquisition functions integrals, enabling gradient based approaches to the optimization of PI, EI, and UCB. This demonstrates particular usefulness in modest to higher dimensions.While EI trades off exploration and exploitation, users do not have a direct control over the balance.To alleviate this, Sobester et al. [32] proposed a weighted EI formulation. An alternative strategy with an explicit explore-exploit trade-off is offered by the UCB acquisition function, a UCB (x) = µ(x) + √ κ • C(x), which directly expresses exploration and exploitation as two terms, traded off by the parameter κ. As we are particularly interested in enabling this direct user control, we focus our primary comparison on q-UCB in the main text, while a more extensive comparison with alternative methods can be found in Appendix B, both theoretically and experimentally.Greedy strategies As mentioned, a popular approach for leveraging single-point acquisition functions is devising batch filling strategies that score candidate points sequentially. Kriging Believer (KB) [33] uses EI to select points and iteratively updates the GP by fantasizing an observation with the posterior mean. Likewise, GP-BUCB [34] uses fantasized observations to update C(x) at each step. Local penalization (LP) [35] introduces a penalization function that repulses selection away from already selected points. Contal et al. [36] propose selecting a single point using UCB and dedicating the remainder of the batch budget for exploration in a restricted region around the believed optimum. GLASSES [37] treats batch selection as a multi-step lookahead problem to overcome the myopia of only considering the immediate effect of selecting a point.Entropy based strategies From an information theory perspective, BO can be interpreted as seeking to reduce uncertainty over the location of optima of the unknown function. This has given rise to entropy-based acquisition functions such as entropy search (ES) [38] , predictive entropy search (PES) [39] and max-value entropy search (MES) [13, 40, 41] . MES is distinct in that it seeks to quantify the mutual information between the unknown f true (x * ) and the observations y|D, rather than the location of x * . General-purpose Information-Based Bayesian OptimizatioN (GIBBON) [42] provides an extension of MES that enables application to batched acquisition as well as other challenges such as multi-fidelity BO. GIBBON proposes a lower bound formulation for the intractable batch MES criterion, which is then optimized using greedy selection. Despite being formulated to handle a large degree of parallelism, Moss et al. [42] reported that GIBBON fails in practice for large batches with Q > 50. Potentially, this behaviour is a consequence of the accuracy of the lower bound approximation. A heuristic scaling of the batch diversity was proposed to improve performance with large batches. GIBBON may also be interpreted as a determinantal point process (DPP) [4, 43] . In Appendix B we provide a detailed discussion of the relationship of the BEEBO acquisition function to GIBBON and DPPs. Note that while we will also make use of the term entropy in BEEBO, the quantity is distinct from the ones leveraged by the aforementioned approaches in the sense that it does not relate to an unknown optimum.Thompson sampling Given the challenges of generalizing acquisition functions to batch mode, Thompson sampling (TS), which was originally adopted from bandit problems [44, 45, 46, 47, 48] , is a popular alternative strategy for guiding batched BO. While being an attractive approach in general, it has been demonstrated that default TS can become too exploitative, motivating the use of alternatives such as Bayesian Quadrature [49] , or advanced strategies on top of TS that ensure diversity [18] . Eriksson et al. [50] demonstrate that overexploration also can be problematic in higher dimensions, and alleviate this using local trust regions in TuRBO. Maintaining such regions with high precision discretization can be memory-expensive, as indicated by [51] , who suggest using MCMC-BO with adaptive local optimization to address this by transitioning a set of candidate points towards more promising positions. 
the beebo acquisition function: 
 Assume f true : X → R is some real output associated with the input and a set of data be given D = {(x i , y i )} N i=1 where y i ∈ R represent some noisy observations of f true (x i ), sayEQUATIONwith ϵ i denoting the measurement noise. Let x = (x 1 , • • • , x Q ) ∈ X Q represent a collection of test points we wish to assign an acquisition value to. In keeping with the BO framework, we assume a given posterior probability distribution over the surrogate function f evaluated at x,EQUATIONThe lack of knowledge we have of the surrogate function at x is quantified by the differential entropy H:EQUATIONThis entropy can be contrasted with the expected entropy of the surrogate function if Q observations y = (y 1 , • • • , y Q ) were acquired at x, i.e. if the training data D would be augmented with D ′ (y) = {(x q , y q )} Q q=1 to form the joint data set D aug (y) = D, D ′ (y) . We refer to this entropy as H aug :EQUATIONwhere P (y | D, x) represent the posterior predictive distribution at x. The expected information gain, I(x), from acquiring observations at x is given by the expected reduction of entropy from this process:EQUATIONWe propose to represent the explore component of the acquisition function, a BEEBO , by I(x). The information gain I(x) is distinct from the quantities exploited by entropy search approaches, as it quantifies global uncertainty reduction, rather than estimating the information over an unknown x * . The information gain is directly applicable to multivariate functions and to heteroskedastic settings where σ 2 = σ 2 (x). Since large measurement uncertainties imply smaller information gain, a BEEBO exhibits risk-averse behaviour [11] by automatically prioritizing regions of small uncertainties from where more precise information of f true can be obtained, everything else being equal.The exploit component of BEEBO relies on taking expectation values of a scalar function of the random variable f (x), Ẽ : R Q → R, that summarizes the optimality properties of a given batch x. Natural choices would be the mean or the maximum of f (x). Of particular interest is expressing the optimality as a softmax-weighted sum over f (x), as this allows us to smoothly interpolate between the two regimes:EQUATION)where β is the softmax inverse temperature. At β = 0, we recover the mean. We scale the expectation with Q so that both I and E scale linearly with increasing batch size. While the mean provides a closed form expression for its expectation, this is not the case for the general softmax-weighted sum of a multivariate normal. Using Taylor expansion, we introduce an approximation of the expectation of the softmax-weighted sum that is fully differentiable and can be computed in closed form. A detailed derivation is provided in Appendix A. At β = 0, all Q points contribute equally to E(x), whereas at β > 0, points that do not compete for optimality are dynamically released. This effect can be quantified as the effective number of points via the entropy of the softmax weights. In the following, we will refer to the (exact) β = 0 limit as meanBEEBO, and the (approximated) general case as maxBEEBO.The BEEBO acquisition function then takes the formEQUATIONwhere T sets the balance between exploitation (small T ) and exploration (large T ). As both E and I scale with the batch size Q, a given choice of T would set the explore-exploit balance in an approximately Q-independent manner. This acquisition function bears a strong similarity to the definition of (negative) free energies in statistical physics, where E and I correspond to respectively the thermodynamic energy and entropy of the system and T corresponds to the temperature.Gaussian processes offer a particular convenient framework for BO, due to the availability of closedform expressions for the inference step [17] . SpecificallyEQUATIONwhere N (• | µ, C) is the multivariate Gaussian distribution with mean µ, and covariance C, x D and y D are the x and y values of the acquired data, σ(x D ) = diag σ 2 1 , • • • , σ 2N is a diagonal matrix with the measurement uncertainties in the diagonal and K(•, •) are matrices derived from the GP-kernel,k(•, •), i.e. K(x, x ′ ) ij = k(x i , x ′ j ).It is worth noting that C(x) only depends on the input location of the test points x and the data points x D with their corresponding measurement uncertainties, σ 2 (x D ), but not on the actual observations, y D . Consequently, the entropy of the posterior distributionEQUATIONis independent of y D as well, with ln det denoting the log determinant. Similarly, the expected entropy of f if observations at x were acquired, simply readsEQUATIONwhereEQUATIONand x aug = x Daug . The BEEBO acquisition function is then given byEQUATIONwhere the expectation is either the mean, 1 Q Q q=1 µ q , or the closed form approximation of the softmax-weighted sum described in Appendix A, andEQUATIONAlgorithm 1: meanBEEBO optimization Input: model GP, initial batch points x, temperature T repeat Calculate µ(x), C(x) from Equation 8using GP E ← -Q q=1 µ q GP aug ← fantasize(GP, x) Calculate C aug (x) from Equation 11 using GP aug I ← 1 2 ln det (C(x)) -1 2 ln det (C aug (x)) a ← -E + T * I x ← x + γ∇a until converged Output: optimized batch points xAll operations needed to compute the acquisition value a BEEBO (x) are analytical. Using automatic differentiation, the batch of points x can therefore be optimized with gradient-based methods, as laid out for meanBEEBO in Algorithm 1, with learning rate γ. In the pseudocode, GP denotes a trained GP model object that holds the training data and the kernel function. Using the kernel's learned amplitude A, we can relate BEEBO's T parameter to the κ of UCB. This allows us to configure BEEBO using a scaled temperature T ′ that ensures both methods have equal gradients at isosurfaces, enabling the user to follow existing guidance and intuition from UCB to control the trade-off. A derivation is provided in Appendix B.1. We benchmark acquisition function performance on a range of maximization test problems with varying dimensions (Table 1) available in BoTorch [52] . Test problems that are evaluated on multiple dimensions support specifying the respective arbitrary d. As a high-dimensional problem with low inherent dimensionality, we embed the six-dimensional Hartmann function in d = 100 [50, 53, 54] . We additionally test on two robot control problems (robot arm pushing and rover trajectory planning) in Appendix D.3 [55, 56] .On each test problem, we perform 10 rounds of BO using q-UCB or BEEBO with a given explore-exploit parameter for direct comparison. We use the scaled temperature T ′ (B.1) to ensure that both methods operate at the same trade-off. In round 0, we seed the surrogate GP with Q random points that were drawn so that each point has a minimum distance of 0.5 to the test problem's true optimum. We perform ten replicate runs for each problem and method, with replicate seeds controlled so that all methods start from the same Q random points in a replicate. As we evaluate performance in a fixed-round, fixed-Q optimization scenario, we set the explore hyperparameter to 0 in the last round (for maxBEEBO, we also set the softmax β to 0). We use Q = 100 for all experiments, which is commonly understood to be a large batch size [50] . Additional results on small batch sizes (5, 10) are provided in Appendix D.2. All experiments use BoTorch's default utilities for acquisition function optimization and GPyTorch [57] GP training (C.1).We investigate performance when optimizing under heteroskedastic noise on the 2D Branin function with three global optima. To construct a heteroskedastic problem, we specify noise so that the noise level is maximal at optima 2 and 3, decaying exponentially with distance to any of the two noised optima (C.4). No noise maximum is added at optimum 1. Therefore, while all three optima share the same f true (x) (Figure A1 ), only optimum 1 is favorable in terms of heteroskedastic risk. We perform BO for ten rounds with β = 0.1 and Q = 10 using a heteroskedastic GP that learns surrogate models for both f true (x) and σ 2 (x). We report results over five replicate runs.Metrics We report the mean best observed objective value after 10 rounds over the five replicates. As test problems have highly varying scales, we normalize the results on each test problem using min-max normalization. Typically, the minimum of a maximization problem is not known explicitly. We therefore set the minimum for normalization to the highest value observed among the random seed points. The maximum is given by the f true (x * ) of the problem. The metric thus directly quantifies how much progress has been made to the true optimum from the random starting configuration on a 0-1 scale.As we are not only interested in identifying a single x with good f true (x), we additionally quantify the overall quality of the final (exploitative) batch. We compute the batch instantaneous regret R = q<Q f true (x * ) -f true (x q ) of the last, exploitative batch. To bring results on all test problems to a similar scale, we divide it by the batch instantaneous regret of a batch of Q random points on each problem. We refer to this metric as the relative batch instantaneous regret, R rel = R t=10 /R random .For BO under heteroskedastic noise, we wish to quantify the preference of a given method for different optima. As optima share the same f true (x * ), metrics operating on f true (x) are inherently unsuitable, and preference needs to be evaluated on x directly. For each acquired point x i , we compute the distances ∥x i -x * j ∥ 2 to the J individual optima. We report the mean distance to each optimum over all points in a batch.  
bo on test problems: 
 We benchmark BEEBO against q-UCB at three rates of κ. Overall, we find that the simpler mean-BEEBO variant outperforms maxBEEBO in terms of mean performance on all but the lowest rate of κ (Table 2 ). As we consider the configuration with the lowest rate to be exploit-dominated, this can be understood as a consequence of maxBEEBO effectively releasing non-contributing points for further exploration, with the low explore rate seeming sufficient to induce the necessary diversity.While results on individual test problems vary, meanBEEBO shows improved performance over q-UCB especially in the medium dimension range up to 50. For d=100, we find mixed performance, with meanBEEBO gradually becoming more competitive with increasing κ. This is due to the fact that with increasing dimensionality, more exploration is beneficial for learning a good surrogate model before an actual BO process becomes effective. As we find that q-UCB inherently performs more random-like sampling at large Q, irrespective of κ, it benefits in such situations.On average over all 33 performed experiments, meanBEEBO improves upon q-UCB for large batches at any of the three rates (Table A3 ). We additionally benchmarked BEEBO against other popular BO strategies without explore-exploit hyperparameters. Interestingly, we found that the Kriging Believer (KB) iterative heuristic [33] can perform very competitively for large batches when using LogEI [58] as the acquisition function (Table A1 ), especially on the two robot control problems (Appendix D.3), but can be slower to optimize than BEEBO (Table A9 ).When evaluating R rel in the ultimate round of BO, we find that BEEBO allows us to effectively acquire a batch with high f true (x), highlighting the controllability of the acquisition function (Table 3 ). The R rel of q-UCB is only slightly better than the R of a random batch in many cases, even though the explore component was explicitly set to 0. We note that this is not due to the surrogate function being unsuitable -the results in known by round 10. Rather, we assume that this a consequence of the challenges of MC-based optimization of the acquisition function at large Q.We compare performance of meanBEEBO and q-UCB on the 3-optimum Branin function. Under heteroskedastic noise, we find that BEEBO preferentially optimizes towards the low-noise optimum 1 at the expense of the noisy optima 2 and 3 (Figure 2 ) and is therefore risk-averse. In the homoskedastic case, BEEBO does not exhibit this preference and optimizes for multiple optima. As expected, q-UCB, which only uses the model posterior variance C(x) instead of quantifying the actual information gain, does not display any preference for low-noise optima, showing similar behaviour under heteroskedastic and homoskedastic noise and remaining risk-neutral. 
discussion: 
 We introduce BEEBO, an acquisition function for BO with GPs that can be optimized analytically and that scales natively to batched acquisition. By exploiting the independence of the information gain I(x) on measurements y when using GP surrogates, BEEBO models the interdependence of unknown points x in a batch and can optimize their positions jointly using gradient descent.BEEBO enables full control of its explore-exploit trade-off using a hyperparameter T that directly balances two terms, akin to UCB. Unlike in the reparametrization-based q-methods, BEEBO's T has predictable behaviour also at increasing batch sizes.The numerical complexity of BEEBO is dominated by the need to compute the inverse of M aug in Equation 11, which in a plain implementation scales as O((N + Q) 3 ). However, this can be reduced to O(N 2 Q); specifically, the Cholesky decomposition of M aug can be expressed as Q rank-1 updates of the pre-computed Cholesky decomposition of M D , where each update will have the complexity of O(N 2 ). The calculation of the energy, E, and the information gain, I, scales as O(N • Q) and O(Q 3 ), respectively, and are thus sub-dominant to the update needed for M -1 aug . For large N this approach may nevertheless become prohibitively slow. To overcome this limitation, methods for scalable GPs and fast predictive covariances such as LOVE [59] can be considered. The LOVE method allows a further reduction of the complexity of the Cholelsky update ofM aug to O(N • r • Q) [60], where r is the rank of the LOVE approximation for M D , typically r ≪ N .As opposed to q-UCB, BEEBO can take heteroskedastic noise into account when computing the information gain, and preferentially acquires more informative low-noise points. We note that when the noise function σ 2 (x) is unknown, and needs to be explored at the same time with f true (x), it is critical that the initial random points sufficiently capture the noise landscape well enough for the information gain component to be useful, as the uncertainty of the surrogate on σ 2 (x) is not used. This would require a fully Bayesian approach that integrates over the distribution of σ 2 (x). The problem does not arise if the heteroskedastic noise of an experiment is known beforehand by e.g. instrument calibration. While not the focus of this work, we note that using the information gain could also be beneficial in sequential single-sample BO on heteroskedastic problems. 
outlook: 
 In our experiments, we have focused on maintaining consistent explore-exploit ratios throughout the optimization rounds to ensure an equitable experimental comparison with q-UCB and demonstrate the effect of the hyperparameter choice. However, a more dynamic approach involving variable ratios could be more effective in real-world applications with a predetermined number of rounds [61] . Adopting a fully Bayesian perspective, one could consider the temperature hyperparameter T as a random variable. This opens up an intriguing avenue for BEEBO, where T could be drawn from a prior distribution that e.g. varies across optimization rounds, depending on the specific application. By tailoring this distribution, one could encourage a high level of exploration in the initial rounds, gradually transitioning towards a more exploitation-focused approach towards the end.In the presented experiments, we have implemented this as a strict constraint, maintaining a fixed T until the final round, at which point we shift to full exploitation, i.e., T = 0.While not explored in this work, we note that the BEEBO expression could naturally be extended to multi-objective optimization problems by capitalizing on GPs that handle vector-valued functions, such as multi-task GPs [12, 62] . Through e.g. the usage of the intrinsic model of coregionalization, we obtain a covariance function k, and thereby a covariance matrix C(x), over all input-task pairs. As the multi-task covariance matrix is jointly Gaussian, the expression of the information gain remains unchanged and can be computed like in the single-task case. The energy E(x) becomes vector-valued, providing an energy term for each of the tasks. This would allow for the introduction of task-specific weights in the acquisition function. As the extension only affects the surrogate model, the scaling remains cubic in the number of input-task observations.Beyond GPs, BEEBO could be generalized to work with any probabilistic model. However, GPs are unique in that H aug is available in closed form and can be used to compute I(x) analytically, without solving the integral over y in Equation 4. Other models may require approximations and sampling-based approaches for computing the information gain. We are free to choose any energy function Ẽ in BEEBO, the only requirement being that we are able to compute an expectation of Ẽ in order to obtain the scalar summaryE = E[-Ẽ(f )].Of particular interest is the softmax weighted sum,EQUATION)where β is the softmax inverse temperature. The softmax weight vector ω computed asω i = exp(βf i ) Q j=1 exp(βf j ) + exp(β y y max ) (A2)where exp(β y y max ) is an optional reference threshold value, as in expected improvement, which we set to 0 if not used (β y is either simply β or a dynamically scaled value that ensures Ẽ(f ) does not become 0, see Equation A28). The parameter β allows us to interpolate between two extreme regimes,Ẽ(f ) = 1 Q Q i=1 f i forβ → 0 (A3) Ẽ(f ) = max(f ) forβ → ∞ (A4)In the first regime, all points of a batch equally contribute to the energy, whereas in the second regime only the single point "responsible" for the maximum is controlling the energy. Note that for numerical reasons, operating towards the β → ∞ limit is impractical, as it will lead to zero gradients for all but one point, preventing optimization. We can set β = A -1/2 , with A being the prior uncertainty scale of the GP kernel. Since A represents the expected energy fluctuations for points far from data, this weighting scheme will reflect a natural compromise between Equation A3 and Equation A4.As opposed to the mean, in the general case, the expectation of the maximum of a Q-dimensional multivariate normal is not available in closed form. To our best knowledge, this is also the case for the softmax weighted sum. In the following, we derive a closed-form approximation of the expectation of the softmax of βf that can be used for gradient-based optimization.Consider the softmax denominatorEQUATIONWe will Taylor expand ln(d) to the second order, usingEQUATIONwhere δ ij is the Kronecker delta. SoEQUATION)where a is the Taylor expansion point, d a is d evaluated at a, ∆f = fa, w is the Q-dimensional vector ω evaluated at a, and W is theQ×Q matrix W = diag(w)-ww T . Inserted into Equation A2 we have ω(f ) i ≈ w i * exp(β∆f i -βw T • ∆f - β 2 2 ∆f T • W • ∆f ) (A10)With this approximation we can calculate the expectation valueEQUATIONSoftmax When y max grows much larger than the softmax input vector f -a situation that can arise easily when initializing with random points for gradient-based optimization -the softmax weights ω can become numerically zero for all "real" points, thus leading to E(x) = 0, and vanishing gradients.As we always wish to preserve a minimal energy contribution from the real points, we parametrize the inverse temperature applied to y max , β y , using a hyperparameter α that denotes the minimal fraction of probability mass pertaining to real points. This parametrization resembles the LogEI version of the expected improvement acquisition function [58] to address the problem of vanishing EI-gradients.Let N denote the softmax denominator excluding y max , N = Q j=1 exp(β∆f i ). We defineexp(β y ∆y max ) = min 1 -α α N, exp(β∆y max ) (A28)We used α = 0.05 as a default in all our experiments.We can interpret the softmax as the number of effective points contributing to the energy of the batch.The entropy H of the softmax is given byEQUATIONand the number of effective points, D eff , is exp(H(ω)), so thatEQUATION)D eff is bounded by 1 (approaching the maximum) and Q (approaching the mean). Note that if we include y max in the softmax denominator, we add ω ymax ln(ω ymax ) to H(ω), and the resulting number becomes bounded by 1 and Q + 1.In the following section, we will discuss how BEEBO is related to UCB, GIBBON, Determinantal Point Processes (DPP), the Local Penalization heuristic and RAHBO. We will base our analysis on meanBEEBO, as the softmax-mediated interdependency of points in maxBEEBO prevents a simple interpretation of the objective in a single-point stepwise manner and does not allow for the same direct analogies to other strategies.BEEBO bears some resemblance to the UCB acquisition function, which in the single particle mode,EQUATIONA31) where the parameter κ controls the balance between exploitation and exploration and µ(x) and C(x) are respectively the mean and variance of the posterior distribution, P (f | x, D), as before. We note that a UCB does not account for the uncertainty of the measurement at x, and therefore remains risk-neutral under heteroskedastic noise [11] . To understand the relationship between BEEBO and UCB, we will therefore limit ourselves to the homoskedastic case and furthermore assume that measurement variances σ 2 are much smaller than the typical prior variance of the GP surrogate, A, of f , e.g. A ≃ N -1 T r(K), so σ 2 ≪ A and M -1 = (K + σ 2 ) -1 ≈ K -1 . In this limit, the variance of f (x) after measurement (indexed at i = n, say) reduces to σ 2 :C(x) = (K -1 + σ -2 I) -1 nn = K(K + σ 2 I) -1 σ 2 nn ≈ σ 2 (A32)and the information gain becomesI(x) ≈ 1 2 ln(C(x)) -log(σ).Consequently, the gradient of the two acquisition functions reads∇a UCB (x) = ∇µ(x) + √ κ 2 • C(x) • ∇C(x) ∇a BEEBO (x) = ∇µ(x) + T 2 • C(x)• ∇C(x).The two gradients will be identical at points x where the posterior uncertainties satisfy C(x) = T √ κ . For comparison, we may desire equal gradients at iso-surfaces corresponding to a given fraction, ν, of the prior uncertainty scale √ A, by setting T accordingly asT = ν • √ A • √ κ.In our experiments, we use ν = 1 2 and configure BEEBO using a dimensionless T ′ explore-exploit parameter, defined asT ′ = T √ A , and set T ′ = 1 2 √κ for a given benchmark experiment.GIBBON [42] approximates the (intractable) General-purpose max-value Entropy Search acquisition function, which quantifies the mutual information M I(f * true ; y|D) of a batch of measurements y and the unknown optimum f * true . It does so using a lower bound on the information gain and MC estimation of the expectation over f * true . It can be written asEQUATION)where R is the correlation matrix with entriesR ij = C(x)ij √ C(x)iiC(x)jj, M is a set of samples for the max-value f * true , and ρ i is the correlation of y i and f true (x i ). ϕ and ϕ are the standard normal cumulative distribution and probability density functions, and γ i (m) = m-µ σ . The definition of BEEBO introduced in Equation 7, with the scalar summarization function set to the expected mean,E(x) = -1 Q Q i=1 µ(x i ), gives α BEEBO (x) = T * 1 2 (ln det (C(x)) -ln det (C aug (x))) + Q i=1 µ(x i ). (A34)From the second formulation of GIBBON, it becomes obvious that although being distinct in their motivation and derivation, BEEBO and GIBBON implement acquisition functions with a similar structure. Taking an information theoretic and multi-fidelity BO standpoint, GIBBON refers to this trade-off as diversity against quality, whereas in BEEBO we follow the intuitions of UCB, and use exploration and exploitation.• Quality -Exploitation: GIBBON employs an MC estimate of the lower bound approximation of the information gain provided by each point, whereas BEEBO directly summarizes the optimality of all points in closed form, either as their mean or an approximated softmax weighted sum. • Diversity -Exploration: In GIBBON, the diversity derived from the differential entropy H(f |D, x) is the entropy of the posterior correlation 1 2 ln det(R). In BEEBO, we employ the reduction of entropy, the information gain I(x). Under homoskedastic noise,I(x) ∝ ln det(C(x)). Since R(x) = diag(C(x)) -1 /2 • C(x) • diag(C(x)) -1 /2 , we have that ln det(R) = ln det(C(x)) - Q i ln(C(x) ii ).Therefore, maximizing the log determinant of R penalizes points that have high variance. Therefore, while GIBBON presents an attractive approximation of max-value Entropy Search for batched acquisition, BEEBO is an alternative that avoids approximating a quality criterion using MC. Moreover, GIBBON's diversity criterion implicitly penalizes points that have high variance, whereas BEEBO's criterion maximizes the reduction of variance. We find that BEEBO is orders of magnitudes faster to compute than GIBBON (Figure A3 ).In the context of large batches (Q >> 10), a modification of GIBBON exists that is further similar to BEEBO. Departing from the strict max-value entropy search derivation, a scaling factor Q -2 is introduced to counteract a growing dominance of the diversity term:EQUATION)This scaling is motivated by the fact that R contains Q 2 elements. However, we note that R is summarized by its log determinant, which scales linearly in Q: As the determinant is the product of the eigenvalues, the log determinant is the sum of the log-eigenvalues. The number of eigenvalues scales linearly with matrix size Q, and so does the log determinant.A Determinantal Point Process [43] specifies a probability over a set of points, or a "configuration of points" drawn from a ground set. Specifically, the probability of a set of Q points x is given byEQUATIONwhere L x is a Q × Q symmetric matrix. Kulesza et al. [43] provide a decomposition of the general DPP kernel L that makes quality and diversity components explicit, so thatEQUATIONwith k being a R d × R d → R + similarity kernel, and q being a unary R d → R scalar quality function. This framework is naturally amenable to batch BO, as we seek to select a collection of points that trade off quality (optimality) and diversity. Note that both k and q are distinct functions that need to be specified by the user, leading to the practical complication that they must be chosen very carefully so that their scales do not dominate each other, which limits the utility of this decomposition in practice [42] .In the following, we show how BEEBO is equivalent to a DPP, and derive the necessary k and q. Again, we consider BEEBOEQUATIONwith the scalar summarization function set toE(x) = - Q i=1 f (x i ).We will first focus on the information gain term I(x), which we can rearrange asI(x) = 1 2 ln det(C(x)) - 1 2 ln det(C aug (x)) = 1 2 ln det C(x) • C -1 aug (x) . (A39)Our similarity kernel k is therefore given by the entries of the matrix S = C(x) • C aug (x) -1 , so that k(x i , x j ) = S ij . Note that due to the augmented covariance term, the implied k also depends on all other currently selected points in x, and L x is not a submatrix of an all-sample L. Therefore, BEEBO does not implement a DPP under heteroskedastic noise. However, if we only consider homoskedastic noise, BEEBO's I(x) simplifies to the posterior entropy [63] , and therefore S = C(x). As C(x) can be accessed as a submatrix of an all-sample C, this permits a DPP.Given the choice of E(f ), we can rewrite BEEBO asα BEEBO (x) = ln det (S) * T * 1 2 + Q i=1 µ i 2 T * α BEEBO (x) = ln det (S) + Q i=1 2 T * µ i 2 T * α BEEBO (x) = ln det (S) + ln det (D)) with D = diag(exp( 2 T * µ)) 2 T * α BEEBO (x) = ln det D 1 2 • S • D 1 2 D 1 2 ii = exp( 2 T µ i ) = exp( 1 T µ i ) α BEEBO (x) = ln det D 1 2 • S • D 1 2 * T * 1 2 α BEEBO (x) = ln det (L) * T * 1 2 (A40)where L is a matrix with entriesL ij = S ij exp( 1 T * µ i ) exp( 1 T * µ j ).BEEBO therefore uses the DPP quality function q(x i ) = exp( 1T * µ i ), and, like proven previously for GIBBON, a batch x with maximal α BEEBO corresponds to the MAP of a DPP.Local penalization (LP) is a greedy batch selection strategy that given any arbitrary single-point acquisition function, ensures diversity by applying a penalization function ψ(x, x ′ ) that downweights the acquisition value of candidate locations x based on their proximity to already selected points. The criterion for selecting x i is given byx i = arg max α(x) i-1 j=1 ψ(x, x j ). (A41)Note that in this formulation, the product includes all previously selected points, not just the current batch. The penalization function ψ may in principle be chosen freely. Gonzalez et al. [35] propose exploiting the fact that f true is Lipschitz continuous in order to bound the position of the unknown optimum and penalize accordingly. The Lipschitz constant L is inferred from the GP surrogate and used to parametrize ψ. In LP, acquisition function optimization proceeds iteratively. After an x i is chosen, the corresponding penalizing multiplier is added to the objective before optimizing for the next x i+1 .While BEEBO enables optimization to proceed in parallel, it is of course possible to also optimize BEEBO greedily (under homoskedastic noise, I is submodular). In this case, it implements an LP strategy where α(x) = µ(x). Rather than a product of individual R d × R d → R function evaluations, the penalizer implied by BEEBO is the information gain I(x) : R i×d → R that we evaluate by concatenating a candidate point to the already acquired x at each iteration. Like in GIBBON, this constitutes an LP strategy that does not require estimation of any properties of f true beyond learning the GP surrogate.Risk-averse Heteroskedastic Bayesian Optimization (RAHBO) [11] is a UCB-derived single-point acquisition function that avoids heteroskedastic risk, preferentially selecting points with low noise. While it is not applicable to batched acquisition directly, we here compare it to single-sample BEEBO to highlight different ways of addressing noise. Given a heteroskedastic surrogate model that learns an additional GP for the noise, the variance proxy, RAHBO readsEQUATIONwhere µ f and σ f are the posterior mean and variance of the surrogate model and β f is the standard UCB trade-off hyperparameter, yielding the standard upper confidence bound UCB f . α is the chosen risk tolerance, and LCB is the lower confidence bound of the variance GP with posterior mean µ var and variance σ var traded off using β var .At Q = 1, BEEBO can be expressed asEQUATION)where the variance proxy at x is considered via the augmented posterior variance σ aug f . While RAHBO penalizes risk on an absolute scale, subject to α, BEEBO optimizes for high uncertainty reduction, quantified as the log ratio of the variance before and after making measurements.Moreover, RAHBO differentiates between known and unknown variance proxies, and uses the LCB var term to discount the predicted variance according to its uncertainty. In its closed-form analytical expression, BEEBO does not permit for the uncertainty of the variance proxy to be taken into account, being more similar to the known variance RAHBOα RAHBO (x) = µ f (x) + β f * σ f (x) -αµ var (x) (A44)where µ var is a noise-free proxy. Either a sampling-based approach, or approximations to I(x) would need to be introduced to handle variance proxy uncertainty in BEEBO.BEEBO was implemented for full compatibility with the BoTorch framework (version 0.9.4) [52] as an AnalyticAcquistionFunction. Standard BoTorch utilities for initializing and training GPs, initializing q-batches and performing gradient descent optimization of the acquisition function are used. We trained GPyTorch (version 1.11) [57] GP models with KeOps [64] Matérn 5/2 kernels (following BoTorch defaults with a separate length scale for each input dimension, and Gamma priors on the length and output scales). Log determinants for the information gain were computed using singular value decomposition for numerical stability.GPyTorch provides a get_fantasy_model method that allows for the efficient augmentation of the training data of a GP with a set of points, as done in BEEBO. However, we observed that GPy-Torch's implementation suffers from GPU memory leaks when used with automatic differentiation enabled. We therefore instantiate augmented models explicitly, not making use of the (more efficient) augmentation strategy.All experiments were performed with double precision. SobolQMCNormalSampler was used for acquisition functions making use of the reparametrization trick. Experiments were run on individual Nvidia RTX 6000 and V100 GPUs. Five replicates for the benchmarking experiments required a total of approx. 5,000 RTX 6000 GPU hours, with the majority of the run time dedicated to the GIBBON baseline, rather than BEEBO itself (Figure A3 , Table A9 ).All methods were benchmarked in BoTorch. For q-EI, we used LogEI [58] . For TS, 10.000 base Sobol samples were drawn and sampled with MaxPosteriorSampling using the Cholesky decomposition of the covariance matrix. GIBBON was optimized using sequential optimization following the BoTorch tutorial. We additionally implemented a custom version of GIBBON that applies the Q -2 scaling factor to the diversity term, as proposed in GIBBON's supplementary material. We used 100,000 random discretized candidates for max-value sampling. In a few iterations, optimizing GIBBON seemed challenging, with BoTorch reporting that no nonzero initialization candidate could be identified. KB was optimized using a custom greedy optimization loop with fantasized observations, using (single-sample) LogEI as the underlying acquisition function. TuRBO-1 was optimized following its BoTorch tutorial. None of the methods use a hyperparameter for controlling their explore-exploit trade-off. The results are therefore based on 10 iterations at defaults.Test functions All test functions were used in their BoTorch implementations. As done in previous work, the embedded Hartmann function was created by appending all-0 dummy dimensions to the original six dimensions [53, 54, 50] .We consider two control problems from previous work: A 14-dimensional parameter tuning task for controlling robot arms pushing two objects to a target location [55] , and a 60-dimensional trajectory planning task for a rover navigating through a maze of obstacles [56] .Instead of converting the problem objectives into rewards as in the original work, we operate on the actual minimization objectives directly (distance to target, navigation loss), and follow BoTorch's approach of simply inverting the objective in order to yield maximization problems. Both problems were adapted from their available implementations in Wang et al. [56] to follow the BoTorch test problem API.The (inverted) Branin function has three global optima f (x * ) = -0.397887 at x * 1 = (9.42478, 2.475), x * 2 = (-π, 12.275) and x * 3 = (π, 2.275). We define heteroskedastic noise so that the variance is maximal at x * 2 and x * 3 . The noise decays exponentially with the distance from any of the two noised optima at a rate λ.σ 2 (x) = σ 2 max * exp(-λ * min(∥x -x * 2 ∥ 2 , ∥x -x * 3 ∥ 2 ) (A45)For our experiments, we set σ 2 max = 100 and λ = 0.05. As the surrogate function, we use a HeteroskedasticSingleTaskGP provided in BoTorch. This model learns two GPs simultaneously, one for the function f (x) and one for the (also unknown) variance function σ 2 (x). When querying the oracle with a batch of points, noised observations of f (x) are provided together with the true σ 2 at each point. The homoskedastic control experiment uses a SingleTaskGP with inferred noise level. The homoskedastic noise is set to σ 2 = 77.5, which is the average noise level of the heteroskedastic function over the whole domain. T ′ =0.05 T ′ =0.5 T ′ =5.0 T ′ =0.05 T ′ =0.5 T ′ =5.0 √ κ=0.1 √ κ=1.0 √ κ=10.0 - - - default scaled - Ackley 2 0.993 ± 0.005 0.985 ± 0.031 0.975 ± 0.035 0.982 ± 0.023 0.980 ± 0.035 0.988 ± 0.013 0.973 ± 0.023 0.967 ± 0.022 0.988 ± 0.011 0.987 ± 0.012 1.000 ± 0.000 0.981 ± 0.014 0.878 ± 0.100 0.951 ± 0.027 0.951 ± 0.027 Levy 2 1.000 ± 0.000 0.999 ± 0.001 0.999 ± 0.001 1.000 ± 0.000 1.000 ± 0.000 0.998 ± 0.002 1.000 ± 0.000 1.000 ± 0.000 0.998 ± 0.002 1.000 ± 0.000 0.999 ± 0.002 1.000 ± 0.000 0.988 ± 0.008 0.993 ± 0.010 0.993 ± 0.010 Rastrigin 2 0.981 ± 0.024 0.989 ± 0.016 0.983 ± 0.016 0.993 ± 0.007 0.983 ± 0.011 0.993 ± 0.006 0.951 ± 0.021 0.983 ± 0.015 0.933 ± 0.025 0.995 ± 0.007 1.000 ± 0.000 0.976 ± 0.021 0.903 ± 0.087 0.944 ± 0.038 0.944 ± 0.038 Rosenbrock 2 0.976 ± 0.045 0.956 ± 0.071 0.955 ± 0.080 0.982 ± 0.032 0.979 ± 0.027 0.938 ± 0.123 0.949 ± 0.074 0.943 ± 0.129 0.962 ± 0.079 0.982 ± 0.029 0.966 ± 0.079 0.976 ± 0.068 0.633 ± 0.355 0.843 ± 0.301 0.843 ± 0.301 Styblinski-Tang 2 0.961 ± 0.072 1.000 ± 0.000 1.000 ± 0.000 1.000 ± 0.000 1.000 ± 0.000 1.000 ± 0.001 1.000 ± 0.000 1.000 ± 0.000 0.999 ± 0.001 1.000 ± 0.000 1.000 ± 0.000 1.000 ± 0.000 0.996 ± 0.003 0.999 ± 0.001 0. maxBEEBO Justification: We include full BO curves with standard deviations over five replicates for all quantitative experiments in the appendix. These detailed curves are referenced in the main text at the appropriate place.T ′ =0.05 maxBEEBO T ′ =0.5 maxBEEBO T ′ =5.0 Problem d q-UCB q-EI TS KB GIBBON GIBBON (s) TuRBO q-UCB q-EI TS KB GIBBON GIBBON (s) TuRBO q-UCB q-EI TS KB GIBBON GIBBON (s)T ′ =0.05 T ′ =0.5 T ′ =5.0 T ′ =0.05 T ′ =0.5 T ′ =5.0 √ κ=0.1 √ κ=1.0 √ κ=10.0 - - -T ′ =0.05 T ′ =0.5 T ′ =5.0 T ′ =0.05 T ′ =0.5 T ′ =5.0 √ κ=0.1 √ κ=1.0 √ κ=10.0 - - -T ′ =0.05 T ′ =0.5 T ′ =5.0 T ′ =0.05 T ′ =0.5 T ′ =5.0 √ κ=0.1 √ κ=1.0 √ κ=10.0 - - - -EQUATIONGuidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We list the used hardware and total GPU hours in the supplement and provide example timings for experiment runtimes.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Justification: The paper does not make use of human participants or datasets. To the best of our understanding, there are no potential harmful consequences and wider negative societal impact expected from the proposed method.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: The paper introduces a method for Bayesian optimization (BO). While BO has widespread applications in the sciences and engineering, there is no direct societal impact expected from this contribution. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: The paper does not introduce any trained models or novel data. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?38th Conference on Neural Information Processing Systems (NeurIPS 2024). 
 
 Appendix:  


Decision: This work introduces a new acquisition function inspired by statistical physics for batch Bayesian optimisation with GPs. The work includes a theoretical analysis and an extensive set of experiments to support the claims. While reviewers raised concerns regarding missing baselines and other experimental aspects, all were addressed during the rebuttal. Overall, this work is well executed and of interest to the community.

Reviews: review: paper_summary: paper_summary: This work introduces a batched acquisition function that balances exploration and exploitation by using a weighted sum of mutual information and expected value, with the weights defining the trade-off. The discussion links the proposed algorithm to UCB and asserts that it naturally addresses heteroskedastic noise.
strength_weakness: strength_weakness: 1. The introduced parameter controlling the trade-off lacks interpretation as in previous methods.
2. The completeness of the related work discussion is concerning. This is potential because the summarization lacks high-level extraction of the design of the algorithm, and the focus of the paper is, to some extent, scattered。
questions: questions: One concrete example of the second aforementioned weakness is the criticism of MC approximation in high-dimensional applications. Recent advancements in applying MCMC address this issue in a principled manner and might be of interest. 

***Reference:***

Yi, Zeji, Yunyue Wei, Chu Xin Cheng, Kaibo He, and Yanan Sui. "Improving sample efficiency of high dimensional Bayesian optimization with MCMC." arXiv preprint arXiv:2401.02650 (2024).
limitations: limitations: Discussed above.
paper_summary: paper_summary: The paper introduces a new approach to batch Bayesian optimization that explicitly trades off between exploration and exploitation via energy and entropy terms. The method is able to efficiently generate large batches for evaluation that outperform comparable methods for Bayesian optimization.
strength_weakness: strength_weakness: I think the method is interesting and will be of value to the field. However, I do not find that the experimental evaluation of the method provides full support for the claims of the paper.

**Issue 1: Evaluation limited to large-batch setting**

The major issue is that the paper claims the method is for general batch Bayesian optimization problems, without any qualifiers that I can see. The experiments all use q=100, which is a large batch size. Smaller batch sizes are often of interest too, e.g. batch sizes of 5, 10, and 50 in the GIBBON paper. The setting q=100 used here is also used in the TURBO paper (Eriksson et al. 2019) where it is described as a "large batch size." 

Given the experiment results in the paper, I don't know if this method will perform well for small- or medium-sized batches. Thus, either the experiments need to be expanded to include experiments with batch sizes such as 5 and 10, or the framing of the paper needs to be adjusted to emphasize that the method is specifically for large-batch Bayesian optimization, not general batch BO problems.

This issue also relates to the choice of baselines. The experiments only explore large-batch settings, where GIBBON (as the paper notes) is known to perform poorly. If the paper wants to claim that it performs better than GIBBON in general, then it needs to make that comparison on batch sizes of 5 and 10. If the paper wants to claim superiority only on large-batch settings, then that's fine to only use q=100, but then it needs to compare to state-of-the-art for large batch. Thompson sampling is a popular method for large-batch settings which is included as a baseline, but to my knowledge the state-of-the-art for large-batch BO is TURBO (Eriksson et al. 2019). In fact the q=100 and the robot pushing and rover trajectory problems are all exactly as in the TURBO paper, so its inclusion as a baseline is pretty obvious and, I think, necessary.

**Issue 2: Lack of statistical significance**

The results of the experiments do not appear to be statistically significant. The main results given in the main text are tables, and these tables do not have confidence intervals. The only place where uncertainty in the results is show are the figures in the appendix, and there the confidence intervals appear to overlap in most cases. This is due to the use of only 5 replicates. I appreciate that these experiments are costly to run since they are using 1000 iterations, nevertheless the lack of statistical significance in most of the results provides weak support for the claim that BEEBO is actually better, vs. what we're seeing just being noise in the experiments. The paper needs some form of statistical analysis to convince the reader that what we're seeing is not just noise in the experiments. The best way to do this would be to include confidence intervals in tables 2 and 3, and then increase the number of replicates as necessary to achieve statistical significance in the differences in means. I do not feel it appropriate to highlight the method as being "best value" when it is possible that the confidence interval for that value contains the values of the other methods.
questions: questions: The issues raised above can be addressed by running more experiments (batch sizes 5 and 10; TURBO; more replicates to get reasonable CIs for the results tables). But it will probably require more experiments that can be run in the rebuttal period. Do we expect the method to work well for Q=5 and Q=10?
limitations: limitations: N/A
paper_summary: paper_summary: Proposing a new acquisition function inspired by statistical physics, which allows explicit control of exploration-exploitation trade-offs in a batch BO setting.
strength_weakness: strength_weakness: **Major Points**
1. **Lack of Unique Selling Point:** 
   The method does not appear to solve any unique cases that other methods cannot. While the related work section outlines many similar approaches, this work only compares itself to q-UCB. Without a theoretical study, such as convergence rate analysis, there is insufficient motivation to adopt another heuristic approach like this. To demonstrate efficacy, a comprehensive empirical study with extensive comparisons to existing works is necessary, given the unclear advantages.

2. **Review of Claimed Selling Points:**
   - **Not Based on MC Integration Like BoTorch:** 
     While this is true, it is unclear if it is beneficial. MC approaches are approximations but have convergence guarantees (refer to the BoTorch paper's appendix). This work lacks such guarantees.
   - **Tight Control of Exploration-Exploitation Trade-off:** 
     The proposed method is not the only solution. UCB theory can bound the feasible region by [max LCB(x), UCB(x)] with 1 - $\delta$ probability (e.g., see [1]). This region can be controlled by $\beta$ hyperparameters, corresponding to $\delta$. Constrained batch BO within this region would yield similar results with theoretical guarantees.
   - **Heteroskedastic BO:** 
     There are no comparisons with existing methods. UCB variants can address this problem. Vanilla UCB theory does not differentiate between epistemic and aleatoric uncertainty. Therefore, UCB with heteroscedastic GP can serve the same purpose. For example, training a heteroscedastic GP with the observed dataset and replacing the inner GP on noise variance with normal isotropic likelihood variance (pure epistemic uncertainty) will yield similar results as risk-averse batch BO, without needing the change in modelling the acquisition functions regardless of hetero-/homo-scedastic cases like this work.

3. **Setting k in Practice:** 
   Setting  $k$ is not an easy task for users. In UCB, $\beta$ presents a similar challenge, but there are theoretical guidelines and algorithm to compute this (e.g., [2]).

4. **Constrained Uncertainty Sampling:** 
   This work can be understood as a variant of constrained uncertainty sampling. As [3] explains, variance reduction can be written similarly to the entropy proposed in this paper (see section 2.2). It also shows that the variance-only approach is inferior to UCB both theoretically and empirically. The batch setting may lead to model misspecification, particularly when hallucination (aka fantasization) is applied. The concerns and approach are notably similar to ([49] in your citation number), making a comparison with their method unavoidable.

5. **Data Augmentation Procedure:** 
   The explanation is unclear. Is it fantasization (aka hallucination) or simply observed points? How does this differ from Eq.(3) in [3]?

- [1] Fengxue Zhang, Jialin Song, James C Bowden, Alexan- der Ladd, Yisong Yue, Thomas Desautels, and Yuxin Chen. Learning regions of interest for bayesian optimiza- tion with adaptive level-set estimation. In International Conference on Machine Learning, pages 41579–41595. PMLR, 2023.
- [2] Kihyuk Hong, Yuhang Li, and Ambuj Tewari. An optimization-based algorithm for non-stationary kernel bandits without prior knowledge. In International Conference on Artificial Intelligence and Statistics, pages 3048–3085. PMLR, 2023.
- [3] Srinivas, N., Krause, A., Kakade, S., & Seeger, M. (2010). Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. In Proceedings of the 27th International Conference on Machine Learning (pp. 1015-1022).
questions: questions: Questions are written in the above weakness section.
limitations: limitations: Limitations are discussed in the discussion section.
paper_summary: paper_summary: This paper introduces a new acquisition function BEBBO for batched BO. BEBBO tries to build (negative) free-energies-like acquisition function, enabling gradient-based optimization, tight exploration-exploitation control, and risk-averse BO under heteroskedastic noise. It tries to improve existing parallel acquisition functions in the following ways:
1.	uses a hyper-parameter T to directly balance exploration and exploitation by separating these two parts clearly;
2.	keeps the behavior predictable by scaling E and I with batch size;
3.	enables the optimization of gradient descent by holding the availability of closed-form expressions for GP.
This paper demonstrates several experimental comparisons and shows its effectiveness.
strength_weakness: strength_weakness: 1.	The idea is straightforward: simply combine two common components—entropy reduction and the weighted sum of values. In my opinion, the novelty is not very strong. 

2.	The article doesn't discuss the situation when BEEBO is used with other surrogate models. It seems that if BEEBO is not used with GP, it loses the advantages of closed-form expressions and gradient descent optimization. 

3.	In control problems shown in Figure A2, the performances of meanBEBBO and maxBEBBO are not outstanding. Especially, these two variants are surpassed by KB in Robot pushing problem.

4.	Although BEEBO performs well on many synthetic test problems, its versatility and effectiveness require more experimental validation in specific applications.

5.	In the experiments in main text, the authors only showed the comparison with q-UCB. It would be better to show comparisons with other batched baselines and provide a thorough analysis. The comparison with q-UCB shows the advantage on the balance between exploration and exploitation. But other advantages emphasized in the paper, such as the benefits of gradient descent optimization and the tight control of the explore-exploit trade-off, are not fully demonstrated. I suggest that the authors can re-organize the paper to move the comparison with other baselines to the main paper. 

6.	The theoretical analysis is not deep. No regret bound is analyzed.
questions: questions: - How to set the parameter T properly? Do you have some instructions?

- In line 158, why does I(x) also scale linearly with batch size? 

- The Equation A5 is not shown.

- What are the advantages of multiplying batch size Q in E(X)?

- How does the behavior of BEBBO change with batch size?
limitations: limitations: Yes

Comments: comment: Thank you for the additional analyses. It is good to see that BEEBO continues to perform well in the small-to-medium batch setting, and that it performs well relative to TurBO. It is interesting that the sensitivity to T' appears to be higher for q=5 and q=10, a phenomenon that probably merits investigation.

As progress was made on my two main concerns, I will raise my score. Although other reviews have raised good points that I hope the authors will consider further.
I appreciate the author's responses to my concerns, especially regarding the interpretability and the related works. Concerning the other potentially related works, I believe my fellow reviewers have proposed sufficient concrete examples. My additional comments on the algorithm's possible advantage over existing theoretical sound algorithms align with the author's recent response. The applicability of the acquisition optimizer is a known problem. Therefore, it is still meaningful to introduce algorithms like BEEBO. I hope the author can address the remaining concerns of other reviewers. At the same time, I maintain my original evaluation.

***Reference:***
In section B by Shahriari et al. 2016: "Unfortunately, these auxiliary optimization techniques can be problematic for several reasons." "First, in practice, it is difficult to assess whether the auxiliary optimizer has found the global maximizer of the acquisition function."

- Shahriari, Bobak, et al. "Taking the human out of the loop: A review of Bayesian optimization." Proceedings of the IEEE 104.1 (2015): 148-175.
Thanks for the response. I keep my evaluation.
Thank you for taking the time to consider our rebuttal!  We appreciate the effort to work together to improve the clarity of our work.

**MC AF optimization**:  
We agree that MC is a powerful paradigm and that different strategies provide different benefits. Optimizing MC AFs uses sample average approximation to optimize towards the true maximizer $\alpha^* $ . As you also mentioned in your comments, Appendix D.3 in the BoTorch paper provides more theoretical background to that approach. As one would expect for MC approximations, there is a critical dependency on the number of samples: The results of theorems 1 and 2 for approximating $\alpha^{*}$ both hold in the **limit $i \to \infty$, an infinite number of samples**.
In practice, one will naturally face a situation with a rather limited number of samples, as also indicated in various BoTorch tutorials (e.g. 256 samples in the $q$-NEI tutorial), so the true AF is by definition *approximated*. 

As mentioned in D.3, this will therefore leave us with an **integration error**, that becomes more critical in large $q$ as the AF is of dimension $\mathbb{R}^{q}$. Our closed-form expression does by design not suffer from an integration error.

Note that BoTorch also uses multi-start L-BFGS-B for MC AFs to mitigate the fact that the optimization approximates the true AF (BoTorch paper section 4.1).

**Tight control**:  
Thank you for clarifying how the BALLET comment was intended. Again, our argument is not that UCB *theory* does not provide control (we are convinced it does), but that the practical limitations of BoTorch's MC methods at large $q$ outlined above do not deliver that control in practice, as we see empirically.

**Heteroscedastic BO**:  
We would really appreciate a reference, as it is not intuitive to us. Our understanding of heteroskedastic noise and UCB is rooted in the analysis of Makarova et al. 2021, which deals with heteroskedastic GPs, and shows how UCB in fact needs to be modified to deliver risk-averse behaviour.
Thank you for your reply, and I acknowledge your effort. Since I had not initially reviewed the appendix, my concern about weak evaluation was misplaced. However, I believe the main paper should be self-contained to convey key points effectively. While it is great to include many experimental results, they should be accompanied by clear explanations of "why". While theory is the easiest to answer, experimental investigation beyond final convergence comparison is also crucial.

I would like to clarify a few questions, as it seems the authors may not have fully read the related papers. Given the limited time, I understand that reviewing all related literature can be challenging. Therefore, I believe another round of conference review would be beneficial, allowing more time to refine the work. Thus, I still recommend rejection but have raised my score to 4 in recognition of the extensive experimental efforts.

**MC Acquisition Function (AF):** To decide the next query point, the AF must be maximised, which involves an "inner" optimisation loop that is non-convex. The MC AF leverages sub-modularity and has proven to converge to the true global maximum of the AF in each round. The closed-form AF requires heuristics, such as multi-start local optimer like L-BFGS-B. In this context, I do not believe a closed-form expression is necessarily better than MC methods, as they offer different benefits.

**Tight Control:** This serves as a counter-example that "UCB theory" can achieve tight control through $\beta$. BALLET is an example of a method that restricts exploration via confidence bounds. This is not about suggesting another baseline comparison; rather, it is about preferring simpler, well-established approaches. When a provably converging method easily comes up, why introduce a new heuristic?

**Heteroscedastic BO:** While I do not have a specific reference, this idea is fairly intuitive and the paper should exist somewhere. It serves as another counter-example.

Good luck with the revision!