Approximation-Aware Bayesian Optimization

Authors: natalie  maus
kyurae  kim
geoff  pleiss
david  eriksson
meta  john
p  cunningham
jacob r gardner
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 01:24:30

Abstract: High-dimensional Bayesian optimization (BO) tasks such as molecular design often require $>10,$$000$ function evaluations before obtaining meaningful results. While methods like sparse variational Gaussian processes (SVGPs) reduce computational requirements in these settings, the underlying approximations result in suboptimal data acquisitions that slow the progress of optimization. In this paper we modify SVGPs to better align with the goals of BO: targeting informed data acquisition over global posterior fidelity. Using the framework of utility-calibrated variational inference (Lacoste–Julien et al., 2011), we unify GP approximation and data acquisition into a joint optimization problem, thereby ensuring optimal decisions under a limited computational budget. Our approach can be used with any decision-theoretic acquisition function and is readily compatible with trust region methods like TuRBO (Eriksson et al., 2019). We derive efficient joint objectives for the expected improvement (EI) and knowledge gradient (KG) acquisition functions in both the standard and batch BO settings. On a variety of recent high dimensional benchmark tasks in control and molecular design, our approach significantly outperforms standard SVGPs and is capable of achieving comparable rewards with up to $10\times$ fewer function evaluations.

Main paper: 
 introduction: 
 Bayesian optimization (BO; Frazier, 2018; Garnett, 2023; Jones et al., 1998; Mockus, 1982; Shahriari et al., 2015) casts optimization as a sequential decision-making problem. Many recent successes of BO have involved complex and high-dimensional problems. In contrast to "classic" low-dimensional BO problems-where expensive black-box function evaluations far exceeded computational costs-these modern problems necessitate tens of thousands of function evaluations, and it is often the complexity and dimensionality of the search space that makes optimization challenging, rather than a limited evaluation budget (Eriksson et al., 2019; Griffiths and Hernández-Lobato, 2020; Maus et al., 2022 Maus et al., , 2023;; Stanton et al., 2022) . Because of these scenarios, BO is entering a regime where computational costs are becoming a primary bottleneck (Maddox et al., 2021; Maus et al., 2023; Moss et al., 2023; Vakili et al., 2021) , as the Gaussian process (GP; Rasmussen and Williams, 2005) surrogate models that underpin most of Bayesian optimization scale cubically with the number of observations. In this new regime, we require scalable GP approximations, an area that has made tremendous progress over the last decade. In particular, sparse variational Gaussian processes (SVGP; Hensman et al., 2013; Quiñonero-Candela and Rasmussen, 2005; Titsias, 2009) have seen an increase in use (Griffiths and Hernández-Lobato, 2020; Maddox et al., 2021; Maus et al., 2022 Maus et al., , 2023;; Stanton et al., 2022; Tripp et al., 2020; Vakili et al., 2021) , but many challenges remain to effectively deploy SVGPs for large-budget BO. In particular, the standard SVGP training objective is not aligned with the goals of black-box optimization. SVGPs construct an inducing point approximation that maximizes the standard variational evidence lower bound (ELBO; Jordan et al., 1999) , yielding a posterior approximation 𝑞 * (𝑓) that models all observed data (Matthews et al., 2016; Moss et al., 2023) . However, the optimal posterior approximation 𝑞 * is suboptimal for the decision-making tasks involved in BO (Lacoste-Julien et al., 2011) . In BO, we do not care about posterior fidelity at the majority of prior observations; rather, we only care about the fidelity of downstream functions involving the posterior, such as the expected utility. To illustrate this point intuitively, consider using the common expected improvement (EI; Jones et al., 1998) acquisition function for selecting new observations. Maximizing the ELBO might result in a posterior approximation that maintains fidelity for training examples in regions of virtually zero EI, thus wasting "approximation budget."To solve this problem, we focus on the deep connections between statistical decision theory (Robert, 2001; Wasserman, 2013, §12) and Bayesian optimization (Garnett, 2023, §6-7) , where acquisition maximization can be viewed as maximizing posterior-expected utility. Following this perspective, we leverage the utility-calibrated approximate inference framework (Jaiswal et al., 2020 (Jaiswal et al., , 2023;; Lacoste-Julien et al., 2011) , and solve the aforementioned problem through a variational bound (Blei et al., 2017; Jordan et al., 1999) -the (log) expected utility lower bound (EULBO)-a joint function of the decision (the BO query) and the posterior approximation (the SVGP). When optimized jointly, the EULBO automatically yields the approximately optimal decision through the minorizemaximize principle (Lange, 2016) . The EULBO is reminiscent of the standard variational ELBO (Jordan et al., 1999) , and can indeed be viewed as a standard ELBO for a generalized Bayesian inference problem (Bissiri et al., 2016; Knoblauch et al., 2022) , where we seek to approximate the utility-weighted posterior. This work represents the first application of utility-calibrated approximate inference towards BO despite its inherent connection with utility maximization.The benefits of our proposed approach are visualized in Fig. 1 . Furthermore, it can be applied to acquisition function that admits a decision-theoretic interpretation, which includes the popular expected improvement (EI; Jones et al., 1998) and knowledge gradient (KG; Wu et al., 2017) acquisition functions, and is trivially compatible with local optimization techniques like TuRBO (Eriksson et al., 2019) for high-dimensional problems. We demonstrate that our joint SVGP/acquisition optimization approach yields significant improvements across numerous Bayesian optimization benchmarks. As an added benefit, our approach can simplify the implementation and reduce the computational burden of complex (decision-theoretic) acquisition functions like KG. We demonstrate a novel algorithm derived from our joint optimization approach for computing and optimizing the KG that expands recent work on one-shot KG (Balandat et al., 2020) and variational GP posterior refinement (Maddox et al., 2021) .Overall, our contributions are summarized as follows:• We propose utility-calibrated variational inference of SVGPs in the context of large-budget BO.• We study this framework in two special cases using the utility functions of two common acquisition functions: EI and KG. For each, we derive tractable EULBO expressions that can be optimized.• For KG, we demonstrate that the computation of the EULBO takes only negligible additional work over computing the standard ELBO by leveraging an online variational update. Thus, as a byproduct of optimizing the EULBO, optimizing KG becomes comparable to the cost of the EI.• We extend this framework to be capable of running in batch mode, by introducing q-EULBO analogs of q-KG and q-EI as commonly used in practice (Wilson et al., 2018) .• We demonstrate the effectiveness of our proposed method against standard SVGPs trained with ELBO maximization on high-dimensional benchmark tasks in control and molecular design, where the dimensionality and evaluation budget go up to 256 and 80k, respectively. 
background: 
 Noisy Black-Box Optimization. Noisy black-box optimization refers to problems of the form: maximize 𝒙∈𝒳 𝐹 (𝒙) , where 𝒳 ⊂ ℝ 𝑑 is some compact domain, 𝐹 ∶ 𝒳 → 𝒴 is some objective function, and we assume that only zeroth-order information of 𝐹 is available. More formally, for some 𝑖 ∈ ℕ >0 , we assume that observations of the objective function (𝒙 𝑖 , 𝑦 𝑖 = F (𝒙 𝑖 )) have been corrupted by independently and identically distributed (i.i.d.) Gaussian noise F (𝒙 𝑖 ) ≜ 𝐹(𝒙 𝑖 ) + 𝜖, where 𝜖 ∼ 𝒩(0, 𝜎 2 n ). The noise variance 𝜎 2 n is also unknown. ❸ After selecting 𝒙 𝑡+1 , F is evaluated to obtain the new datapoint (𝒙 𝑡+1 , 𝑦 𝑡+1 = F (𝒙 𝑡+1 )). This is then added to the dataset, forming 𝒟 𝑡+1 = 𝒟 𝑡 ∪ (𝒙 𝑡+1 , 𝑦 𝑡+1 ) to be used in the next iteration.Utility-Based Acquisition Functions. Many commonly used acquisition functions, including EI and KG, can be expressed as posterior-expected utility functionsEQUATIONwhere 𝑢 (𝒙, 𝑓; 𝒟) ∶ 𝒳 × ℱ → ℝ is some utility function associated with 𝛼 (Garnett, 2023, §6-7) . In statistical decision theory, posterior-expected utility maximization policies such as 𝛿 𝛼 are known as Bayes policies. These are important because, for a given utility function, they attain certain notions of statistical optimality such as Bayes optimality and admissibility (Robert, 2001, §2.4; Wasserman, 2013, §12) . However, this only holds true if we can exactly compute Eq. ( 2) over the posterior. Once approximate inference is involved, making optimal Bayes decisions becomes challenging.Sparse Variational Gaussian Processes. While the 𝒪(𝑛 3 ) complexity of exact Gaussian process model selection and inference is not necessarily a roadblock in the traditional regression setting with 10,000-50,000 training examples, BO amplifies the scalability challenge by requiring us to sequentially train or update many large scale GPs as we iteratively acquire more data.To address this, sparse variational GPs (SVGP; Hensman et al., 2013; Titsias, 2009) have become commonly used in high-throughput Bayesian optimization. SVGPs modify the original GP prior from 𝑝(𝑓) to 𝑝(𝑓 | 𝒖)𝑝(𝒖), where we assume the latent function 𝑓 is "induced" by a finite set of inducing values 𝒖 = (𝑢 1 , … , 𝑢 𝑚 ) ∈ ℝ 𝑚 located at inducing points 𝒛 𝑖 ∈ 𝒳 for 𝑖 = 1, … , 𝑚. Inference is done through variational inference (Blei et al., 2017; Jordan et al., 1999) , where the posterior of the inducing points is approximated using 𝑞 𝝀 (𝒖) = 𝒩 (𝒖; 𝝀 = (𝒎, 𝑺)) and that of the latent functions with 𝑞 (𝑓 | 𝒖) = 𝑝 (𝑓 | 𝒖). Here, the variational parameters 𝒎 and 𝑺 are defined as the learned mean and covariance of the variational distribution 𝑞 𝝀 (𝒖). It is standard practice to define 𝝀 = (𝒎, 𝑺) so that 𝝀 can be used as shorthand to represent all of the trainable variational parameters. As is typical in the BO literature, we use the subscript 𝝀 ∈ Λ to denote that the distribution denoted as 𝑞 contains trainable parameters in 𝝀.For a positive definite kernel function 𝑘 ∶ 𝒳 × 𝒳 → ℝ >0 , the resulting ELBO objective, which can be computed in a closed form (Hensman et al., 2013) with k𝒙𝒙 ≜ 𝑘 (𝒙, 𝒙) -𝒌 𝒙𝒁 𝑲 -1 𝒁𝒁 𝒌 ⊤ 𝒁𝒙 , the vector 𝒌 𝒁𝒙 ∈ ℝ 𝑚 is formed as [𝒌 𝒁𝒙 ] 𝑖 = 𝑘 (𝒛 𝑖 , 𝒙), and the matrix 𝑲 𝒁𝒁 ∈ ℝ 𝑚×𝑚 is formed as [𝑲 𝒁𝒁 ] 𝑖𝑗 = 𝑘(𝒛 𝑖 , 𝒛 𝑗 ). Additionally, the GP likelihood and kernel contain hyperparameters, which we denote as 𝜽 ∈ Θ, and we collectively denote the set of inducing point locations as 𝒁 = (𝒛 1 , … , 𝒛 𝑚 ) ∈ 𝒳 𝑚 . We therefore denote the ELBO as ℒ ELBO (𝝀, 𝒁, 𝜽; 𝒟 𝑡 ). 
approximation-aware bayesian optimization: 
 When SVGPs are used in conjunction with BO (Maddox et al., 2021; Moss et al., 2023) at iteration 𝑡 ≥ 0, acquisition functions of the form of Eq. ( 2) are naïvely approximated as 𝛼 (𝒙; 𝒟) ≈ ∫ 𝑢 (𝒙, 𝑓; 𝒟 𝑡 ) 𝑞 𝝀 (𝑓) d𝑓, where 𝑞 𝝀 (𝑓) is the approximate SVGP posterior given by Eq. ( 4). The acquisition policy implied by this approximation contains two separate optimization problems: (Matthews et al., 2016; Moss et al., 2023) , without regard for how the resulting model performs at selecting the next function evaluation in the BO loop. For an illustration of this, see Figure 1 . Instead, we propose a modification to SVGPs that couples the posterior approximation and data acquisition through a joint problem of the form:𝒙 𝑡+1 =EQUATIONThis results in 𝒙 𝑡+1 directly approximating a solution to Eq. ( 2), where the expected utility lowerbound (EULBO) is an ELBO-like objective function derived below.Consider an acquisition function of the form of Eq. ( 2), where the utility 𝑢 ∶ 𝒳 × ℱ → ℝ >0 is strictly positive. We can derive a similar variational formulation of the acquisition function maximization problem following Lacoste-Julien et al. (2011) where 𝑍 is a normalizing constant. A restriction on 𝑢 comes from the inequality in Eq. ( 7), where the utility needs to be strictly positive. This means that non-strictly positive utilities need to be modified to be incorporated into this framework. (See the examples by Kuśmierczyk et al., 2019. ) Also, notice that the derivation is reminiscent of expectation-maximization (Dempster et al., 1977) and variational lower bounds (Jordan et al., 1999) . That is, through the minorize-maximize principle (Lange, 2016) , maximizing the lower bound with respect to 𝒙 and 𝝀 approximately solves the original problem of maximizing the posterior-expected utility.Expected Utility Lower-Bound. Up to a constant and rearranging terms, maximizing Eq. ( 7) is equivalent to maximizing which is the joint objective function alluded to in Eq. ( 6). We maximize EULBO to obtain (𝒙 𝑡+1 , 𝝀 * ) = arg max 𝒙∈𝒳,𝝀∈Λ ℒ EULBO (𝒙, 𝝀), where 𝒙 𝑡+1 corresponds our next BO "query".From Eq. ( 8), the connection between the EULBO and ELBO is obvious: the EULBO is now "nudging" the ELBO solution toward high utility regions. An alternative perspective is that we are approximating a generalized posterior weighted by the utility (Table . 1 by Knoblauch et al., 2022; Bissiri et al., 2016) . Furthermore, Jaiswal et al. (2020 Jaiswal et al. ( , 2023) ) prove that the resulting actions satisfy consistency guarantees under assumptions typical in such results for variational inference (Wang and Blei, 2019) . We emphasize here that the SVGP-associated parameters 𝝀, 𝜽, 𝒁 have gradients that are determined by both terms above. Thus, the expected log-utility term 𝔼 𝑓∼𝑞 𝝀 (𝑓) log 𝑢 (𝒙, 𝑓; 𝒟 𝑡 ) simultaneously results in acquisition of 𝒙 𝑡+1 and directly influences the underlying SVGP regression model.The EI acquisition function can be expressed as a posterior-expected utility, where the underlying "improvement" utility function is given by the difference between the objective value of the query, 𝑓 where ReLU (𝑥) ≜ max (𝑥, 0). Unfortunately, this utility is not strictly positive whenever 𝑓 (𝒙) ≤ 𝑦 * . Thus, we cannot immediately plug 𝑢 EI into the EULBO. While it is possible to add a small positive constant to 𝑢 EI and make it strictly positive as done by Kuśmierczyk et al. (2019) , this results in a looser Jensen gap in Eq. ( 7), which could be detrimental. This also introduces the need for tuning the constant, which is not straightforward. Instead, we define the following "soft" EI utility:𝑢 SEI (𝒙, 𝑓; 𝒟 𝑡 ) ≜ sof tplus ( 𝑓 (𝒙) -𝑦 * 𝑡 ) ,where the ReLU in Eq. ( 9) is replaced with sof tplus (𝑥) ≜ log (1 + exp(𝑥)). sof tplus(𝑥) converges to the ReLU in both extremes of 𝑥 → ±∞. Thus, 𝑢 SEI will behave closely to 𝑢 EI , while being slightly more explorative due to positivity.Computing the EULBO and its derivatives now requires the computation of 𝔼 𝑓∼𝑞 𝝀 (𝑓) log 𝑢 SEI (𝒙, 𝑓; 𝒟 𝑡 ), which, unlike EI, does not have a closed-form. However, since the utility function only depends on the function values of 𝑓, the expectation can be efficiently computed to high precision through one-dimensional Gauss-Hermite quadrature. Crucially, the expensive 𝐾 -1 𝑧𝑧 𝑚 and 𝐾 -1 𝑧𝑧 𝑆𝐾 -1 𝑧𝑧 solves that dominate both the asymptotic and practical running time of both the ELBO and the EULBO are fixed across the log utility evaluations needed by quadrature. Because quadrature only depends on these precomputed moments, the additional work necessary due to lacking a closed form solution is negligible: Gauss-Hermite quadrature converges extremely quickly in the number of quadrature sites, and only requires on the order of 10 or so of these post-solve evaluations to achieve near machine precision.Although non-trivial, the KG acquisition is also a posterior-expected utility, where the underlying utility function is given by the maximum predictive mean value anywhere in the input domain after conditioning on a new observation (𝒙, 𝑦) ∈ 𝒳 × 𝒴: 𝑢 KG (𝒙, 𝑦; 𝒟 𝑡 ) ≜ max KG; Frazier, 2009; Garnett, 2023) Note that the utility function as defined above is not non-negative: the maximum predictive mean of a Gaussian process can be negative. For this reason, the utility function is commonly (and originally, e.g. Frazier, 2009, Eq. 4.11) written in the literature as the difference between the new maximum mean after conditioning on (𝒙, 𝑦) and the maximum mean beforehand:EQUATION𝑢 KG (𝒙, 𝑦; 𝒟 𝑡 ) ≜ max 𝒙 ′ ∈𝒳 𝔼 [ 𝑓(𝒙 ′ ) | 𝒟 𝑡 ∪ {(𝒙, 𝑦)} ] -𝜇 + 𝑡 ,where𝜇 + 𝑡 ≜ max 𝒙 ′′ ∈𝒳 𝐸 [ 𝑓(𝒙 ′′ ) | 𝒟 𝑡 ]. Note that 𝜇 + 𝑡 plays the role of a simple constant as it depends on neither 𝒙 nor 𝑦. Similarly to the EI acquisition, this utility is still not strictly positive, and we thus define its "softplus-ed" variant: 𝑢 SKG (𝒙, 𝑦; 𝒟 𝑡 ) ≜ sof tplus ( 𝑢 KG (𝒙, 𝑦; 𝒟 𝑡 ) -𝑐 + ) .Here, 𝑐 + acts as 𝜇 + 𝑡 by making 𝑢 KG positive as often as possible. This is particularly important when the GP predictive mean is negative as a consequence of the objective values being negative. One natural choice of constant is using 𝜇 + 𝑡 ; however, we find that simply choosing 𝑐 + = 𝑦 + 𝑡 works well and is more computationally efficient. Here, 𝑦 + 𝑡 is the highest value of 𝑦 𝑡 (the highest objective value observed so far).One-Shot KG EULBO. The EULBO using 𝑢 SKG results in an expensive nested optimization problem.To address this, we use an approach similar to the one-shot knowledge gradient method of Balandat et al. (2020) ] 𝑆-times. Notice that we only need to compute an updated posterior predictive mean, and can ignore predictive variances. For this, we can leverage the online updating strategy of Maddox et al. (2021) . In particular, the predictive mean can be updated in 𝒪(𝑚 2 ) time using a simple Cholesky update. The additional 𝒪(𝑆𝑚 2 ) cost of computing the EULBO is therefore amortized by the original 𝒪(𝑚 3 ) cost of computing the ELBO.The EULBO can be extended to support batch Bayesian optimization by using the Monte Carlo batch mode analogs of utility functions as discussed e.g. by Balandat et al. (2020) ; Wilson et al. (2018) ) .As done in the BoTorch software package (Balandat et al., 2020) , we observe that fixing the set of base samples 𝜖 1 , ..., 𝜖 𝑆 during each BO iteration results in better optimization performance at the cost of negligible q-EULBO bias. Now, optimizing the q-EULBO is done over the full set of 𝑞 candidates ( 𝒙 1 , ..., 𝒙 𝑞 ) jointly, as well as the GP hyperparameters, inducing points, and variational parameters. resulting in a similar analog to q-KG as described by Balandat et al. (2020) .Optimizing the ELBO for SVGPs is known to be challenging (Galy-Fajou and Opper, 2021; Terenin et al., 2024) as the optimization landscape for the inducing points is non-convex, multi-modal, and non-smooth. Naturally, these are also challenges for EULBO; we found that care must be taken when implementing and initializing the EULBO maximization problem. In this subsection, we outline some key ideas, while a detailed description with pseudocode is presented in Appendix A.Initialization and Warm-Starting. We warm-start the EULBO maximization procedure by solving the conventional two-step scheme in Eq. ( 5): At each BO iteration, we obtain the "warm" initial values for (𝝀, 𝒁, 𝜽) by optimizing the standard ELBO. Then, we use this to maximize the conventional acquisition function corresponding to the chosen utility function 𝑢 (the expectation of 𝑢 over 𝑞 𝝀 (𝑓)), which provides the warm-start initialization for 𝒙.Alternating Maximization Scheme. To optimize ℒ EULBO (𝒙, 𝝀, 𝒁, 𝜽), we alternate between optimizing over the query 𝒙 and the SVGP parameters 𝝀, 𝒁, 𝜽. We find this block-coordinate descent scheme to be more stable and robust than jointly updating all parameters, though the reason why this is more stable than jointly optimizing all parameters requires further investigation. 
experiments: 
 We evaluate EULBO-based SVGPs on a number of benchmark BO tasks, described in detail in Section 4.1. These tasks include standard low-dimensional BO problems, e.g., the 6D Hartmann function, as well as 7 high-dimensional and high-throughput optimization tasks.Baselines. We compare EULBO to several baselines with the main goal of achieving a high reward using as few function evaluations as possible. Our primary point of comparison is ELBO-based SVGPs. We consider two approaches for inducing point locations: 1. optimizing inducing point locations via the ELBO (denoted as ELBO), 2. placing the inducing points using the strategy proposed by Moss et al. (2023) at each stage of ELBO optimization (denoted as Moss et al.). The latter offers improved BO performance over standard ELBO-SVGP in BO settings, yet-unlike our method-it exclusively targets inducing point placement and does not affect variational parameters or hyperparameters of the model. In addition, we compare to BO using exact GPs using 2, 000 function evaluations as the use of exact GP is intractable beyond this point due to the need to repeatedly fit models.Acquisition Functions and BO algorithms. For EULBO, we test the versions based on both the Expected Improvement (EI) and Knowledge Gradient (KG) acquisition functions as well as the batch variant. We test the baseline methods using EI only. On high-dimensional tasks (tasks with dimensionality above 10), we run EULBO and baseline methods with standard BO and with trust region Bayesian optimization (TuRBO) (Eriksson et al., 2019) . For the largest tasks (Lasso, Molecules) we use acquisition batch size of 20 (𝑞 = 20), and batch size 1 (𝑞 = 1) for all others.Implementation Details and Hyperparameters. Code to reproduce all results in the paper is available at https://github.com/nataliemaus/aabo. We implement EULBO and baseline methods using the GPyTorch (Gardner et al., 2018) and BoTorch (Balandat et al., 2020) packages. For all methods, we initialize using a set of 100 data points sampled uniformly at random in the search space. We use the same trust region hyperparameters as in (Eriksson et al., 2019) . In Appendix B.1, we also evaluate an additional initialization strategy for the molecular design tasks. This alternative initialization matches prior work in using 10, 000 molecules from the GuacaMol dataset Brown et al. (2019) rather than the details we used above for consistency across tasks, but does achieve higher overall performance.Hartmann 6D. The widely used Hartmann benchmark function (Surjanovic and Bingham, 2013) .Lunar Lander. The goal of this task is to find an optimal 12-dimensional control policy that allows an autonomous lunar lander to consistently land without crashing. The final objective value we optimize is the reward obtained by the policy averaged over a set of 50 random landing terrains. For this task, we use the same controller setup used by Eriksson et al. (2019) .Rover. The rover trajectory optimization task introduced by Wang et al. ( 2018) consists of finding a 60-dimensional policy that allows a rover to move along some trajectory while avoiding a set of obstacles. We use the same obstacle set up as in Maus et al. (2023) .Figure 3 : Ablation study measuring the impact of EULBO optimization on various SVGP parameters. At each BO iteration, we use the standard ELBO objective to optimize the SVGP hyperparameters, variational parameters, and inducing point locations. We then refine some subset of these parameters by further optimizing them with respect to the EULBO objective.Lasso DNA. We optimize the 180-dimensional DNA task from the LassoBench library (Šehić et al., 2022) of benchmarks based on weighted LASSO regression (Gasso et al., 2009) .We select four challenging tasks from the Guacamol benchmark suite of molecular design tasks (Brown et al., 2019) In Figure 2 , we plot the reward of the best point found by the optimizer after a given number of function evaluations. Error bars show the standard error of the mean over 20 replicate runs. EULBO with TuRBO outperforms the other baselines with TuRBO. Similarly, EULBO with standard BO outperforms the other standard BO baselines. One noteworthy observation is that neither acquisition function appears to consistently outperform the other. However, EULBO-SVGP almost always dominates ELBO-SVGP and often requires a small fraction of the number of oracle calls to achieve comparable performance. These results suggest that coupling data acquisition with approximate inference/model selection results in significantly more sample-efficient optimization.While the results in Fig. 2 demonstrate that EULBO-SVGP improves the BO performance it is not immediately clear to what extent joint optimization modifies the posterior approximation beyond what is obtained by standard ELBO optimization. To that end, in Fig. 3 we refine an ELBO-SVGP model with varying degrees of additional EULBO optimization. At every BO iteration we begin by obtaining a SVGP model (where the variational parameters, inducing point locations, and GP hyperparameters are all obtained by optimizing the standard ELBO objective). We then refine some subset of parameters (either the inducing points, the variational parameters, the GP hyperparameters, or all of the above) through additional optimization with respect to the EULBO objective. Interestingly, we find that tasks respond differently to the varying levels of EULBO refinement. In the case of Lasso DNA, there is not much of a difference between EULBO refinement on all parameters versus refinement on the variational parameters alone. On the other hand, the performance on Median Molecules 2 is clearly dominated by refinement on all parameters. Nevertheless, we see that EULBO is always beneficial, whether applied to all parameters or some subset. 
related work: 
 Scaling Bayesian Optimization to the Large-Budget Regime. BO has traditionally been confined to the small-budget optimization regime with a few hundred objective evaluations at most. However, recent interest in high-dimensional optimization problems has demonstrated the need to scale BO to large data acquisition budgets. For problems with ∼10 2018) turn to ensembles of subsampled GPs. For problems with ≫ 10 3 acquisitions, SVGP has become the de facto approach to alleviate computational complexity (Griffiths and Hernández-Lobato, 2020; Maus et al., 2022 Maus et al., , 2023;; Stanton et al., 2022; Tripp et al., 2020; Vakili et al., 2021) . As in this paper, many works have proposed modifications to SVGP to improve its performance in BO applications. Moss et al. (2023) proposed an inducing point placement based on a heuristic modification of determinantal point processes (Kulesza and Taskar, 2012) , which we used for initialization, while Maddox et al.(2021) proposed a method for a fast online update strategy for SVGPs, which we utilize for the KG acquisition strategy.Utility-Calibrated Approximate Inference. The utility-calibrated VI objective was first proposed by Lacoste-Julien et al. (2011) , where they used a coordinate ascent algorithm to maximize it. Since then, various extensions have been proposed: Kuśmierczyk et al. (2019)  
limitations and discussion: 
 The main limitation of our proposed approach is increased computational cost. While EULBO-SVGP still retains the 𝑂(𝑚 3 ) computational complexity of standard SVGP, our practical implementation requires a warm-start: first fitting SVGP with the ELBO loss and then maximizing the acquisition function before jointly optimizing with the EULBO loss. Furthermore, EULBO optimization currently requires multiple tricks such as clipping and block-coordinate updates. In future work, we aim to develop a better understanding of the EULBO geometry in order to develop developing more stable, efficient, and easy-to-use EULBO optimization schemes. Nevertheless, our results in Section 4 demonstrate that the additional computation of EULBO yields substantial improvements in BO dataefficiency, a desirable trade-off in many applications. Moreover, EULBO-SVGP is modular, and our experiments capture a fraction of its potential use. It can be applied to any decision-theoretic acquisition function, and it is likely compatible with non-standard Bayesian optimization problems such as cost-constrained BO (Snoek et al., 2012) , causal BO (Aglietti et al., 2020) , and many more.More importantly, our paper highlights a new avenue for research in BO, where surrogate modeling, approximate inference, and data selection are jointly determined from a unified objective. Extending this idea to GP approximations beyond SVGP and acquisition functions beyond EI/KG may yield further improvements, especially in the increasingly popular high-throughput BO setting.We will now provide additional details on the implementation. For the implementation, we treat the SVGP parameters, such as the variational parameters 𝝀, inducing point locations 𝒁, and hyperparameters 𝜽, equally. Therefore, for clarity, we will collectively denote them as 𝒘 = (𝝀, 𝒁, 𝜽) such that 𝒘 ∈ 𝒲 ≜ Λ × 𝒳 𝑚 × Θ, and the resulting SVGP variational approximation as 𝑞 𝒘 . Then, the ELBO and EULBO are equivalently denoted as follows:ℒ ELBO (𝒘; 𝒟) ≜ ℒ ELBO (𝝀, 𝒁, 𝜽; 𝒟) ℒ EULBO (𝒙, 𝒘; 𝒟 𝒙 , 𝒟 𝒘 ) ≜ 𝔼 𝑓∼𝑞 𝒘 (𝑓) log 𝑢 (𝒙, 𝑓; 𝒟 𝒙 ) + ℒ ELBO (𝒘; 𝒟 𝒘 ) .Also, notice that the ℒ EULBO separately denote the dataset to be passed to the utility and the ELBO.(Setting 𝒟 𝑡 = 𝒟 𝒘 = 𝒟 𝒙 retrieves the original formulation in Eq. ( 8).)Alternating Updates We perform block-coordinate ascent on the EULBO by alternating between maximizing over 𝒙 as 𝒘. Using vanilla gradient descent, the 𝒙-update is equivalent to𝒙 ← 𝒙 + 𝛾 𝒙 ∇ 𝒙 ℒ EULBO (𝒙, 𝒘; 𝒟) = 𝒙 + 𝛾 𝒙 ∇ 𝒙 𝔼 𝑓∼𝑞 𝒘 (𝑓) log 𝑢 (𝒙, 𝑓; 𝒟) ,where 𝛾 𝒙 is the stepsize. On the other hand, for the 𝒘-update, we subsample the data such that we optimize the ELBO over a minibatch 𝑆 ⊂ 𝒟 of size 𝐵 = |𝑆| as𝒘 ← 𝒘 + 𝛾 𝒘 ∇ 𝒘 ℒ EULBO (𝒙, 𝒘; 𝑆, 𝒟) = 𝒘 + 𝛾 𝒘 ∇ 𝒘 ( 𝔼 𝑓∼𝑞 𝒘 (𝑓) log 𝑢 (𝒙, 𝑓; 𝒟) + ℒ ELBO (𝒘; 𝑆) ) ,where 𝛾 𝒘 is the stepsize. Naturally, the 𝒘-update is stochastic due to minibatching, while the 𝒙-update is deterministic. In practice, we leverage the Adam update rule (Kingma and Ba, 2015) instead of simple gradient descent. Together with gradient clipping, this alternating update scheme is much more robust than jointly updating (𝒙, 𝒘).Input: SVGP parameters 𝒘 0 = (𝝀 0 , 𝒁 0 , 𝜽 0 ), Dataset 𝒟 𝑡 , BO utility function 𝑢, Output: BO query 𝒙 𝑡+1 Overview of Pseudocode. The complete high-level view of the algorithm is presented in Algorithm 1, except for the acquisition-specific details. AdamStep 𝛾 (𝒙, 𝒈) applies the Adam stepsize rule (Kingma and Ba, 2015) to the current location 𝒙 with the gradient estimate 𝒈 and the stepsize 𝛾. In practice, Adam is a "stateful" optimizer, which maintains two scalar-valued states for each scalar parameter. For this, we re-initialize the Adam states at the beginning of each BO step.Initialization. In the initial BO step 𝑡 = 0, we initialize 𝒁 0 with the DPP-based inducing point selection strategy of Moss et al. (2023) . For the remaining SVGP parameters 𝝀 0 and 𝜽 0 , we used the default initialization of GPyTorch (Gardner et al., 2018) . For the remaining BO steps 𝑡 > 0, we use 𝒘 from the previous BO step as the initialization 𝒘 0 of the current BO step.Warm-Starting. Due to the non-convexity and multi-modality of both the ELBO and the acquisition function, it is critical to appropriately initialize the EULBO maximization procedure. As mentioned in Section 3.5, to warm-start the EULBO maximization procedure, we use the conventional 2-step scheme Eq. ( 5), where we maximize the ELBO and then maximize the acquisition function. For ELBO maximization, we apply Adam (Kingma and Ba, 2015) with the stepsize set as 𝛾 𝑤 until the convergence criteria (described below) are met. For acquisition function maximization, we invoke the highly optimized BoTorch.optimize.optimize_acqf function (Balandat et al., 2020) .As commonly done, we use the reshuffling subsampling strategy where the dataset 𝒟 𝑡 is shuffled and partitioned into minibatches of size 𝐵. The number of minibatches constitutes an "epoch." The dataset is reshuffled/repartitioned after going through a full epoch.For both maximizing the ELBO during warm-starting and maximizing the EULBO, we continue optimization until we stop making progress or exceed 𝑘 epochs number of epochs. That is if the ELBO/EULBO function value fails to make progress for 𝑛 f ail number of steps. Hyperparameters. The hyperparameters used in our experiments are organized in Table 1 . For the full-extent of the implementation details and experimental configuration, please refer to the supplementary code.We provide additional results and plots that were omitted from the main text.In Fig. 4 , we provide plots on additional results that are similar to those in Fig. 2 . On three of the molecule tasks, we use 10,000 random molecules from the GuacaMol dataset as initialization. This is more consistent with what has been done in previous works and achieves better overall optimization performance.Figure 4 : Additional optimization results on three molecule tasks using 10,000 random molecules from the GuacaMol dataset as initialization. Each line/shaded region represents the mean/standard error over 20 runs. We count oracle calls starting after these initialization evaluations for all methods.In this section, we provide additional plots separating out BO and TuRBO results to make visualization easier.Figure 5 : BO-only optimization results of Fig. 2 . We compare EULBO-SVGP, ELBO-SVGP, ELBO-SVGP with DPP inducing point placement (Moss et al., 2023) , and exact GPs. These are a subset of the same results shown in Fig. 2 . Each line/shaded region represents the mean/standard error over 20 runs.Figure 6 : TuRBO-only optimization results of Fig. 2 . We compare EULBO-SVGP, ELBO-SVGP, ELBO-SVGP with DPP inducing point placement (Moss et al., 2023) , and exact GPs. These are a subset of the same results shown in Fig. 2 . Each line/shaded region represents the mean/standard error over 20 runs.For the results with approximate-GPs in Section 4, we used 𝑚 = 100 inducing points. In Fig. 7 , we evaluate the effect of using a larger number of inducing points (𝑚 = 1024) for EULBO-SVGP and ELBO-SVGP. As in Fig. 2 , we compare running TuRBO with EULBO-SVGP and with ELBO-SVGP using 𝑚 = 100 inducing points used for both methods. We add two additional curves for TuRBO with EULBO-SVGP and TuRBO with ELBO-SVGP using 𝑚 = 1024 inducing points. Each line/shaded region represents the mean/standard error over 20 runs.Fig. 7 shows that the number of inducing points has limited impact on the overall performance of TuRBO, and EULBO-SVGP outperforms ELBO-SVGP regardless of the number of inducing points used.The results in Section 4 used a standard SVGP objective. In this section, we evaluate the effect of using an alternative objective: the parametric Gaussian process regressor (PPGPR; Jankowiak et al., 2020) objective. PPGPR differs from the standard SVGP objective in that the variational approximation is optimized to maximize the predictive accuracy instead of matching the posterior. We compare the choice of objective (PPGPR vs SVGP) in Fig. 8 and observe that the objective has limited impact on the overall performance of TuRBO. In particular, EULBO-EI outperforms ELBO-EI regardless of the GP objective. Completing all of the runs needed to produce all of the results in this paper therefore required roughly 2680 total GPU hours.Compute Resources Used During Preliminary Investigations. In addition to the computational resources required to produce experimental results in the paper discussed above, we spent approximately 500 hours of GPU time on preliminary investigations. This was done on the aforementioned internal cluster shown in Table 2 .In Table 3 , we provide average wall-clock run times of different methods on the Lasso DNA optimization task. 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: [No] . Justification: The paper is methodological, where the considered algorithm does not immediately pose societal risks.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] .Justification: The paper does not introduce new assets.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] .Justification: The paper does not involve human participants.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] .Justification: The paper does not involve live participants.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix:  


Decision: This work scales up Bayesian Optimisation to large datasets. The paper is clear and claims supported. All reviewers were enthusiastic about this work and all concerns raised were adequately addressed during the rebuttal. This is a solid piece of work that is relevant to the community.

Reviews: review: paper_summary: paper_summary: The paper proposes a new approach for scaling Bayesian Optimisation to large datasets. Contrary to previous approaches, which fitted a sparse GP and optimised acquisition function independently, the paper proposes a method, which jointly optimises the variational parameters of sparse GP and searches for the next point to query. This is done by adding a utility term to ELBO. Authors discuss how this objective can be connected to generalised Bayesian Inference. Authors evaluate the proposed method on a number of benchmark with high number of datapoints and show it outperforms other baselines.
strength_weakness: strength_weakness: - As mentioned by the authors, optimisation of the EULBO objective is a bit cumbersome and seems this process is also slower than the optimisation of standard ELBO. However, it might be worth to sacrifice a bit of computational efficiency for an increase in the sample complexity.

- While authors cite previous work, which proved that the selected action satisfies convergence guarantee, it does not directly prove anything about the performance of optimisation process choosing action in such a way (e.g. regret bound or convergence to optimum). The paper could be made much stronger if authors managed to connect those notions in some way. However, I also do not think it is critical, as the main contribution of the paper is empirical (although it would be greatly appreciated).
questions: questions: - Can the authors provide the running time of the proposed method (and baselines) ? It is ok if the proposed method is slower, but it would be good to quantify exactly how much slower it is.

- Do authors have any expectations on how would the method perform in comparison to an exact GP on a smaller dataset? While I understand that the method is particularly designed for large datasets, it would be interesting to see how the utility VI-based acquisition strategy compares to standard acquisition functions like EI. It would be nice to have at least such experiment in the camera-ready version to judge whether the improvement delivered by the method comes from better acquisition strategy or better modelling
limitations: limitations: Authors admit that the increased difficulty of optimising EULBO is a limitation of their work.
paper_summary: paper_summary: This paper proposes a modification to sparse variational Gaussian processes (SVGPs) used in Bayesian optimization (BO) to better align the SVGP posterior approximation with the goal of optimizing an acquisition function. The key idea is to jointly optimize the SVGP and the acquisition function using a unified objective called the expected utility lower bound (EULBO). This approach ensures the posterior approximation is well-suited for the downstream decision-making task. The authors derive efficient EULBO objectives for the expected improvement (EI) and knowledge gradient (KG) acquisition functions and demonstrate improved performance over standard SVGPs on several high-dimensional optimization benchmarks.
strength_weakness: strength_weakness: The paper focuses specifically on SVGPs, but it would be interesting to explore whether the proposed approach can be extended to other sparse GP approximations, even those without a tractable ELBO. This would broaden the applicability of the method and strengthen the contribution.

In Figure 2, the exact GP performs worse than the proposed method and some other baselines in the BO setting. This is counterintuitive, as one would expect the exact GP to be the gold standard when computationally feasible. The authors should provide more discussion and insights into this unexpected result.

The paper presents several variations of the proposed method (e.g., EI-EULBO, KG-EULBO, batch versions), but it lacks a clear recommendation on which variant to use in practice. A more in-depth discussion of the trade-offs between these variants and guidance on when to use each one would enhance the practical value of the work.
questions: questions: Can the proposed EULBO approach be extended to other sparse GP approximations beyond SVGPs, even if they do not have a tractable ELBO? If so, what challenges would need to be addressed, and how might the optimization problem be formulated?

Why does the exact GP perform worse than the proposed method and some other baselines in the BO experiments (Figure 2)? Is this due to the limitations of the exact GP in high dimensions, or are there other factors at play?

What are the key factors to consider when choosing between the different variants of the proposed method (e.g., EI-EULBO, KG-EULBO, batch versions)? In what scenarios would one variant be preferred over another?
limitations: limitations: The authors have addressed the limitations of their work in Section 6, acknowledging the increased computational cost of the EULBO approach compared to standard SVGPs. They also mention the need for multiple optimization tricks and the potential instability of the EULBO optimization problem. These are important limitations that users should be aware of when considering this method.

However, the authors do not discuss any potential negative societal impacts of their work. While the proposed method is primarily methodological and does not pose immediate societal risks, it would be valuable for the authors to briefly comment on any broader implications of improved BO in high-dimensional spaces, such as the potential for misuse or unintended consequences in certain application domains.
paper_summary: paper_summary: The paper proposes a new method for training sparse Gaussian processes (GPs) for large-budget Bayesian Optimization (BO), designed to facilitate the sequential decision tasks inherent in BO. This is achieved by introducing the “expected utility lower bound (EULBO),” which formulates sparse GP training as a joint optimization problem of BO queries and posterior approximation. The proposed framework is compatible with other significant advancements in the field, such as Expected Improvement (EI) and Knowledge Refinement (KR) acquisition functions, TuRBO, and batch optimization. This compatibility was demonstrated across several BO benchmarks.
strength_weakness: strength_weakness: While the paper is generally well-written, many derivations were abbreviated and some of the notation was not properly explained, making the paper more difficult to read than necessary. These issues, however, are easily amendable. I provide some suggestions for improvement in the following subsections.
questions: questions: Introduction:  

(line 51): '\emph{minorize-maximize}'

Section 2:

(Fig. 1): It would be beneficial to show the SVGP model fit for EULBO as well.

(line 110): Could you define \lambda, m and S?

(Eq. between 114-115): q_\lambda is initially presented with one argument, then two. If the intention is for q_\lambda to represent all approximate distributions  (in this case, p(f, u) ), please state this explicitly. It is natural to assume that q_\lambda(.) is just the Gaussian pdf (as defined in line 110), which creates confusion.

(line 116): k(.,.) as a function also needs to be explicitly defined.

Section 3:

(line 123): This should probably refer to the equation between lines 114-115, not Equation 3.

(Eq 4): Please define \Lambda.

(Eq 6): Please use double integrals to indicate nested integration. Additionally, I don't understand the emergence of l(D_t |f) and the normalizing constant Z. I was under the impression that l(.|.) is the Gaussian likelihood as defined in line 113. Notice that D_t, as defined in line 92, contains both x and y.

(line 147): 'corresponds to'

(line 165): Please clarify where exactly the requirement for the utility to be strictly positive plays a role or comes from.

(lines 173-175): Would using Gauss-Hermite quadrature for EI make EULBO significantly slower than ELBO?

(line 189): What is y^{+}_t and what its relation to y_t?

Question:

Because you formulate the BO queries and GP approximation as a joint optimization problem without necessarily reducing the dimensionality of the problem, how much slower is this method compared to ‘approximation-ignorant’ BO?
limitations: limitations: The authors are open about the limitations of their work, which they leave for future exploration. Specifically, the increased computational complexity seems to be a primary limitation. It would be beneficial if the authors could quantify this increased complexity, perhaps through reporting empirical runtimes or other relevant metrics.
paper_summary: paper_summary: The authors present an extension of the sparse Variational Gaussian Process framework (SVGP) that is better suited for Bayesian optimization (BO) tasks. This is achieved by a new optimization criterion (EULBO) which aims to optimize the parameters of SVGPs s.t. the training data is fit well while achieving a better performance in selecting new solution candidates. Two concrete instantiations of EULBO (with Expected Improvement and Knowledge Gradient) are derived. An extensive empirical evaluation is conducted to demonstrate the effectiveness of the proposed approach.
strength_weakness: strength_weakness: ## Methodology
- there are some details unclear/missing regarding figure 1:
	- which points were used as inducing points?
	- why is the uncertainty relatively constant at the right side of the plot (left plot)? Depending on the choice of the inducing points, I'd expect the predicted variance (and thus the EI) to increase towards the right edge of the plot
	- what was the exact problem definition leading to Fig. 1? It would be good to see the exact parameterizations for reproducibility as it would help to gauge whether the motivational example is rather cherry-picked or whether the shown issue manifests across different parameterizations
- missing $t$ in Eq. 2? shouldn't it be $\pi(f | \mathcal{D}_t)$?
- in the equation after line 221 it seems that the sum does not depend on $j$. shouldn't it be $\mathbf{x}_j$?

## Experiments
- why were the baselines only applied in conjunction with EI? (see line 259)
- it seems that, depending on the choice of the objective and method, there are quite significant differences in the variance of EULBO. For example, in Fig. 2 "Media Molecules 1", TuRBO + ELBO EI has much higher variance than TuRBO + EULBO EI. Do you account this observation to EULBO?
questions: questions: see weaknesses
limitations: limitations: the authors discussed the limitations of their approach in a detailed waysssssssssssssssssssssss

Comments: title: Answer to question about motivation
Rebuttal Acknowledgement
Regarding your question about motivation
Answer to Rebuttal
comment: Thank you for your response, which resolves my concerns and raises some very interesting and insightful discussion on the svgp.
Thank you very much for responding to my review. I am happy with the response and I believe the paper is worthy of publication.
Thank you for the detailed answer, this clarifies my question.
I would like to thank authors for their response and addressing my concerns. I am happy to keep the assessment of their work the same.
This is a good question that we're not sure we have enough data to confidently answer. The benchmark problems we consider in this paper are a mix of problems that some were originally introduced in BO papers using non variational GPs, and others using them. For the ones that have already used variational GPs, there is obviously some survivorship bias where, if SVGPs catastrophically failed on those problems, we wouldn't have seen them. The other problems that originally did not use variational GPs like Rover give us at least some evidence about the performance degradation when switching to SVGPs with ELBO, and it at least does not seem to be catastrophic. We haven't encountered a problem where SVGPs with the ELBO have totally failed compared to known exact GP results, in part because they are often used in situations where evaluation budgets mandate them.
Thank you very much for addressing the points I raised!

The rebuttal clarified my points.

There is only one question left from my side regarding motivation: Could you provide intuition on whether there are certain circumstances in which using ELBO instead of EULBO heavily leads to performance deterioration of BO with SVGPs? In the experiments, one can see that there are problem instances in which EULBO leads to significantly better results than ELBO. Still, in other problem instances, the difference is not that large.

Overall, I'm happy to keep my score.