Scalable Bayesian Optimization via Focalized Sparse Gaussian Processes

Authors: yunyue  wei
vincent  zhuang
saraswati  soedarmadji
yanan  sui
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 03:27:24

Abstract: Bayesian optimization is an effective technique for black-box optimization, but its applicability is typically limited to low-dimensional and small-budget problems due to the cubic complexity of computing the Gaussian process (GP) surrogate. While various approximate GP models have been employed to scale Bayesian optimization to larger sample sizes, most suffer from overly-smooth estimation and focus primarily on problems that allow for large online samples.  In this work, we argue that Bayesian optimization algorithms with sparse GPs can more efficiently allocate their representational power to relevant regions of the search space. To achieve this, we propose focalized GP, which leverages a novel variational loss function to achieve stronger local prediction, as well as FocalBO, which hierarchically optimizes the focalized GP acquisition function over progressively smaller search spaces. Experimental results demonstrate that FocalBO can efficiently leverage large amounts of offline and online data to achieve state-of-the-art performance on robot morphology design and to control a 585-dimensional musculoskeletal system.

Main paper: 
 introduction: 
 Bayesian Optimization (BO) is a powerful approach for solving black-box optimization problems, demonstrating notable success in hyperparameter tuning [1] , reinforcement learning [2, 3] , and scientific discovery [4] . The efficacy of BO is attributed to its ability to model the unknown objective function using a surrogate model and to strategically select the next sample position by optimizing an acquisition function. Among the surrogate models, Gaussian Processes (GPs) [5] are usually favored due to their flexibility and robust uncertainty quantification. However, the computation of the posterior GP covariance matrix scales as O(n 3 ) with the number of data points n, which can severely restrict the applicability of BO in handling large datasets. This poses a significant challenge for real-world applications with high-dimensional and heterogeneous function landscapes such as those in robot control, which often necessitate a substantial amount of data to adequately explore the vast search space. To extend the scope of BO to accommodate larger datasets (from long-horizon online trials and/or pre-collected offline datasets [6] ), it is imperative to employ surrogate models that offer enhanced computational efficiency.Using sparse GP models is a popular method for reducing the computational cost of BO. Sparse GPs accomplish this by learning an approximation of the full GP, either by using a subset of data [7] , ensemble of local models [8] , or variational inference [9] . However, classical sparse GP models are typically tailored for regression tasks, and therefore are designed to fit to the entire function landscape. Given limited representational resources, the resulting posterior is likely to be overly smooth, which may negatively impact the performance of BO. This issue is exacerbated in the high-dimensional setting, in which accurately fitting the entire domain is a far more challenging task. As such, several works have proposed strategies to improve BO performance with sparse GP models by focusing promising regions [10, 11] or advanced sparse GP models [12] . However, most of their empirical evaluations are only conducted under large online sample setting in low-dimensional problems with fewer than 20 variables. It is unclear whether existing methods can be generalized to large offline data or high-dimensional setting.In this work, we explore the application of sparse Gaussian processes for optimizing highdimensional problems with large offline (and optionally large online) datasets. We argue that by iteratively identifying key sub-regions of the input space and focusing the modeling capacity on these areas, we can enhance the modeling fidelity of the sparse GP in regions that are most relevant, thereby improving the overall performance of the Bayesian optimization algorithm. To this end, we propose a novel loss function to train a variational sparse GP model (focalized GP) that emphasizes the fitting of local functional landscapes through weighting the training data. Along with focalized GP, we design a hierachical algorithm, FocalBO, to propose sample points via acquisition function optimization across varying scales of the search space. Experimental results demonstrate that FocalBO can improve upon commonly used acquisition functions in optimizing heterogeneous functions and can effectively utilize large offline datasets for efficient high-dimensional optimization. Furthermore, we showcase that FocalBO can efficiently optimize a policy with 585 parameters to control a musculoskeletal system, leveraging both offline and online data. To the best of our knowledge, FocalBO is the first sparse GP-based Bayesian optimization algorithm capable of efficiently optimizing high-dimensional problems under both large online sample and large offline data settings.Our main contributions: 1) We design FocalBO, which employs a hierarchical acquisition optimization strategy to achieve efficient optimization over high-dimensional problems with heterogeneous structure with limited representation capability. 2) Experimental results demonstrate the superior performance of FocalBO in leveraging large offline datasets for online optimization, and its capability to optimize high-dimensional musculoskeletal system control problems involving over 500 variables. 
sparse gaussian processes: 
 Scaling Gaussian processes to large datasets is an important topic [13] . It can be broadly divided into global approximation strategies and local approximation strategies.Global sparse GPs perform distillation over the whole dataset to approximate the expensive full covariance matrix with a sparse representation. Several methods aim to choose a subset of representative training points from the whole dataset, and use the corresponding covariance matrix in place of the full covariance [14, 7, 15, 16] . Sparse kernels aim at removing uncorrelated entries in the full covariance to obtain a compact matrix [17, 18, 19, 20] . Sparse approximation methods use inducing variables to learn a low-rank representation of full covariance matrix [21, 22, 23, 24, 8, 9, 25, 26, 27] . Stochastic variational Gaussian process (SVGP) is a popular sparse GP method which employs variational inference to learn inducing variables and kernel hyperparameters jointly and enable training using stochastic gradient descent from mini-batch data [25] . Recently, nearest neighbor information has also been used to further improve the scalability of sparse GP over massive amount of data [28, 29] .In contrast, local sparse GPs divide the entire dataset and employ local GPs trained from different data subsets to approximate the full GP. For a given test set, the prediction can be extracted from one of the local GPs [30, 31] , mixture of GPs [32, 33] or product of GPs [34, 35] .Recent works have proposed modifications to sparse GPs for Bayesian optimization. Sparse GP has been used to determine the search region where local GPs are used to determine the next samples [36] . Weighted-update online Gaussian processes (WOGP) was developed to select a subset of training points to approximate high performing regions of the input space [10] . IMP-DPP is motivated by a similar observation and uses a weighted Determinantal Point Process to select training points as inducing variables for the SVGP [11] . However, their proposed selection strategies require sequentially evaluating every training point, which can be computationally very expensive with large offline datasets. Combining SVGP with Thompson sampling has the same order of regret as standard Thompson sampling method [37] . Online variational conditioning (OVC) was proposed to efficiently conditioning SVGPs in an online setting, enabling using look-ahead acquisition functions [38] . Vecchia approximation of GP was also applied [39] for Bayesian optimization, with improved performance compared to prior works [12] . A concurrent work [40] aims at improving the acquisition optimization performance based on target-aware Bayesian inference [41] .Besides sparse GPs, Neural network [42, 43] and random forest [44] can also be used as BO surrogate model to circumvent the cubic complexity of GP. Ensemble Bayesian optimization utilizes the addictive function structure and uses ensembles of addictive GPs in parallel to achieve scalability [45] . Trust Region Bayesian optimization (TuRBO) and its variants uses exact GP to optimize over local regions, and employs a restart mechanism to achieve large number of evaluation, which is a representative line of works in high-dimensional Bayesian optimization [46, 47, 48] . TuRBO can also be combined with sparse GP models to further enhance the scalability [49, 50] . 
bayesian optimization: 
 For an unknown objective function f , Bayesian optimization aims to solve max x∈X f (x) over input space X ∈ [0, 1] d . BO mainly consists of two components: a surrogate model to approximate the objective function, and an acquisition function a to decide the next sample position based on surrogate model.Gaussian process is a commonly used surrogate model. Consider a given dataset D = (X, y) where X = (x 1 , ..., x t ) are input locations and y = (y 1 , ..., y t ) are associated noisy observations of f (X). We assume the observation noise to be independent Gaussian, i.e. y i = f (x i ) + η, η ∼ N (0, σ 2 ). Using GP with kernel function K, the function distribution f * at test positions X * = (x * ,1 , . . . , x * ,t * ) T is a multivariate Gaussian:EQUATIONwhere K is the covariance matrix between subscript inputs. With the posterior distribution given D, the next sample point is the maximum position of the acquisition function:x t+1 = max x∈X a(x|M t ), where M t is the GP model fitted on dataset collected at time step t.Common-used choice of a includes upper confidence bound(UCB, [51] ), expected improvement (EI, [52] ) and Thompson sampling (TS, [53] ). The inner optimization problem is usually solved by grid search, evolutionary algorithms [54] , or gradient-based methods [55] . When the online sample budget is large, batch optimization is commonly used to evaluate multiple inputs in parallel [56] .For predictive distribution conditioned on given dataset of size t, the computational complexity of exact Gaussian process is O(t 3 ) for each test position due to the inverse of the covariance matrix K XX , which is expensive for large scale datasets with more than a few thousand points. A common used strategy is to approximate full GP regression using sparse GPs. In sparse GP, m ≪ t inducing variables u = (u 1 , . . . , u m ) T characterized by inducing inputs Z = (z 1 , . . . , z m ) are introduced to approximate the covariance matrix of the full GP. In this section, we focus on sparse GP derived from variational inference.Variational GP [9] considers the joint latent priorEQUATIONwhere f = (f (x 1 ), . . . , f (x t )) T . A variational distribution q(u) = N (u | m, S) is used to approximate the posterior over inducing variables using the exact conditional distribution of f given u, that is, q(f , u) = p(f | u)q(u). The posterior of f can be computed by marginalizing u with analytic form:EQUATIONwhereA = K -1 ZZ K ZX .The variational parameters Z, m, S are optimized by maximizing the Evidence Lower Bound (ELBO) which can be written in the following formulation [25] :EQUATION)where KL[• ∥ •] is the KL divergence between two distributions. The ELBO breaks into a data likelihood term which factorized over training data, and a KL divergence term which can be computed in closed form. The factorization over data allows optimization via stochastic gradient descent (SGD), reducing the computational complexity to O(m 3 ). Prior studies about variational sparse GPs are mainly designed for regression tasks, where the goal is to fit global training data distribution. In Bayesian optimization, the next sample is determined by the predictive function distribution over test positions. 
focalized gaussian process for bayesian optimization: 
 Gradient-based and evolutionary-based acquisition function optimization methods employ local search from random starting points to find a local optimal of the acquisition function. Recent works also scale grid search-based optimization to high dimensional space by restricting the search space within local sub-regions [46, 47] . All the above procedure would benefit from an accurate estimation over sub-region of the input space. Therefore, a sensible way to improve BO performance is to allocate limited computational resources to obtain better prediction over specific search regions instead of the entire input domain.We define the search region as the region where the acquisition function is optimized on, which is an axis-aligned hypercube with length l = (l 1 , • • • , l d ) T centered at c:EQUATION)When l = (1, • • • , 1) T and c = (0.5, • • • , 0.5) T , the acquisition optimization is performed over the entire input space X , as commonly-used in vanilla BO algorithms. In the rest of this section, we first present the derivation of focalized loss function to improve GP prediction over the search region. Then we demonstrate how to incorporate our proposed GP model into Bayesian optimization.We recall eq.1 and rewrite the mean estimation µ t (x * ) and variance estimation σ t (x * ) for each test position x * :EQUATION)where kij is the (i, j)-th entry of [K XX + σI] -1 . From eq. 6 we can observe that the mean estimation at x * is a linear combination of observation y multiplied by k(x * , x), and the reduction of variance is a quadratic form of the covariance between x * and training points.Both estimation can be written as linear summations of constant values with kernel function as weight. As mentioned in prior works [57] , data points far from the test positions have a vanishingly small influence on the predictive distribution with commonly used kernel functions. Utilizing this observation, we propose to weight the data likelihood term using the kernel function to focus training over points that contribute to the prediction of the search region:EQUATION)We use the maximum covariance of x i to positions in the search region as the corresponding weight to filter out points that have marginally influence to the search region during GP training. In this way, the model can selectively utilize the training data to achieve good local prediction. When using a popular kernel functions such as RBF or Matern kernel, the maximum kernel value is equivalent to finding the nearest point in the search region, which can be easily calculated when the region boundary is axis-aligned as defined in eq. 5.We additionally regularize the sum of weights to make the model focus on improving prediction over search region:EQUATIONwhere |X ∈ S c,l | = t i=1 1 xi∈S c,EQUATION)Compared to the original ELBO loss in SVGP, our proposed function maintains the same computational complexity and does not introduce additional hyperparameters. Our ELBO also reproduces eq. 4 when considering to predict the entire input space X . During the model training, both GP hyperparameters and variational parameters are jointly optimized to obtain focalized GP for Bayesian optimization. Figure 1 shows a comparison of focalized GP and SVGP over 1d functions sampled from GP. While SVGP can only able to vaguely predict the function, focalized GP accurately delineate the function landscape within search region by training with the focalized loss. Our proposed GP model is sensitive to high-performing positions within the search space which contribute to better acquisition optimization. We also systematically compare the GP prediction performance in Appendix B. 3 , where our GP model trained from focalized ELBO consistently achieves good prediction on small size of search space.Theoretical implications of focalized ELBO. Our focalized ELBO can be interpreted as a soft variant of training a local approximation over datapoints that lie within the search region. Here, we illustrate how local approximations can substantially reduce the KL divergence of the approximate posterior over the search region, and discuss the effects of tighter approximations on BO regret bounds. We focus on providing general theoretical intuition rather than deriving precise bounds due to the lack of existing convergence guarantees for ELBO maximization in the general setting.Suppose that we know the optimal point lies in some small sub-region of X that contains N ′ << N training points. Corollary 19 in [58] shows that given a squared exponential kernel and some assumptions on the inducing point selection, for a fixed number of inducing points the KL-divergence upper bound scales super-quadratically in the number of training points. Hence, fitting locally can yield much tighter approximations than fitting globally (e.g. SVGP).Next, we consider the impact of the KL approximation error on the optimization regret. Proposition 1 in [58] states that the gap between the means of the approximate and exact posteriors is upper bounded by O(σ √ γ), where γ is an upper-bound on the approximation KL-divergence. This has an immediate impact on the regret -for example, when GP-UCB [51] is combined with sparse GPs, the confidence bounds must be enlargened by an additive √ γ factor to account for the approximation error. Because the regret bound scales with √ β T where β T is the maximum confidence interval coefficient, having a large approximation error can arbitrarily scale the regret incurred by the algorithm. In order to achieve no additional regret order, the additional approximation error noise must be uniformly bounded (Assumption 4 in [37] ). Although focalized GP cannot guarantee a constant bound, it still directly reduces the regret of the algorithm, where we empirically investigate in Appendix B.1.Input Initial Dataset D 0 , Inducing Vari- able Size m, Batch Size B 1: H ← 1 2: for t = 1, 2, • • • do 3: {x t,i } B i=1 , {h t,i } B i=1 ← FocalAcq(D t-1 , H, m, B) 4: Observe {y t,i } B i=1 = {f (x t,i ) + η} B i=1 5: D t+1 ← D t ∪ {(x t,i , y t,i )} B i=1 6: i best ← argmax i∈1,••• ,B y i 7: if h i best < H then 8: H ← H -1 9: else H ← H + 1 10:end if 11: end for One advantage of focalized GP is that it can be easily integrated into existing BO algorithms. To further leverage the strong local modeling properties of focalized GP, we design FocalBO, a hierachical acquisition optimization framework described in Algorithm 1.At each BO iteration, FocalBO iteratively optimizes the acquisition function over a progressively smaller search region via focalized acquisition function (FocalAcq) as shown in Algorithm 2. The first depth of acquisition optimization starts with the entire input space X with l = (1, • • • , 1) T and c = (0.5, • • • , 0.5) T (line 1). We train specific focalized GP base on the search region at each round of acquisition optimization (line 4-5). Our framework is compatible with any acquisition function that extracts instant posterior information from the GP and is optimized within pre-defined search region. After one round of acquisition function optimization, the search space length l is halved to focus on a smaller search region centered at current best position x best (line 6-7). In this way we can obtain a more accurate model for decision making, and also relieve the over-exploration problem when the problem dimension is high [59] . One batch of inputs is proposed at each round of optimization, and the final decision is sampled from all proposed inputs via Softmax distribution over their corresponding acquisition function values (line 9). Our hierarchical optimization strategy enables collecting candidates from both global sparse estimation and local focalized prediction, achieving balance between exploration and exploitation with constrained computation power.The optimization depth H in FocalAcq controls the degree of utilizing local information from current best position, where the GP estimate variance decreases with the shrinkage of search space. The best-performing optimization depth is likely problem-dependent (e.g. highdimensional functions may require higher optimization depths). Therefore in FocalBO, we propose to automatically adjust the optimization depth according to the instant optimization performance. At the beginning of the optimization, we initialize the optimization depth as 1, indicating global search of the input space (Algorithm 1, line 1). Then we keep track of the depth where the proposed positions are sampled from. If the depth of the best point in this round is less than the current optimization depth H, we reduce H to encourage exploration of the input space, otherwise we increase H for better exploitation of x best (line 6-10).Input Dataset D t-1 , Optimization Depth H, Inducing Variable Size m, Batch Size B 1: l ← (1, • • • , 1) T , c ← (0.5, • • • , 0.5) T 2: Select current best point x best from D t-1 3: for h = 1, • • • , H do 4:Train M h t using L 2 given S c,l 5:{x h t,i } B i=1 ← argmax x∈X * a(x|M h t ) 6: l ← l/2 7: c ← x best 8: end for 9: return {x t,i } B i=1 , {h t,i } B i=1 ∼ P (i = i ′ ) ∝ exp a(x h ′ t,i ′ |M h ′ t ) H h=1 B j=1 exp a(x h t,j |M h t )Our proposed framework is orthogonal to TuRBO-M [46] , but bears some similarities in searching over multiple sub-regions and adaptively adjusting the search region. Our algorithm differs in that TuRBO-M constructs equal-sized trust regions and fits independent Exact GP using separated dataset, aiming at searching for different local optima in the search space. By contrast, the search region in FocalBO is constructed with different sizes to make decision based on both global and local information. Our framework allows data sharing across search regions, and the use of focalized GP helps to accurately estimate local region with limited representation. Additionally, FocalBO does not introduce extra hyperparameters. Finally, we demonstrate in Section 5 that TuRBO is complementary to FocalBO in optimizing high-dimensional problems. 
experiments: 
 In this section, we extensively evaluate FocalBO over a variety of tasks. We first use synthetic functions to showcase the compatibility of FocalBO in improving commonly-used acquisition functions. Next, we consider the online optimization of robot morphology design that is additionally given a large offline dataset. We also show that FocalBO is able to optimize very high-dimensional musculoskeletal system control with both a large offline dataset and a large number of online budget. Finally we dig deeper into FocalBO to analyze how each of its components contributes to superior optimization performance.We compare FocalBO with representative sparse GP models used for Bayesian optimization, including SVGP [25] , WOGP [10] , and Vecchia GP [12] . We only run WOGP on synthetic functions due to its extremely low speed in dealing with the datasets in the remaining tasks.The number of inducing variables in sparse GP models is set as 50 for synthetic functions and as 200 for other tasks. The optimization performances are shown as mean ± 1 standard error for all considered problems over 10 independent trials.We select Shekel and Michalewicz as the test functions, which are heterogeneous with both smooth and rigid regions. We also sample functions directly from Gaussian processes to evaluate algorithm performance under full BO assumption. For each function, we choose to use different acquisition functions to optimize: TS optimized by grid search, EI optimized by analytic gradient, and probability of improvement (PI) optimized by Monte Carlo gradient [55] .Optimization performances are shown in Figure 2 . We observe that FocalBO significantly improves the performance of all acquisition functions compared to SVGP, and is able to consistently achieve top-tier performance over all problems. In Michalewicz function where a large fraction of the input space is flat, all baselines tend to increase the noise estimation to maintain a stationary prediction, while focalized GP is able to focus on the local search region and successfully optimize the function. Additional experiment with online samples as major data source is shown in Appendix B.2, where FocalBO still maintains comparable or better performance against baselines. We compare FocalBO to several baselines over robot morphology design task from Design-Bench, which provides large offline dataset with an exact function oracle [6] . The goal of the task is to optimize the morphological structure of D'Kitty robot [60] to improve the simulation performance under RL controller. While the benchmark is initially designed for offline model-based optimization (MBO), it can also be used as an offlineto-online BO benchmark. In this task, we use the training dataset with 10,000 points and additionally evaluate 128 points on-thefly with batch size of 4. EI is used as the base acquisition function for better optimizing with small batch size. We also try to combine FocalBO with TuRBO to optimize over the high-dimensional space, with the results shown in Figure 3 . We observe that FocalBO achieves significant improvement from the initial data while other baselines struggle to obtain performance gain, even combined with TuRBO. FocalBO with TuRBO effectively extracts information from large offline dataset and is the first GP-based method to achieve top-tier performance reported by prior MBO works [61] .We further apply FocalBO to control a human arm musculoskeletal system [62] for the task of pouring liquid into a cup, as shown in Figure 4 (a). To control the musculoskeletal system, we optimize a linear policy π ∈ |A| × |O|, where |A| = 5 and |O| = 117 are the corresponding action and observation dimensions. The action dimension has been reduced from individual muscles to synergetic groups of muscles by applying principled component analysis to sampled action data from an RL agent (Appendix A.6). Although the original control dimension is reduced, the remaining 585-dimensional input space is still very high for existing highdimensional BO algorithms. Therefore we consider a large offline-online setting, where we randomly sample 2000 points from the input space to serve as the offline dataset, and set the online budget as 3000 with batch size of 100. We use Thompson sampling as the base acquisition function. Figure 4 (b) demonstrates that FocalBO outperforms other baselines, achieving higher maximum reward and faster convergence speed. Our supplementary video shows that the optimized policy is able to perform well on the task, demonstrating the successful application of FocalBO to high-dimensional control problems. To understand the reasons behind FocalBO's superior optimization performance, we investigate the optimization depth in FocalBO, which is the central component of the method. Figure 5 (a) shows the evolution of optimization depth over different problems, where FocalBO is able to adapt the optimization depth according to different function structure. For Shekel and musculoskeletal model control where the promising regions are distinct, the optimization exhibits an increasing trend to exploit current best points, while for other problems the depth tends to converge at a fixed level. Figure 5 (b) shows the sources of proposed batches during the optimization of musculoskeletal system control. Overall the samples exhibits clear trend from exploration to exploitation over high-dimensional input space. Our hierarchical optimization strategy enables flexibility between exploration and exploitation.  
conclusion: 
 In this paper, we propose FocalBO, which uses a hierarchical acquisition optimization strategy equipped with focalized GP model to scale Bayesian optimization to problems with large offline datasets and/or a large number of online samples. Despite limited representation capability, FocalBO consistently improves various acquisition functions in optimizing heterogeneous functions, and adeptly leverages large offline dataset for efficient optimization over robot morphology. Under the large offline-to-online optimization setting, FocalBO achieves stable high-dimensional control of human musculoskeletal model with over 500 parameters. Ablation studies over the algorithm components further verify the principled design of FocalBO. Future work may include theoretically analyzing FocalBO, and applying the method to more complex problems, such as large-scale parameter tuning and whole-body human musculoskeletal system control. r = 50r pos * 10r ori + 10r reach + r lift -r act -5r done (10) where r pos encourages the bottle near the target position, r ori encourages the bottle near the target orientation, r reach encourages the hand to grab the bottle, r lift encourages the hand to lift the bottle, r act penalize the overall muscle activation, r done penalize the early ended episode due to dropped bottle or hand outside of pre-defined range.We trained a Soft Actor-Critic (SAV) [66] agent for 6M timesteps to collect task-related muscle activation data, and use principled component analysis to reduce the action dimension from 81 to 5. 
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)? 
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?Answer: [NA]Justification: The paper does not use existing assets.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset. • The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• 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. 
 
 Appendix: B.1 Theoretical implications of sparse GP approximation.In Figure 6 , we also empirically measure our claim that Focalized GP can significantly reduce approximation error on the search region. We sampled 8000 training points from 2d GP functions to train focalized GP and SVGP. Over different size of the search region, we compare the KL divergence of the GP posterior prediction over search region between sparse GPs and the exact GP. We observe that the KL divergence between focalized GP and exact GP is consistently smaller than that between SVGP and exact GP, implying tighter approximation to the exact GP over local region.Figure 6 : KL divergence between sparse GPs and exact GP. Results shows the mean and one standard error, averaged over 50 independent trials.While a rigorous regret bound is hard to derive, we conduct an empirical study where we directly compare the optimization performance between focalized GP and SVGP when combining with TuRBO. In this way we can eliminate the influence of hierarchical acquisition optimization. The optimization performances are shown in Figure 7 . We observe that focalized GP outperforms SVGP on both high-dimensional problems, which empirically demonstrates our theoretical implications that Focalized GP contributes to reducing regret.We empirically investigate this in Figure 8 (a), which compares different ways of selecting the search region center by measuring the distance from the search region center to the global optima. We observe that current best point consistently is the closest to the global optimum, which validates this design choice.For the experiment above, we sampled 2d functions from GPs with Matern 5 2 kernel and lengthscale of 0.05 (representing rigid functions), and selected the best point over unifromly sampled 10,000 points as the global optima.A sparse GP is already more explorative than using the full GP, since the smaller representational capacity leads to smoother posteriors. In Figure 8 (b), We demonstrate this empirically below, where we measure the pair-wise distance of 100 Thompson sampling points under exact and SVGP (with 50 inducing points). We observe that sparse GP actually samples more diverse sets compared to exact GPs, i.e. exhibiting more exploration. Therefore, We choose Ackley and Hartmann, which are common-used test functions for BO community. We use the similar optimization setting in [12] . The optimization performances are shown in Figure 9 , where FocalBO is still able to achieve comparable or better performance when online samples dominates the data source. We use two common-used synthetic functions, Ackley and Rastrigin, to analyze the the GP predictive performance of focalized GP compared with Exact GP and SVGP under different search region size l and different inducing variables number m. We show the negative log likelihood (NLL) and root mean squared error (RMSE) in Figure 10 . The results shows that focalized GP outperforms both Exact GP and SVGP in terms of both NLL and MSE when the search space size is lower than 0.5. In Rastrigin function where Exact GP achieves similar performance as SVGP, focalized GP is still able to accurately predict the local search region over different choice of inducing variable numbers. We also show in Figure 11 that the regularization term L reg is indispensable to the training of focalized GP to achieve good local prediction. We run the original TuRBO implementation (with exact GP and Thompson sampling) and TuRBO with nearest neighbor GO model on both robot morphology design and human musculoskeletal system control task (Figure 12 ). We observed that FocalBO outperforms TuRBO on both tasks with smaller computational cost. The reason of TuRBO's poor performance may be that it cannot quickly adapt over the search space when the online evaluation budget is small. [58] ). Let k be a squared exponential kernel. Suppose that N real-valued (onedimensional) covariates are observed, with identical Gaussian marginal distributions. Suppose the conditions of Theorem 13 are satisfied for some R > 0. Fix any γ ∈ (0, 1]. Then there exists an M = O(log(N 3 /γ)) and an ϵ = Θ(γ/N 2 ) such if inducing points are distributed according to an ϵ-approximate M-DPP with kernel matrix K f f , Lemma 2. (Proposition 1 in [58] ). Suppose 2KL[Q ∥ P ] ≤ γ ≤ 1 5 . For any x * ∈ X , let µ 1 denote the posterior mean of the variational approximation at x * and µ 2 denote the mean of the exact posterior at x * . Similarly, let σ 2 1 , σ 2 2 denote the variances of the approximate and exact posteriors at x * . Then,EQUATIONLemma 3. (Assumption 4 in [37] ). (quality of the approximate prediction). For the approximate μt , the exact µ t and σ t , and for all x ∈ X ,EQUATION)where 0 ≤ c t ≤ c for all t > 1 and some constant c ∈ R• 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). 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).Justification: There is no risk about the possible misuse.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.https://botorch.org/https://gpytorch.ai/https://github.com/ermongroup/bayes-opthttps://github.com/feji3769/VecchiaBO/tree/master/code/pyvecchhttps://botorch.org/tutorials/turbo_1https://github.com/brandontrabucco/design-bench 


Decision: This paper proposes an approach for Bayesian Optimization in the large data regime that uses an SVGP with modified ELBO loss to improve quality of the surrogate in “promising” regions of the search space. Reviewers generally agree that the problem is relevant, and that the empirical performance of the approach compared to baselines is compelling (based also on the additional results provided by the authors during the rebuttal period). One of the main criticisms is that the methodology feels somewhat ad-hoc and that the theoretical analysis is underdeveloped. The former criticism is valid, but the approach is also simple, reasonably intuitive, and effective. There are some connections to other literature that I think are missing, in particular on the work of target-aware Bayesian inference (https://www.jmlr.org/papers/v21/19-102.html) whose goal is precisely to modify the inference procedure of the Bayesian model to improve its performance on the downstream task (in this case acquisition function optimization). I also should note that the following contemporary work on Approximation-Aware Bayesian Optimization using sparse GP approaches is very related https://arxiv.org/abs/2406.04308 (as this is only a recent preprint I am considering this as a miss of the authors, this is rather a heads up to help with improvements to the paper).

Overall, I think the quality of this paper is somewhat borderline - However, it’s one of the few that considers modifying the surrogate model for Bayesian Optimization in a way that targets the downstream task of better optimization performance, which in my opinion is under-studied and so this makes it a relevant contribution to the field. Hence I recommend it be accepted, but ask the authors to incorporate the constructive reviewer feedback and the additional empirical results.

Reviews: review: paper_summary: paper_summary: This paper focuses on scalable Bayesian optimization, where the training data size is large and the search space is high dimensional. This work proposes a new kind of sparse GP model by designing a variational loss function that allows for adaptively focusing on interesting regions throughout the optimization process. By coupling this GP model with common acquisition functions such as EI, TS, PI, the authors show in experiments the effectiveness of this newly developed GP model. The authors also provide comprehensive algorithm analysis to further decipher the exploration vs. exploitation behavior of the algorithm.
strength_weakness: strength_weakness: Minor error:
- L130, "more that" -> "more than"
questions: questions: - One concern is that the search region is always centered around the current best position. In the case when the best point from initially sampled points is far away from the true global optimum, and the landscape of objective function is not very smooth, then the algorithm will have a hard time to attain the global optimum quickly, because Focalized GP reduces to regular sparse GP in the broader search region and suffers inaccurate prediction. This issue might explain in Figure 3 why "FocalBO TuRBO" has a better performance than "FocalBO EI" (because TuRBO divided up the search space). I wonder if there is any better idea in setting the center of search region to encourage more exploration.
limitations: limitations: The discussion of limitations is absent in the paper, and it is better for the authors to discuss limitations of this work.
paper_summary: paper_summary: Operating in the context of Bayesian optimization, the authors propose to train a surrogate model which focuses on a specific sub-region of the input space by weighting the log-likelihood contributions of each datapoint relative to their distance to that region. They also propose an algorithm for choosing this subregion over the course of optimization, which they demonstrate in numerical experiments.
strength_weakness: strength_weakness: The Lreg term of equation 8 seems to come out of left field, and the implications of adding it to the cost in equation 9 are insufficiently discussed; I say more about this in the Questions section. 

Figure 1 is not super convincing: it's clear that with a different hyperparameter/inducing point setting (perhaps with admittedly more inducing points), that SVGP would be able to interpolate as well.

The theoretical implications section is underdeveloped, containing neither new math nor computational experiments verifying that the claims made there hold up empirically.

Regarding novelty level, I would assess it as medium as it falls within these existing frameworks of trust-region style optimization spawning from Erikkson et al as well as work about tailoring inducing point methods to consider specific parts of the input space (a different implementation but a similar spirit is [1]), though the combination to Bayesian optimization is in my view novel.  

[1] Cole, D. Austin, Ryan B. Christianson, and Robert B. Gramacy. "Locally induced Gaussian processes for large-scale simulation experiments." Statistics and Computing 31.3 (2021): 33.


Only for the benefit of the authors, I am attaching here some grammar/clarity issues which did not affect my scoring of this article:
201) shows an comparison -> shows a comparison
232) directly reduces
246-248) sentence spanning these two lines needs attention.
251) from GP -> from the GP
298) both large offline dataset and large number -> both a large offline dataset and a large number
304) inducing variable number -> number of inducing variables
374) "consistently augments" does not make sense in context; I'm not sure what you're trying to say.
questions: questions: I understand that Lreq (Eq 8) is given such that it is divided by the number of points in the subregion and has one subtracted to it such that its minimum possible value is zero, and since the w's for any points inside S is constant, we are effectively minimizing the sum of kernel functions evaluated at extra-S training points with the projection of that point onto S. Why is this the right thing to do? It's going to bias lengthscale estimates to be small relative to a GP fit to all the data or even one fit only to data inside the subregion. Why is this desirable, and why is Lreg the right way to operationalize this? 

Am I correct that none of FocalBO Turbo, SVGP Turbo, or Vecchia TurBO represent the actual algorithm originally proposed my Eriksson et al, but rather are re-implementations of TuRBO's strategy in other approximate GP frameworks? 

I think there are some interesting ideas here, but I regretfully am presently suggesting rejection as a result of the limited scope of the competitors in the numerical experiments which in my view is required of a methodology paper with this medium level of conceptual novelty appearing in Neurips. I think this is a solid paper but does not yet meet the high bar of a Neurips article.
limitations: limitations: You say that you don't introduce any additional hyperparameters, but it might be more accurate to say that you provide defaults for additional hyperparameters: you set the relative strength of Lreg to the inverse of the number of points in the search region (which makes sense for making the minimum of this objective zero, but it's not clear what implications the scaling has) and you have a specific heuristic for choosing the optimization depth. Why is it reasonable to believe this would still be the right thing to do on some other black-box test function outside of the five you've presented here? 

I think the problems in Section 5.2 and 5.3 are really interesting problems.
But to help compare against prior work, it would be really important to compare against TuRBO itself: though TuRBO uses an "exact" GP, they use numerical analysis tricks to scale to large datasets; (glancing at the paper,they do 20,000 evaluations on a 60D problem). So it seems like TuRBO as originally presented by Eriksson et al should have been a competitor in this simulation study.
paper_summary: paper_summary: The paper presents a new method for scaling Bayesian Optimisation to large datasets by employing sparse GP in a way that facilitates BO.  The proposed method of Focal BO fits a sparse GP at different scales, optimises the acquisition function for each of them and selects the next point to query by sampling with probability proportional to its acquisition function value at its scale. The proposed algorithm is evaluated at a number of benchmarks, where a large number of data points is available.
strength_weakness: strength_weakness: - Since the paper is mostly empirical I would expect a bit more baselines. The simplest way to scale TuRBO to large datasets is simply to keep $N$ points closest to the optimal solution. It would be great to see how this naive "keep closest N TuRBO" baseline compares to FocalBO TuRBO.

- "FocalBO with TuRBO (...) is the first GP-based method to achieve top-tier performance reported by prior MBO works." - this is a strong statement without much justification. Please either directly compare with MBO baselines or quote the exact numbers achieved by prior work on the exact same setting as considered in the paper.
questions: questions: -  In Equation 8, why is there a $ - 1$ at the end? Since it is a constant it should not affect the optimisation process at all. 

- I cannot understand why authors refer to ELBO as "ELBO Loss", since due to how it is defined, I imagine you would rather want to maximise it (as it corresponds to expected likelihood). Are you maximising the loss function? In that case, I do not think it is proper to call it that name.
limitations: limitations: Within the paper authors do not discuss the limitations (or at least I cannot find them).

Comments: title: Rebuttal (4/4)
Rebuttal (3/4)
Rebuttal (2/4)
comment: Thank you for your feedback. Below, we clarify the novelty and impact of our work.

**The novelty of our work**

**We design the first sparse GP model that improves acquisition function optimization**. Our design of focalized GP enables **better local estimation, aiming at improving acquisition function optimization during Bayesian optimization**. We adopt variational inference to derive focalized GP, which is **capable of performing joint posterior inference and sampling given a set of test points just like exact GP**. Therefore our proposed GP model is automatically compatible with any acquisition function that is used for exact GPs. (Figure 2). Previous work which allocate inducing variables based on test points can not directly be incorporated into BO framework. Reviewer SUei also recognized that focalized GP is **fundamental different** from previous works.



**We propose the first scalable BO algorithm that is capable of utilizing large offline dataset for optimization.** Our design of hierarchical acquisition function optimization  with focalized GP **searches promising positions over different scale of search region during one BO iteration, enabling making decision based on both global and local information under restricted computational budget**. Besides, FocalBO also demonstrates **superior optimization performance in large online setting** in our human musculoskeletal system control task. 

Regarding the Reviewer SUei's concerns of the scalibility, we emphasize that **our used focalized GP in FocalBO share the same $\mathcal{O}(m^3)$ computational complexity as SVGP**, while the complexity of TuRBO is $\mathcal{O}(n^3)$ for using exact GP. We addressed this point by investigating the restart mechanism of the original TuRBO and demonstrating a failure case of TuRBO under low computational resources. 
Our numerical experiment also demonstrates that TuRBO needs significantly more GPU memory than FocalBO even under approximated linear algebra computation (78.9 GB v.s. 7.8 GB). Our empirical evidences demonstrate that FocalBO is more scalable than orginal TuRBO. 

**The impact of our work**

We argue that our addressed problem setting is not niche. **The problems we tackle in the experiment section includes both large offline and large online setting, where FocalBO consistently demonstrates the superior performance.** In the human musculoskeletal system control task, FocalBO outperforms TuRBO with 3000 online evaluations, which we consider to be a substantial number. We emphasize the importance of the large offline setting as no previous scalable BO methods have been able to effectively handle it. We believe our setting is important because the problems addressed by BO are typically expensive to evaluate. Instead of collecting data online from scratch, many problems already have data available from various sources. FocalBO’s ability to utilize larger datasets enables it to tackle challenging problems where the online evaluation budget is limited.
**Continue on scalibility**

We argue that a direct comparison between the original TuRBO and FocalBO is not entirely fair. While focalized GP need to be optimized at every depth, it **share the same $\mathcal{O}(m^3)$ computational complexity as SVGP**. The $H\ll m$ re-optimization rounds do not alter this order. The computational complexity of TuRBO $\mathcal{O}(n^3)$ for using exact GP, and our empirical results demonstrate that using approximate linear algebra computations in GpyTorch does not reduce the computational cost to the same level as sparse GPs. We tracked GPU memory usage during the optimization of the robot morphology design task. The maximum GPU memory usage for the original TuRBO reached 78.9 GB, whereas FocalBO, with an inducing size of 200, used only 7.8 GB. Our emprirical evidence indicates that FocalBO is more scalable than original TuRBO.


**Impact of our work**

We argue that our addressed problem setting is not niche. **The problems we tackle in the experiment section includes both large offline and large online setting, where FocalBO consistently demonstrates the superior performance.** In the human musculoskeletal system control task, FocalBO outperforms TuRBO with 3000 online evaluations, which we consider to be a substantial number.
We emphasize the importance of the large offline setting as no previous scalable BO methods have been able to effectively handle it. We believe our setting is important because the problems addressed by BO are typically expensive to evaluate. Instead of collecting data online from scratch, many problems already have data available from various sources. FocalBO’s ability to utilize larger datasets enables it to tackle challenging problems where the online evaluation budget is limited.
Thanks for providing empirical evidence for the design choice of search region center. I still think this work has adequate contribution to the field, however, regarding SUei's comment on novelty / impact of this work, the author did not offer more evidence, I would revise the score to 7 (Accept).
Thank you for providing detailed feedback and raising the score! 
Below, we address your concerns regarding the scalability of FocalBO.

As mentioned in our novelty rebuttal section, the scalability of TuRBO is primarily achieved through a restart mechanism, which discards all previously sampled data and re-initializes the optimization process. In our reproduction of the experiments from the original TuRBO paper, we found that the average number of evaluations in the rover planning problem (originally reported as 20,000 total evaluations) before restarting was $3198.7\pm649.8$ (mean$\pm$1 std, averaged over 77 restarts). This suggests that the actual data size used for GP in TuRBO is much smaller than data size in our addressed problem. 

Regarding the GP implementation in the original TuRBO, we agree that the conjugate gradient method and Lanczos decomposition in GPyTorch help to accelerate and parallelize computation via GPU. 
To further demonstrate the computational efficiency of our method, we ran FocalBO and the original TuRBO on an NVIDIA GeForce RTX 2080Ti with 11 GB of memory, instead of an NVIDIA A100 with 80 GB of memory. We observed that FocalBO was still capable of optimizing the robot morphology design task, while TuRBO encountered an out-of-memory error during GP training, even when using approximate computing methods. We believe that data size remains an issue in TuRBO, and this may be one of the reasons why TuRBO introduces a restart mechanism.
I’ll reply to all 4 of your rebuttal messages here. 

Regarding Figure 1: Thanks for challenging me on this and I think you’re right. 

Regarding Lreg: Thanks for pointing out this alternative characterization of that cost function. I agree that heuristically this sort of makes sense as a way of dealing with some of the problems that arose. But what I’m suggesting is that it’s not convincing that this will work on general problems. The modification appears ad-hoc and appears not to be motivated by a statistical model, numerical approximation, or otherwise a precise idea of how big that term should be in order to balance against the expected log-likelihood and KL terms.

Regarding the theoretical implications: Thanks for these simulations, they are interesting and help to build understanding of the proposed methodology. 

Regarding novelty of the focal approximation: I agree that there are fundamental differences between LIGP and your work, and I do not know of prior work which uses a “goal oriented” inducing point approximation for BO. I reference LIGP to point out that the idea of teleologically determining inducing points is not novel; it is rather its application to BO which is novel. 

Regarding Scalability: I must for now push back on the idea that the proposed method’s scalability is novel. The approach used in your article for scaling the GP falls into the meta-approach of approximating somehow the GP, and then doing the linear algebra exactly. It’s true that TuRBO does not use this meta-approach. However, an alternative meta-approach used for GPs is to avoid approximating at the statistical/GP level, but then to instead approximate the linear algebra directly. This latter approach is the one taken by Eriksson et al in their original implementation, which leverages pytorch and fancy linear algebra to avoid direct cholesky decomposition. Consequently, looking back at Eriksson et al, it seems like their largest dataset is of size 20,000. If I am not mistaken, the largest study in the proposed article in terms of total sample size is the Robot Morphology problem with a total problem size of 10,128, which is actually less. Furthermore, all but the last 128 of these are acquired offline. A priori, I would imagine that the fact that the ELBO has to be re-optimized at every level of depth makes for some amount of overhead. However, you did mention in your rebuttal that your algorithm had a smaller computational cost than TuRBO, though without providing the numerical results; I understand you may not have been able to fit those results onto the single rebuttal page. But on the basis of the observed experiments, we actually find that TuRBO’s original article contains the larger problems. (Please let me know if I missed or misinterpreted something).

Regarding the comparison with the original TuRBO: Thanks very much for running these experiments, those are very helpful. It’s impressive that your method is able to keep ahead of it on those important applications. 

Regarding The rationality of FocalBO design: As I explained above, I still don’t understand why it would be the case that the Lreg term will always have the right magnitude relative to the Lwll and Lkl terms. But regarding the optimization depth criterion, I think I agree that what you’re doing makes sense.

In conclusion: I thank the authors for the additional experiments and for engaging with my review. I think the direct comparisons against TuRBO are convincing so I’ll raise my score to a 4/borderline reject. Again, I think this is an interesting paper which clearly makes incremental progress on this challenging problem. Regretfully, my opinion is still that this paper does not have the novelty/impact necessary for a 2024 neurips publication (particularly because the large-offline setting, though certainly not contrived, I would argue is somewhat niche) but nor do I feel strongly that this article should be excluded from these proceedings.

PS, upon reviewing the article, I noticed on line 188 that it’s not a stationary kernel that you’re looking for; stationary kernels which are periodic form counterexamples to your claim. Rather, it is kernels which decrease monotonically wrt some distance. (Just pointing it out for your benefit; no need for us to litigate this minor point and it does not impact my score; it's indeed clear that this projection will usually be efficiently computed in practice).
Thank you for your feedback on the new experiments and the novelty of our work! We will improve the presentation of the manuscript by incorporating the new experimental results and highlight our novelty in the revised version.
Thank you very much for your rebuttal. Given the new experiments, I am willing to increase my score to 7 (accept).

Regarding the novelty issues mentioned by Reviewer SUei, I believe the TuRBO algorithm became universally accepted as the go-to solution when it comes to high-dimensional Bayesian Optimisation. As such, in my opinion, the fact that FocalBO is merged with TuRBO to produce a stronger baseline is not really a problem, given that it severely outperforms other TuRBO-based baselines. It is clear that the main component critical to success is the FocalBO strategy. I also believe that the FocalBO itself is sufficiently different from the previously proposed sparse-GP frameworks.
**The rationality of FocalBO design**

FocalBO has two main components: the focalized GP and the hierarchical acquisition optimization. In the clarification of the regularization term, we have explained that the scaling is used for maintaining similar scale for similar size of search regions, and our emprical evaluation demonstrate the rationality of our loss function design. Below we give some empirical results and the intuition in the design of hierarchical acquisition optimization.

After each BO iteration, we use a simple yet effective heuristic to adjust the optimization depth for the next round, without introducing additional hyperparameters. We encourage more exploitation (increase the depth) when finding good points in the smallest search region, and encourage more exploration (decrease the depth) when good points are found in boarder search regions. On our set of diverse benchmarks - including synthetic and commonly-used test functions, as well as two challenging high-dimensional problems - this simple heuristic worked without any tuning required.

**Comparison with original TuRBO**

Our current experiments with TuRBO are all re-implementations to make sure baselines have the same computational cost. We run the original TuRBO implementation (with exact GP and Thompson sampling) on both robot morphology design and human musculoskeletal system control task (**Figure R4**). We observed that FocalBO outperforms TuRBO on both tasks with smaller computational cost. The reason of TuRBO's poor performance may be that it cannot quickly adapt over the search space when the online evaluation budget is small.

**Typos and unclear statements**

Thanks for pointing out the typos. We will correct them and improve our paper in the new version.
**The novelty of our work**

Here we highlight the novelty of FocalBO in both GP and BO aspects.

**We design the first sparse GP model that improves acquisition function optimization**. Our design of focalized GP enables **better local estimation, aiming at improving acquisition function optimization during Bayesian optimization**. We adopt variational inference to derive focalized GP, which is **capable of performing joint posterior inference and sampling given a set of test points just like exact GP**. Therefore our proposed GP model is automatically compatible with any acquisition function that is used for exact GPs. (Figure 2). 

The locally induced GP (LIGP) in the cited work [2] is designed for regression tasks, which assigns different inducing points to every test point, aiming at improving the point estimation measured by MSE. The inability of joint posterior computation and discontinuous prediction prevents LIGP from easily being incorporated into the BO framework.

**We propose the first scalable BO algorithm that is capable of utilizing large offline dataset for optimization.** Our design of hierarchical acquisition function optimization  with focalized GP **searches promising positions over different scale of search region during one BO iteration, enabling making decision based on both global and local information under restricted computational budget**. 

The overall idea of TuRBO [3] is to restrict the search space to a fixed size during one BO iteration, and adjust the search region size based on the optimization results. It achieves scalability by discarding previous samples and restarting when the search region reduce to certain threshold. Their use of exact GP still faces the scalability issue and would be computationally intensive when handling large offline datasets. As mentioned in the last paragraph of section 4.2, our proposed framework is orthogonal to TuRBO, and find that FocalBO and TuRBO can have a complementary effect in sections 5.2 and 5.3. 

[2] Cole, D. Austin, Ryan B. Christianson, and Robert B. Gramacy. "Locally induced Gaussian processes for large-scale simulation experiments." Statistics and Computing 31.3 (2021): 33.

[3] Eriksson, David, et al. "Scalable global optimization via local Bayesian optimization." Advances in neural information processing systems 32 (2019).
**Theoretical implications of sparse GP approximation**

We agree that the theory section can be fleshed out more, and will do so by (1) providing more precise formulae for the additional regret incurred due to too-sparse approximations, and (2) experimental and numerical evidence for the increased approximation fidelity of Focalized GP. 

For the former, we do not aim to propose any new theorems due to the restrictive nature of common assumptions, lack of guarantees for ELBO-based methods, and difficulty for controlling the effective number of inducing points within the local region. However, we will better motivate our method by providing examples of high KL error under reasonable assumptions, e.g. those used in SVGP-TS[1].

In **Figure R2**, we also empirically measure our claim that Focalized GP can significantly reduce approximation error on the search region. We sampled 8000 training points from 2d GP functions to train focalized GP and SVGP. Over different size of the search region, we compare the KL divergence of the GP posterior prediction over search region between sparse GPs and the exact GP. We observe that the KL divergence between focalized GP and exact GP is consistently smaller than that between SVGP and exact GP, implying tighter approximation to the exact GP over local region.

While a rigorous regret bound is hard to derive, we conduct an empirical study where we directly compare the optimization performance between focalized GP and SVGP when combining with TuRBO. In this way we can eliminate the influence of hierarchical acquisition optimization. The optimization performances are shown in **Figure R3**. We observe that focalized GP outperforms SVGP on both high-dimensional problems, which empirically demonstrates our theoretical implications that Focalized GP contributes to reducing regret. 

[1] Vakili, Sattar, et al. "Scalable Thompson sampling using sparse Gaussian process models." Advances in neural information processing systems 34 (2021): 5631-5643.