Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization
Authors: Zhitong Xu
Haitao Wang
Jeff M. Phillips
Shandian Zhe
Venue: ICLR 2025 Oral (2025)
Keywords: Gaussian Process
Bayesian Optimization
High Dimensional Bayesian Optimization
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We identified a critical failure mode of standard Bayesian Optimization in high-dimensional optimization and proposed a simple yet effective solution

Abstract: A long-standing belief holds that Bayesian Optimization (BO) with standard Gaussian processes (GP) --- referred to as standard BO --- underperforms in high-dimensional optimization problems. While this belief seems plausible, it lacks both robust empirical evidence and theoretical justification. To address this gap, we present a systematic investigation. First, through a comprehensive evaluation across twelve benchmarks, we found that while the popular Square Exponential (SE) kernel often leads to poor performance, using Mat\'ern kernels enables standard BO to consistently achieve top-tier results, frequently surpassing methods specifically designed for high-dimensional optimization. Second, our theoretical analysis reveals that the SE kernel’s failure primarily stems from improper initialization of the length-scale parameters, which are commonly used in practice but can cause gradient vanishing in training. We provide a probabilistic bound to characterize this issue, showing that Mat\'ern kernels are less susceptible and can robustly handle much higher dimensions. Third, we propose a simple robust initialization strategy that dramatically improves the performance of the SE kernel, bringing it close to state-of-the-art methods, without requiring additional priors or regularization. We prove another probabilistic bound that demonstrates how the gradient vanishing issue can be effectively mitigated with our method. Our findings advocate for a re-evaluation of standard BO’s potential in high-dimensional settings.

Published as a conference paper at ICLR 2025
STANDARD GAUSSIAN PROCESS IS ALL YOU NEED FOR HIGH-DIMENSIONES BAYESIAN OPTIMIZATION
Zhitong Xu, Haitao Wang, Jeff M. Phillips, Shandian Zhe Kahlert School of Computing University of Utah, Salt Lake City, UT 84112, USA
(u1502956, haitao. wang) (anng))uuaheedu (jeffp, zhe) www.esschh............e..eeaeeae edu
ABSTRACT
A long-standing belief holds that Bayesian Optimization (BO) with standard Gaus- sian processes (GP) referred to as standard BO underperforms in high- dimensional optimization problems. While this belief seems plausible, it lacks both robust empirical evidence and theoretical justification. To address this gap, we present a systematic investigation. First, through a comprehensive evaluation across twelve benchmarks, we found that while the popular Square Exponential (SE) kernel often leads to poor performance, using Matérn kernels enables stan- dard BO to consistently achieve top-tier results, frequently surpassing methods specifically designed for high-dimensional optimization. Second, our theoreti- cal analysis reveals that the SE kernel's failure primarily stems from improper initialization of the length-scale parameters, which are commonly used in prac- tice but can cause gradient vanishing in training. We provide probabilistic bound to characterize this issue, showing that Matérn kernels are less suscep- tible and can robustly handle much higher dimensions. Third, we propose a simple robust initialization strategy that dramatically improves the performance of the SE kernel, bringing it close mate-offthe-eee methods, without requiring additional priors regularization. We prove another probabilistic bound that demonstrates how the gradient vanishing issue can be effectively mitigated with our method. Our findings advocate for re-evaluation of standard BO's potential in high-dimensional settings. The code is released at atttpps:/ github. com/ XZTC08/ Standa rd- GP- is -you- need foreed-for-tae Foo-HDRO
1 Introduction
Many applications require optimizing complex functions, allowing queries only for the function val- ues, possibly with noises, yet without any gradient information. Bayesian Optimization (BO) (Snoek al. 2012 Mockus 2012) popular approach for addressing such challenges. BO typically employs Gaussian process (GP) (Williams & Rasmussen] 2006) as a probabilistic surrogate. iteratively approximates the target function, integrates the posterior information to maximize an acquisition function for generating new inputs at which to query, then updates the GP model with new examples, and concurrently approach the optimum.
Despite numerous successes, there has been widespread belief that BO with standard GP regression, referred to as standard BO, is limited
diw-dimnnioonll optimization problems (Frazier 2018
Nayebi et al. 2019; Briksson & Jankowiak] [021 Moriconi al.) 2020: Letham et al. 2020 Wang et al. 2016
et al. 2016). is commonly thought that the number of optimization variables should not exceed 15 or 20, beyond this threshold, BO is prone to failure (Frazier 2018 Nayebi et al.) 2019). A continuous line of research is dedicated to developing novel high high-dimensional BÔ methods. The fundamental strategy is to impose strong structural assumptions into surrogate modeling avoid directly dealing with high-dimension structure within the functional space (Kandasamy et al. 2013 Rolland et al. [018]] Han et al. 202T Ziomek & Ammar 2023), where the target function is expressed as the summation of group of low-dimensional functions, each operating over small number of variables. Another family of methods assumes that the inputs are intrinsically low-rank. These methods either project the original
to inputs. One class methods assumes a decomposition
"Corresponding author

---

Published as a conference paper at ICLR 2025
inputs into a low-dimensional space (Wang al. 2016 Nayebi al. 2019: Letham et al. 2020 or use sparse-inducing priors to trim down massive input variables (Erksson & Jankowiak, 2021). Subsequently, GP surrogates are built with the reduced input dimensions.
While the aforementioned concerns may seem plausible, there is a lack of both strong empirical evidence and theoretical justification to confirm and explain why standard BO would be ineffective in high-dimensional optimization. To bridge this gap, we systematically investigated standard BO in this paper. The major contributions of our work are summarized as follows:
Empirical Results. We investigated BO with standard GP across eleven widely used benchmarks and novel benchmark, encompassing six synthetic and six real-world digendimensional optimization tasks. The number of variables ranged from 30 to 1,003, We compared standard BO with nine state-of-the-art nigh-dimensional BO methods, and performed extensive evaluations. Surprisingly, while the popular SE kernel often led to poor performance, switching to the ARD Matérn kernel enabled standard BO to nearly always achieve the best (or near-best) optimization performance. When successful, standard BO seems more flexible in accommodating various structures within the target functions. Theory. Through analyzing the gradient structure of the GP likelihood, we identified the primary failure mode of standard BO, in particular when using the SE kernel. Specifically, the commonly used initializations for the length-scale parameters, such setting them to one, is improper in high-dimensional settings and can easily cause gradient vanishing. preventing effective training. We proved a probabilistic tail bound characterize this issue under mild conditions, Applying the same analytical framework to the Matérn kernel, we show through comparison of the tail bounds that the Matérn kernel is less prone to gradient vanishing with the same length-scale initialization, making it more effective for handling high-dimensiones problems.
Simple Robust Initialization. Based on our theoretical results, we proposed a simple yet robust length-scale initialization method, without requiring any additional priors or regularization. We proved a probabilistic bound showing that the probability of gradient vanishing decreases exponentially with the increase dimensionality, implying that the gradient vanishing issue can be effectively mitigated in high- nigh-dimensional settings. Em- pirical evaluations demonstrate that our initialization method dramatically improves the performance of standard BO with SE kernels, enabling it to achieve, come close to, the satte-of-the-ass art performance in high-dimensional optimization.
2 Standard Bayesian Optimization
Consider maximizing a d-dimensional objective function f C Rd + R, where the function form is unknown. We query the function value any input X, possibly with noises, but no gradient information is accessible. We aim to find xt argmaxxex//xxx To achieve this, Bayesian Optimization (BO) employs probabilistic surrogate model to predict across the domain X, while also quantifying the uncertainty of the prediction. This information is then integrated to compute an acquisition function, which measures the utility of querying at any new input location given the current function estimate. By maximizing the acquisition function, BO identifies a new input location which to query. ideally closer to the optimum. Concurrently, the acquired new example is incorporated into the training dataset, and the surrogate model is retrained to improve the accuracy. BO begins by querying at a few (typically randomly selected) input locations, and trains an initial surrogate model. The iterative procedure repeats until convergence or stopping criterion is met.
The standard BO adopts GP regression (Williams & Rasmussen 2006) for surrogate modeling. GP prior is placed over the target function, ~~ GP (m(xx, <(x, )), where m(x) is the mean function, which is often set as a constant, and s((--) -) is the covariance function, which is usually chosen as a Mercer kernel function. One most popular kernel used in BO is the Square Exponential (SE) kernel with Automatic Relevance Determination (ARD),
KisE x') = a (exp(=p²),
where a > 0 is the amplitude, V(x x') diag ( (x- x'), and l ea) 0 are the length-scale parameters. We refer this kernel as ARD because each input dimension has
2

---

Published as a conference paper at ICLR 2025
distinct length- scale parameter. An alternative choice is the ARD Matérn-5// 5/2 kernel,
KMatém X, x') - V5p + 5/2 exp (-v5p).
(2)
Given training inputs X = X1, ,,,...,,NN] xN|" and (noisy) outputs y = [e/...... yN VN let us denote the function values at X as f [f(x)] ..... ((x)) multi-variate Gaussian distribution, p(f) = N (f|m, K), where m = [m((x))..., (xx)]" K is
which according to the GP prior. follow
the kernel matrix on X, and each [K]; = (x,,,, Xy). Then one can employ a Gaussian likelihood for the observations and
x (y K, o²I),
where o² is the noise variance. To estimate the GP parameters, e.g., the length-scales L, one can maximize the marginal likelihood (3). The predictive distribution at any new input x", is condi- tional Gaussian, p ((x))ly] = (xx()) u(x")) )) where u(x*) = m(xx) + =(x", XK(K
oy m). and v(x") =(x* x²) <(x", (K + I²1) *x(x x*).
In each iteration, given the current GP surrogate model, an acquisition function is maximized to identify the next input location at which to query. One commonly used acquisition function is the upper confidence bound (UCB), defined as UCB(x) represents the exploration level. There have been other popular acquisition functions, such as Expected Improvement (EI) (Jones et al. 1998). Thompson sampling (TS) (Russo al. 2018), the
μ(x) + Av/u(x),
E X, where
recently proposed log-E! EI (Ament et al. 2024), among others.
3 High Dimensional Bayesian Optimization
An enduring and widespread belief is that when the input dimension d is high, e.g., a few hundreds, the standard BO prone to failure. This belief might partly arise from an intuition that commonly used kernels, such as (I). could encounter challenges in dealing with high- sigh--imensional inputs, making GP struggle in capturing the target function. A dedicated line of research is developing novel high- dimensional BO methods. The key idea of these methods is to impose some structural assumptions in surrogate modeling to sidestep directly handling the high-dimensional inputs in kernels.
Structural Assumption in Functional Space. The first class of methods assumes a decomposition structure within the functional space. Specifically, the target function f is modeled as
f(x) = Σ fx(x),, fo ~ GP (m,(x'), (z)))))))
where each x' C x is a small group of input variables, and x = x U U U There can be variety of choices for the group number M and each variable group x². In (Kandasamy al., 2015), the input is partitioned into non-overlapped groups After every few BÒ iterations, new partition selected from set of random partitions. The selection based on the model evidence. the maximum
In (Rolland et al. 2018), the variable groups can overlap, and are represented cliques on a dependency graph. A Gibbs sampling method was developed to learn the structure of the dependency graph. (Han et al. (2021) proposed Tree-UCB, which restricts the dependency graph to be tree-structured so as to boost the efficiency of structure learning. However, the more recent work (Ziomek & Ammar 2023) points out that, learning the group structure through relatively small data can be misleading. A wrongly learned structure can cause queries to be stuck at local optima. This work proposes RDUCB that randomly decomposes into a tree structure, and surprisingly works better than those learned structures through various methods.
Structural Assumption in Input Space. Another class of methods assumes a low-rank structure in the input space. Many of them introduce Low-dimensional Embeddings. [Wang et al. (2016) proposed a method named REMBO, building the GP surrogate in a low-dimensional embedding space, 2 C Re where demb < d. The acquisition function is optimized in the embedding space. When querying the function value, the input is recovered by x Az where z Z, and A is a random matrix. If the recovered x is out of the domain X. it will be clipped boundary point. To avoid clipping, Nayebi et al. (209)) proposed HESBO, which randomly duplicates z (and/or flip the signs) to form x. This way. the optimization of can respect the same bounding constraints in X. This essentially imposes further constraints on the embedding space. A more flexible solution, ALEBO, was proposed in (Letham et al. 2020). ALEBO also random matrix A recover the
3

---

Published as a conference paper at ICLR 2025
inputx: But when maximizing the acquisition function, ALEBO incorporates a constraint that Az must be in X. thereby avoiding the clipping issue. Another direction is to impose Partial Variable Dependency by triming down the imput variables. The recent start-of-the-an SaasBO (Eriksson & Jankowiak 2021), uses horse-shoe prior (Carvalho et al. 2009) over the length-seale parameter for each variable. As the horse-shoe prior induces strong sparsity, massive number of variables can be pruned, substantially reducing the input dimension.
Theoretical Analysis
To understand why standard BO can fail in high-dimensional cases, in particular with the SE kernel (see Section 6), we look into the gradient structure during the GP training. Specifically, we consider at the beginning, all the length-scale parameters are set to the same initial value, E =... ld = lg. From the marginal likelihood
(3) we derive the gradient w.r.t each length-scale,
2 alg
OK
tr
A
O2k
Olx
where = aa T (K + - and a = (K + I²1) -0 (y m). When using the SE kernel as
specified in (T.. each entry of @K as follows: 86
@KK
and for j.
88 ii
OK 06
ÓNSE Xi, x1)
2a
||x x,,l
(Nin xzx)²
2a
0²
exp(
(6)
Ck
c
where Vin and Zizk are the k-th element in x, and Xj, respectively, and
to
Intuitively, as the input dimension d increases, the squared distance between input vectors grows rapidly. Consequently, the factor - in (6) can quickly fall below the machine epsilon &, causing the gradient vanishing issue (Bengio et al. 1994 Hanin 2018). That is, the length-scale cannot be effectively updated according to ÖesE x,, x due to the limitations of the numerical precision.
o² Proposition 4.1. Given any & > 0, Gor < when P TSE 1
+ log
When using the double-precision floating-point format (float64), the rounding machine epsilon is & 2 2-33 and the above threshold is TsE 6.58. To analyze the likelihood of gradient vanishing, provide a probabilistic bound under uniform data distribution.
Lemma 4.2. Suppose the input domain is [0, 1]d and the input vectors are sampled independently from the uniform distribution, then for any constant threshold O, when d >> 66272
p(p 2 1 (d 2) (d
(7)
18d
The lower bound grows exponentially with the input dimension d and converges to one. This implies that given any fixed choice of to relevant to (e.g., fo = 1), as d increases, the probability p(p 2 TSE) will rapidly approach one, leading vanishing gradient for every (sE,(,, (XxXxx,,, Xj). As consequence, each length scale la in 5 cannot be effectively updated, preventing meaningful training The inferior prediction accuracy of the GP model further contributes to poor performance the BO procedure. Since TSE is small, gradient vanishing can occur before d becomes extremely large. For example, for & = 0.5 and to = 1, the probability p(p > TSF) exceeds 0.99 when d' > 205 and 2 473, respectively.
Why Matérn kernels perform much better? In our practical evaluations (see Section 6), standard BO with Matérn kernels @ often perform much better than with the SE kernel in high-dimensional problems, frequently achieving state-of-the-an results. To understand the reasons behind this im- provement, we apply the same framework analyze the behavior of GP training with the Matérn
'A union bound can be directly derived from 7 to lower bound the joint probability that all gradients (Ōrsa (xi, / Ji, across N training instances are below machine epsilon. However, this bound is too Joose to reflect our practical observation (see Fig. and and Table Therefore, we focus our discussion on single pair of inputs.
4

---

Published as a conference paper at ICLR 2025
Lower bound
0.95
0.99
0.995
0.999
0.9995
0.9999
SE (d)
172
205
219
250
264
294
Matérn (d)
980
1040
1064
1116
1137
1185
Table I: The input dimension d vs. lower bound probability of gradient vanishing under &u 0.5.
kernel. First, we obtain the gradient,
dKMaiem Xi, x,)
3 (1 v5p) v5p) (rik - =m²k)² 5(1 551+ P²
3
ev5p
As the distance increases, the factor p² /evop will also converge zero, implying that with the growth of d, the gradient will eventually fall below machine epsilon. However, since this factor p² / V5p decreases much slower than the corresponding factor to the SE kerne,, namely p² /ee in Matérn kernel able to robustly handle much higher dimensions.
66 the
Proposition 4.3. Given VA > O,
_ ev3p
< & when p > TMatérm = (1 VI + log 1/5 log { /v5.
With the machine epsilon & = 2- 2-53, the threshold TMatem = 21.98. more than 300% larger than TSE = 6.58. Therefore, to achieve the same probability bound in 77. a much higher d is required for Matérn kernel compared to the SE kernel under the same fo. For concrete examples, we list the dimension d versus the probability bound in Table [ implies that using the same length-scale initialization, the Matérn kernel is able to robustly handle higher dimensions.
Numerical Verification. To verify whether gradient vanishing occurs during practical GP training and to access whether our analysis aligns with the empirical observations, we conducted numerical experiments using two synthetic functions, Hartmann6(d,, 6) and Rosenbrock(d.d) (see Section for definitions). We varied the input dimension d across (50, 100, 200, 300, 400, 500, 600), with the domain set to [0, 1]d_ For each value of d, we uniformly sampled 500 training inputs and 100 test inputs. We trained GP with both SE and Matérn kernels, using length-scale initialization to from {0.1, 0.5, 0.693, 1.0, vd). Note la 0.693 SoftPlus(O) is popular initialization choice, used as the default choice in influential GP/BO libraries GPyTorch (Gardner et al. 2018) and BOTorch (Balandat al. 2020). Training and testing were repeated 20 times for each and to. We then evaluated (A) the average relative Lg difference between the length-scale vector before and after training, namely, |Iaindd Lmaill fmmll//|[eee e|iill (B) the average gradient norm at the first training step, and (MSE).
(C) the average test Mean-Squared-FFree
As shown in Fig. D gradient vanishing indeed occurred, leading to training failure. When using to = 0.1 for GP training with the SE and Matérn kernels, both the length-scale relative L2 difference and the gradient norm remained close to zero across all choices of d (those values are not displayed because they are too small). Consequently, the MSE hovered around 1.0, indicating consistent training failure. When we increased ln to 0.5 and 1.0 for training with the SE kernel, both the length-scale difference and gradient norm became significantly larger and non non-trivial for d = 50 and d 100, resulting in MSE below 0.2. This demonstrates that gradient vanishing did not occur and training was successful. However, as d increased, gradient vanishing reappeared at d 200 for & 0.5 and at d = 400 for to = 1.0, for both Hartmann6 and Rosenbrock. In these cases, the length-scale difference and gradient norm dropped near zero, while the MSE returned to around 1.0. A similar trend was observed with the Matérn kernel. However, GP training with the Matérn kernel can succeed with higher dimensions under the same initialization conditions. For example, with lu 0.5, for both Hartmannó and Rosenbrock, training with the Matérn kernel started to failed until d = 400, compared to d = 200 with the SE kernel. When using t 1.0, the Matém kernel succeeded across all d values, while the SE kernel started to fail at d = 400. Table [ summarizes the input dimensions at which gradient vanishing and training failure began to happen.
Confirmation in BO Running. We next investigated whether gradient vanishing occurs during the BO process, potentially leading to poor optimization performance. To this end, we examined two benchmark BO tasks: SVM (D=388) and Rosenbrock (300,100), as detailed in Section 6 As shown in Fig. 3 and 4] standard BO with the Matérn kernel (SBO-Matérn) achieved top-tier results, while BO with the SE kernel (SBO-SE) performed poorly. Both methods used the initialization & 0.693. We then extracted the BO trajectories for both methods and analyzed the relative L2 difference in the length-scale vectors before and after training, and the gradient norm at the first training step.
5

---

Published as a conference paper at ICLR 2025
As illustrated in Fig. 2 during the BO iterations, GP-SE (f, 0.693) consistently experienced
gradient vanishing and so training failure, whereas GP-Matérn (&) = 0.693) maintained non-trivial
gradients. This difference directly correlates with their respective optimization performances poor for SBO-SE and excellent for SBO-Matérn.
0.1
0.5
0.693
1.0
sqrt D
Hartmann6 (Matérn)
Hartmann6 (SE)
Rosenbrock (Matém)
Rosenbrock (SE)
1.0
1.0
1.0
1.0
0.8
MSE
0.2
0.2
0.2
0.2
50100
200
390
400
500
600
50,00
200
300
400
500
500
54000
200
200
409
500
900
50300
200
300
400
500
600
10a
102
10a
10a
Relative Length Scale Diff
109
10
100
109
10-2
10-2
10
10-2
10-3
10-3
10
10-2
S0ide
200
300
409
500
600
50.00
200
300
400
500
600
50.00
200
300
400
500 600
50000
200
300
400
500
600
10
10
10
10
Grad Norm
10-6
10-6
10-6.
10
10-8
10-8
10-8.
10-1
10-12
10 -12
10-12
10-12
Figure 1: Mean Square Error (MSE), Relative L2 difference between the length-scale vectors before and after training, and L2 norm of the length-scale gradient at the first training step, across different input dimensions (r-axis) and initializations (legend). The results were averaged from 20 runs. For to 0.1, the length-scale relative difference and gradient are displayed, they fall much below the minimum values the y-axis.
GP SE(L, 0.693)
GP Matérn 0.693)
GP SE vd)
SVM
Rosenbrock (300.100)
SVM
Rosenbrock (300,100)
Relative Length-San Scale DIf
10
101-
10-3
30-1.
10-1
Grad Norm
30
20-3
10-1
10-13.
30-5
10-5
10-4
20-7
10-1
10-11
10-9
10-1
Nurbe"
Nimber ot Itmmioons
Nimber Encyinee
Numitor Femitions
Figure 2: The relative La differencc of the length-scan scale vectors before and after training and the gradient norm at the first training iteration across every step of BO. The results for GP-SE (l, 0.693) are not displayed, as they fall much below the minimum values the y-axis.
Robust Initialization
From the analysis in Section 4 one straightforward approach to improving learning with the SE (and also Matérn) kernel is initialize each length-scale with larger fo. which loosen the lower

---

Published as a conference paper at ICLR 2025
Hartmann6
to 0.1
0.5
0.6931
1.0
vd
SE:
50
200
200
400
v
Matérn
50
400
600
v
v
Rosenbrock
tu 0.1
0.5
0.6931
1.0
Vd
SE
50
200
200
400
Matérn
50
400
600
Table 2:
The input dimension from which gradient vanishing started to occur and training started to fail; means gradient vanishing occurred.
bound in 7 and allow the GP model better adapt higher dimensions. However, this approach is not robust enough, as d - 6f@r² will continue to increase with d, and the bound will still rise toward one rapidly, causing the gradient to vanish. To address this issue, we propose setting to growing with d rather than keeping it constant,
fo =vddd c>0.
Lemma 5.1. Suppose the input domain is [0, 17th and each input vector is independently sampled from the uniform distribution. Given any constant threshold T > O, we set to = cvd. such that /or then
p(p 2 T) 5 2exp -2(62--
(9d) x exp(-O(d)).
With this initialization, the increase exponentially reduces the upper bound on the probability of exceeding the given threshold As a result, the gradient vanishing issue can be fundamentally mitigated the high-dimensional settings.
In our experiments, tested standard BO using SE kernel with our proposed length-scale initializa- tion as specified in (8) where we set c= 1. As shown in Fig. and Fig. 4 the performance across all the benchmarks is dramatically improved, as compared to using popular initialization & 0.693. The performance matches that of using Matérn kernels, achieving state-of-the-ar (or near-best) results. In addition, shown in Fig- aand Table our numerical experiments confirm that gradient vanishing never occurred with our robust initialization method. Fig. 2 further demonstrates that our method eliminates gradient vanishing during practical BO runs, enabling effective GP training and much better optimization results. Together this shows that our method can effectively adapt SE kernels high-dimensional settings, and mitigate the failure mode of standard BO.
Alternative method. Recent work (Hvarfner et al. 2024) proposes an intuitive approach to address the learning challenge of the SE kernel the high-dimensional spaces by constructing a log-norma prior over each length scale, p((x) = LoggNrmaal(4Aler + log(d) oo). This prior regularizes the length- scale to be the order of Vd during training. Using theoretical framework, this method can be justified as an alternative strategy to alleviate the gradient vanishing issue. However, our method, directly motivated by analysis, is simpler- requiring additional prior construction or associated parameters, and has also shown effective in addressing the gradient vanishing issue. The comparison with (Hvarfner et al. 2024) is given in Section 6
2
Comprehensive Evaluation
6.1 Experimental Settings
For a comprehensive evaluation the standard BO, we employed twelve benchmarks, of which six synthetic and six are real-world problems. For clarity, all the tasks aim for function maximization.
Synthetic Benchmarks. We considered four popular synthetic functions: Ackley, Rosenbrock, Hartmann6, and Stybtang. Definitions of these functions provided in Section B. of the Appendix. Each task is represented by "Fun(d, d')", where is the input dimension that BO methods optimize for, and d' is the number of effective variables (d' < d). The ground-truth of the target function
7

---

Published as a conference paper at ICLR 2025
values is computed with the first d' variables. The benchmarks and the structures within the target functions are summarized in Appendix 3
Real-World Benchmarks. We employed the following real- world benchmark problems: Mopta08 (124) (Jones 2008), a vehicle design problem. The objective is to minimize the mass of vehicle with respect to 124 design variables, subject to 68 performance constraints. We followed (Eriksson & Jankowiak] 2021) to encode the constraints as a soft penalty, which is added into the objective. SVM (388) (Eriksson & Jankowiak] 202Ip. hyper-parameter tuning problem that optimizes regularization parameters and 385 length-scale parameters for a kernel support vector machine regression model. Rover (60), a rover trajectory planning problem from (Wang et al. 2018). The goal is to find an optimal trajectory determined by the location of 30 waypoints in a 2D environment, DNA (180) (vSehic et al. 2022). a
optimization problem that optimizes 180 regularization parameters for weighted lasso regression on DNA dataset (Mills) 2020). Prior analysis (vSehic al., 2022) shows that only around 43 variables are relevant to the target function. NAS201 (30) Dong & Yang 2020), neural architecture search problem CIFAR-100 dataset. Humanoid Standup (1003): A novel trajectory optimization benchmark based on a humanoid simulator that uses the MuJoCo physics engine (odorov et al. 2012). The problem dimension is 1,003. The details are provided in Appendix Section B.2
Methods. We tested with four versions of standard BO:
SBO-Matém: GP regression with ARD Matérn-572 kernel as specified in (2) the length- scale initialization Eo = SoftPlus(O)- = 0.693, the default choice in the influential GP/BO libraries GPyTorch and BOTorch.
SBO-SE: GP regression with ARD SE kernel as specified in ① the length-scale initialization fn = (oftPlus(O) = 0.693.
SBO-Matérn (RI): GP regression with ARD Matémn-522 kernel, using our proposed robust initialization to = cvd d where d is the input dimension and
SBO-SE (RI): GP regression with ARD SE kernel, using our robust initialization with 1.
We used UCB as the acquisition function where the exploration level was set to 1.5. We used GPyTorch for GP training and BOTorch for Bayesian optimization. To ensure efficiency, the GP was trained via point estimation of the length-scale and noise variance parameters, with the optimizer selected from L-BFGS, Adam, or RMSProp. Further details are provided in Appendix Section B.3
We compared with six state-of of-the-art dighediimessioall BO methods, including Tree UCB (Han et al. 2021), RDUCB (Ziomek& & Ammar 2023). HESBO (Nayebi et al.] 22119) ALEBO (Letham et al. 2020), SaasBO (Eriksson & Jankowiak 2021). and TURBO (Eriksson et al. 2019). The first five have been introduced in Section TURBO is designed for special scenario where the target function can be extensively evaluated, e.g., tens of thousands evaluations. This is uncommon in BO applications, as the function evaluation is typically deemed expensive, and one aims to evaluate the target function as few as possible. TURBO searches for trust-regions in the input domain, learns a set of local GP models over each region, and then uses implicit bandit strategy to decide which local optimization runs to continue. In addition, we tested the most recent method VBO (Vanilla BO) (Hvarfner et al.) (022)) that levels up the performance of the standard BO by constructing a log-normal prior over the length-scale parameters (see Section We used the original implementation and default (recommended) settings for the competing methods. Besides, we compared with two BO methods based on Bayesian neural networks (BNN). One uses Laplace approximation and the other Hamiltonian Monte-Carie (HMC) sampling for posterior inference. We used high-quality open- implementation for each BNN based method. Appendix A the structural assumptions made by all the methods. The details of all the competing methods are given in Appendix Scction B.3
For each optimization task, we randomly queried the target function 20 times to collect the initial data, except for Humanoid-Standug where we collected 50 initial data points. We tested TURBO with a single trust region and five trust regions, denoted TURBO-1 and TURBO-5 respectively. For HESBO and ALEBO, we varied demb E {10. 20) and denote the choice as HESBO/ALEBO- (10, 20). For SaasBO, we performed tests using both NUTS and MAP estimation for surrogate training, denoted by SaasBO-NUTS and SaasBO-MAP respectively. For a reliable comparison, we ran each method ten times, ensuring that all methods used the same randomly collected initial dataset in each run (these datasets varied across different runs). Additionally, we conducted ten extra runs for all the standard BO methods and VBO further validate their performance.
8

---

Published as a conference paper at ICLR 2025
SBD Matern
580 SE
TURBO-5
HESBO 20
ALEBO-20 20
SaasBO NUTS
RDUCB
5B0-Matern (RI)
(R)
TURB0-1
HESBO-10
ALFBO-10
SaasBO-MAP
V80
Rosenbrock (000,
Ackley (150,150)
StybTang (200,200)
40003
-20002
-40038
10000
-14
-40000
330000
59302
2226038
-70000
-246028
Rasenbrock (300,100)
Ackley (300,150)
Hartmann6 (300,
-20038
40000
2.5
-40030
.3000
330000
320030
-340038
3.5
Figure 3: Optimization performance in synthetic benchmarks.
SBD Matern
S80 SE
TURBD-5
HESBO-20 20
ALEBO-20
SaasBO NUTS
RDUCH
Matern
(RI)
ALEBO-10
UCB
VB0
Mopta08 (124)
Rover (60)
NAS201 (30)
-220
6:12
Rnal value
-240
DNA (180)
SVM (388)
Humanoid-Stande (1003)
-6.3 -75
-0.31
120
-4.34
46
Figure 4: Optimization performance in real-world problems.

---

Published as a conference paper at ICLR 2025
6.2 Optimization Performance
Ultimate Outcomes. We first examined the ultimate optimization outcomes of each method. The results are reported in Fig. 3 for synthetic benchmarks and in Fig. for real-world benchmarks. Due to space limit, supplement the BNN-based BO results in Appendix Section C.T We used violin plots that illustrate the average of the obtained maximum function value over the runs (the middle bar), the best (the top bar), the worst (the bottom bar), and the distribution of the obtained maximum (the shaded region).
average, SBO-Matérn (RI), SBO-SE (RI) and VBO consistently deliver top-tier results, often outperforming the other methods by a large margin, e.g., on benchmarks Ackley (300, 150), DNA, and SVM. It is noteworthy that 3BO-Matém, despite using a small length-scale initialization (4, 0.693), also achieved top-tier performance across all benchmarks except for Hlumanidd-Standup where the optimization dimension exceeds 1,000. In contrast, SBO-SE failed when the input dimension due to the gradient vanishing issue, e.g., Ackley(150, 150), Rossenbrock(300, 100), DNA, and SVM (see our thcoretical analysis and numerical experiments in Section 4p. The result verifies that the Matérn kernel is far less susceptible to gradient vanishing and performs much more robustly for high-dimensional problems than the SE kernel. However, when the dimension is extremely high (e.g., over IK) the small length-scale initialization can still cause vanishing gradients and degraded optimization performance. After applying our robust initialization, both kernels achieved top-tier performance consistently, closely matching VBO.
150
Runtime Performance. Next, we examined the runtime performance. Appendix Fig. 5 and 6 the average maximum function value obtained at each step and the standard deviation across all runs. Since we have 16 methods for comparison, showing the standard deviation for every method makes the figure cluttered. To ensure clarify, we show the curves and standard deviation of SBO, VBO, and the best, worst, and median performed remaining baselines so far at each step.
report
After an initial stage, SBO-Matérn (RI), SBO-SE (RI), and VBO consistently produced superior queries, achieving larger function values. So did SBO-Matérn in all the cases except Humanoid- Standup. Their performance curves are generally above those of competing methods, reflecting a more efficient optimization of the target function. This trend is particularly noticeable in benchmarks such as Ackley (300, 150), Rosenbrock (300, 100), and SVM. The performance of SBO-SE degraded when d 2 150, with its curve positioned toward the bottom due the gradient vanishing issue. In benchmarks such as Rosenbrock (100, 100) and DNA (180), the curve for VBO is predominantly above those of SBO-Matérn (RI) and SBO-SE (RI), indicating even more efficient optimization. However, on Ackley (150, 150), Ackley (300. 150) and Humanoid-SSanduup (1003), VBO's curve is mostly below that of SBO-Matérn (RI) and SBO-SE (R/), and exhibits worse performance. For the remaining benchmarks. their curves are close and often overlap. Overall, in all cases, they tend to converge to similar function values.
These results collectively demonstrate that standard BO can excel in high high-ddmeesiionl optimization problems. The lack of additional structural assumptions GP modeling may make more flexible in capturing various correlations within nigh-dimensional spaces. The primary challenge for standard BO arises from the gradient vanishing issue, which is often due to improper initialization of the length-scales. By using our robust initialization method, this risk can be significantly mitigated.
Additional Results. We conducted extensive additional evaluations to examine exploration parame- ters in UCB, and other acquisition functions SBO performance in high-dimensional settings. We found that UCB notably robust to variations in the exploration parameter. Additionally, log-El facilitates strong performance in high-dimensional optimization, whereas El and Thompson Sampling (TS) often lead to subpar results. Detailed results and discussions are provided in Appendix Section
7 Conclusion
We conducted thorough investigation of standard BO in high-dimensional optimization, both empirically and theoretically. Our analysis identified gradient vanishing as a major failure mode of standard BO, and we proposed simple yet robust initialization method to address this issue. Our empirical evaluation demonstrates that, once this problem is mitigated, standard BO can achieve oate-of-the-art performance in high-dimensional optimization.
10

---

Published as a conference paper at ICLR 2025
Acknowledgement
This work has been supported by NSF CAREER Award IIS 204490, NSF OAC-2311685, and Margolis Foundation.
References
Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Unex- pected improvements to expected improvement for bayesian optimization. Advances in Neural Information Processing Systems, 36, 2024.
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wil- son, and Eytan Bakshy. Botorch: A framework for efficient monte-carlo bayesian optimization. Advances in neural information processing systems, 33:21524 -25388, 2020.
Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE transactions on neural networks, 5(2):157-16,, 1994.
Carlos M Carvalho, Nicholas G Polson, and James G Scott. Handling sparsity via the horseshoe. In Artificial intelligence and statistics, pp. 73-80. PMLR, 2009.
Xuanyi Dong and Yi Yang. Nas-bench-201: Extending the scope of reproducible neural architecture search. arXiv preprint DeXiv:2001.003220 2020.
David Eriksson and Martin Jankowiak. High- dimensional bayesian optimization with sparse axis- aligned subspaces. In Uncertainty in Artificial Intelligence, pp. 493-503, PMLR, 2021.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan Turner, and Matthias Poloczek. Scalable global optimization via local bayesian optimization. Advances in neural information processing systems, 32, 2019.
Peter Frazier. A tutorial bayesian optimization. arXiv preprint arXiv: 000.:00077001111 2018.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrin gaussian process inference with gpu acceleration. Advances in neural information processing systems, 31, 2018.
Eric Han, Ishank Arora, and Jonathan Scarlett. High-dimensional bayesian optimization via tree- structured additive models. In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pp. 7630-7638, 2021.
Boris Hanin. Which neural net architectures give rise exploding and vanishing gradients? Advances in neural information processing systems, 31, 2018.
Carl Hvarfner, Erik Orm Hellsten, and Luigi Nardi. Vanilla bayesian optimization performs great in high dimension. arXiv preprint arXiv:2402.02222, 2024.
Donald R Jones. Large-scale multi-disciplinary mass optimization in the auto industry. In MOPTA 2008 Conference (20 August 2008), volume 64, 2008
Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13:455-492, 1998. hetLh:ttss//// /api
Remanticschol. crrcchoola... org/ Cowpus/rpuuusaaaa ID: 2638644111
Kirthevasan Kandasamy, Jeff Schneider, and Barnabás Póczos. High dimensional bayesian opti- misation and bandits via additive models. In International conference on machine learning, pp. 295-304. PMLR, 2015.
Ben Letham, Roberto Calandra, Akshara Rai, and Eytan Bakshy. Re-examining linear embeddings for high-ding dineh-iimnnionall bayesian optimization. Advances in neural information processing systems, 33:15 13.15615558. 2020.
11

---

Published as a conference paper at ICLR 2025
Chun-Liang Li, Kirthevasan Kandasamy, Barnabas Poczos, and Jeff Schneider. High dimensional bayesian optimization via restricted projection pursuit models. In Arthur Gretton and Chris- tian C. Robert (eds.). Proceedings of the 19th International Conference and Statistics, volume 51 of Proceedings of Machine Learning Research, pp. 8924892,-00 Cadiz. Spain, 09-11 May 2016. PMLR. httl:stts:/// ptt.ttpsss//pooeeeiini wttpgg/.........//55551111151111 html
Artificial Intelligence
ml ptt.s//51/11111111
Peter Mills. Accelerating kernel classifiers through borders mapping. Journal of Real-Time Image Processing, 17(2):313-327, 2020.
Jonas Mockus. Bayesian approach to global optimization: theory and applications, volume 37.
Springer Science & Business Media, 2012.
Riccardo Moriconi, Marc Peter Deisenroth, and KS Sesh Kumar. High-dimensional bayesian
optimization using low-dimensional feature spaces. Machine Learning, 109:1925-1943, 2020.
Amin Nayebi, Alexander Munteanu, and Matthias Poloczek. A framework for Bayesian optimization in embedded subspaces. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 4752 4761. PMLR, 09-15 Jun 2019. URL httpp:://pooceccinnnng m1 wwwppppsssss..ee wtwp////////iiii/// 19m.html
Paul Rolland, Jonathan Scarlett, Ilija Bogunovic, and Volkan Cevher. High-dimensional bayesian optimization via additive models with overlapping groups. In International conference on artificial intelligence and statistics, pp. 298-307, PMLR, 2018.
Daniel Russo, Benjamin Van Roy, Abbas Kazerouni, Ian Osband, Zheng Wen, al. A tutorial on Thompson sampling. Foundations and Trends@ in Machine Learning, 11(1):1-96, 2018.
Kenan vSehié, Alexandre Gramfort, Joseph Salmon, and Luigi Nardi. Lassobench: A high- dimensional hyperparameter optimization benchmark suite for lasso. In International Conference on Automated Machine Learning, pp. 2-1. PMLR, 2022.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. In Advances in neural information processing systems, pp. 2951-2959, 2012.
yl'ya Meerovich Sobol'. On the distribution of points in cube and the approximate evaluation of integrals. Zhurnal Vychislitel 'nii Matematiki Matematicheskoi Fiziki, 7(4):784-800 802, 1967.
Surjanovic and D. Bingham. Virtual library of simulation experiments: Test functions and datasets. Retrieved January 23, 2024, from homm/ttp:://
Emanuel Todorov, Tom Erez, and Yuval Tassa. Mujoco: A physics engine for model-based control. In 2012 TEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 5026-5033, 2012. doi: 10. 109/ROOS SAT.R/S...12.668886 2012.6386199.
Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale bayesian optimization in high- nigh-dimensional spaces. In International Conference on Artificial Intelligence and Statistics, pp. 745-754. PMLR, 2018.
Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando De Feitas. Bayesian optimization in billion dimensions via random embeddings. Journal of Artificial Intelligence Research, 55:361-387, 2016.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Juliusz Krzysztof Ziomek and Haitham Bou Ammar. Are random decompositions all we need in high dimensional bayesian optimisation? In International Conference on Machine Learning, pp. 43347-43368. PMIRR, 2023.
12

---

Published as a conference paper at ICLR 2025
A Proofs
A.1 Proof of Proposition 4.1
Proof. We intend to solve p²/e" <,, which is equivalent to solving
logr lo
where r = p². We then leverage the fact that log < /r when r > Since 1/e > E. we know the solution of must be larger than 1, implying that the solution for is larger than 1. To meet the inequality, we can solve
logr - vF logg.
It then converts to quadratic form,
log Note that p= Vr. It is straightforward solve this inequality and we can obtain that the inequality holds when
(10)
> TE
logg.
2 4
A.2 Proof of Lemma 4.2
Proof. First, since each Tik and Zzk independently follow Uniform 1), straightforward to obtain [[ - 2zxx = Let define 7 = ||x; -xx
x,, 112 Then we have B[y] dd.. According to
Hoeffding's 's inequality, for all t > 0,
2t
p(h- σd
22 +) < 2cxp
(11)
Since p we have p(p > = pl7 > 2()) Since d > 668-2, we set tdd GJ2"
0, and
apply to obtain
2 ( -2--2 (5
(d - 2)
p(h a 1- 15-1 s exxp
2expp
6
18d
which is equivalent to
6
1
p(h-
d d 6T- > 1- 1
(12)
6
6
18d
Since p(> > B72)
2 pl8
<9 <
d 22-2 = p(1 - 1/ll < 6
f272)
combining with
12
3
obtain the bound
2exp (_ (d (d- [d-66f)>)
2
p(p> >
pl > 23
18d
A.3 Proof of Propusition 4.3
Proof. To solve p² 2/ < f it is equivalent to solving
1 2
< E, r= V5p,
5 er
which is further equivalent
1
log
+ 2logr r - r < log
5
13

---

Published as a conference paper at ICLR 2025
Fun(d, d')
Structure
Ackley(300, 150)
Partial variable dependency
Partial variable dependency
(0senbrok(300, 100)
Nonoverlap additive decomposition
H4r0mamm6((00,6)
Partial variable dependency
Rosenbrock(101. 100)
Nonoverlap additive decomposition
Stybtang (200, 200)
Acklcy(150, 150)
None
Table 3: Synthetic benchmarks and the structures within the target function: d is the input dimension, and d' is the number of effective variables compute the target function.
leverage logr < V for 1, and solve the upper bound,
log + 2Vr r < log
which gives
1 (1 / + log 1/5 log
P Matérr
(13)
v5
A.4 Proof of Lemma 5.1
Proof. First, similar to the proof in Section A.2 Hoeffding's 's inequality, we have for all > O,
under the uniform distribution on [0, 1J", from
1 d 21
212
p(h
2exp
(14)
When we set lo = c such that > 1/Vōr, we have r3 2 d/6 > 0, and apply the Hoeffding's ' inequality
d/6. Let us choose 17-2 d/6 =
p - d/66 er2 d/6) < 2exp
2(E2r2 d/6)2
exp -2
1 -2(6672--/3d (15)
6
Since p(p 2 p(7 2) S p(l- d/6]
lg d/6), obtain
-22(c272- 1
-
6 (2d)
p(p 2 7) s 2exp
B More Experiment Details
B.1 Definition of Synthetic Functions
Stybtang. We used the following hightly-mdifiiee Stybtang function.
f(x)
Σ (( - -c, 16(x; - c, + 5(x; - -ci)),
2
where x E [-5,5] 5)P, and we set (cl...., cp) as an evenly spaced sequence in [0, 7.5] (c) = 0 and
cD 7.5). The optimum is xt |c 2.903534.....44 2.903534).
Rosenbrock. We used the following Rosenbrock function,
D-1
f(x) = Σ
100 (1x441 - G+1) - (xx - -c,)2)2 + (1 (x; - c;))
14

---

Published as a conference paper at ICLR 2025
SBO-Matern
VBO
SBO-SE (RI)
Median Baseline
SBO-Matern (RI)
SBO-SE
Best Baseline
Worst Baseline
Rosenbrock (100,100)
Ackley (150,150)
StybTang (200,200)
0
0
-7.5
50000
-5000
10.0
-100000
50000
12.5
-150000
15.0
-55000
-200000
-250000
17.5
1000000
-300000
-00.0
125000
200
400
200
400
0
200
400
Rosenbrock (300,100)
Ackley (300,150)
Hartmann6 (300, 6)
-50000
10.0
3-
1000000
12.5
2
150000-
1550
200000
17.5
1
-250000
-00.
-300000
0
200
400
200
200
400
Figure 5: Runtime performance: Maximum Function Value Obtained vs. Number of Steps in all synthetic benchmarks. For cleaner view, we plotted the full results only for SBO and VBO. the remaining methods. we show best, median, and worst result at each step.
where [-2.048, 2. We set [ei.... cp] as an evenly spaced sequence in [-2, 2). where 22 and cp = 2. The optimum is at xt (cr + 1..., cD + 1].
Ackley. We used the Ackley function defined in (Surjanovic Bingham
d cos(2
1
f(x) = 20 exp
0.22
d
x
exp
d
20 exp((),
where each v, 02.32.7788 32.768) and the (global) optimum at x 10.....00 0]. Hartmann6. The function is given by
f(x)
Σ 1=1
a,exp
Σ 1=1
Agla, Pij)2
where each E [0, 1J. A = [Ag]] and P = [Ptj] are two pre-defined matrices in (Surjanovic &
Bingham). The global optimum is at xt = (0.20169. 0.150011, 0.4768747 0.275332, 0.6573).
B.2 Humanoid-Standuy Benchmark
We created novel trajectory optimization benchmark in Humanoid Standup task based on the MuJoCo physics engine (Todorov et al. 2012).. The action for the environment determined by 17 parameters (corresponding to 17 motors in the humanoid). We set the trajectory length to 59, Therefore, the dimension of the problem 59 x 17 = 1, .003. The goal is to find a trajectory (aj..... asg) of motor actions that maximize the reward, for given initial state. This is an instance of optimal control and planning, which is a classical problem in reinforcement learning. noteworthy that the recent work of Hvarfner al. 2024) also created BO benchmark (named
15

---

Published as a conference paper at ICLR 2025
SBO-Matern
VBO
SBO-SE (RI)
Median Baseline
MBO-Mn (RI)
SBO-SE
Best Baseline
Worst Baseline
Mopta08 (124)
Rover (60)
NAS201 (30)
-220
2.5
0.74
240
0.0
-60
0.72
2.5
-280
5.0
0.70
-300
775
-320
10.00
0.68
-340
112.5
100
200
300
0
100
200
300
100
200
300
DNA (180)
SVM (388)
Humanoid-Standup (1003)
-0.05
140
-0.30
120
0.10
100
-0.2
80
00.15
60
-0.20
-0.6
100
200
300
0
200
400
600
800
100
200
300
Figure 6: Runtime performance: Maximum Function Value Obtained vs. Number of Steps all real-world benchmarks. For cleaner view, we plotted the full results only for SBO and VBO. For the remaining methods, we show best, median, and worst result at each step.
"Humanoid") based on the MuJoCo engine, but the problem setting is different. Their benchmark seeks to optimize linear policy, represented by 376 x 17 parameter matrix that linearly maps the humanoid state to an action at each step. In contrast, our benchmark makes no assumptions about an underlying policy and directly optimizes the entire trajectory, making it policy-free.
B.3
Implementation
ALEBO. We used ALEBO implementation shared by the Adaptive Experimentation (AX) Platform (version 0.2.2). The source code at lattps://githu/ htt.atuuuu.cmm//ccebbbooll
Ax/blob/ wwt.loo/maiin//m aataainnaa/iiieeaa http////aaaaaaaa//////e///////// http:/mmoells/tt/m
HESBO. We used HESBO implementation of the original authors (attps://githu)
com/ htt////nnnnybieeeee))
TURBO. We used the original TURBO implementation (attps://github.o)
uber iber-research/TuRB:
SaasBO-NUTS, We used the original implementation of the authors (attps://githun conn/naatiinjjnnoowiakk/ sat.:/oooiiaksaassooooooooooo
SaasBO-MAP. The SaasBO implementation available to the public does not include the version using MAP estimation. We therefore implemented this method based the original paper (Eriksson & Jankowiak] 2021). All the hyperparameter settings follow exactly the same the original paper.
RDUCB and Tree UCB. The implementation of both methods is publicly available at https::/github. www/huaww........ www.uaaweinohhhaaeee httpaaa//...eeeeeeee//eeeeeeee/e
BNN Laplace and BNN HMC. For BNN Laplace, we used the implementation at
bttps::/giitani com/ wi http://oo//pppppeeeeaaeeeeeeeeee
and for BNN HMC,
used the implementation at hitps:://///......... /tttps::/qithtee http:ttwww..cccc.ucaa////maoaaae the architecture of the neural network, perform leave-one-out cross-validation the
To identify
16

---

Published as a conference paper at ICLR 2025
SBO-ARD
BNN Laplace
Best Baseline
Worst Baseline
S80-Single-1s
BNN HMC
Median Baseline
Mopta08 (124)
Rover (60)
NA5201 (30)
DNA (180)
-220
2.5
0.74
240
0.0
-0.30
260
280
-0.32
5.0
0.70
300
-75
-0.34
320
-10.0
0.68
340
100
200
300
-2255
100
200
300
100
200
300
-0.36
100
200
300
SVM (388)
Rosenbrock (100,100)
Ackley (150,150)
StybTang (200,200)
0
-7.5
-0.14
-50000
-10.0
25000
-0.16
1000000
12.5
500000
150000
15.0
-50000
0.18
2000000
-17.5
100000
250000
-300000
-20.0
-125000
100
200
330
0
200
400
0
200
400
0
200
400
Rosenbrock (300,100)
Ackley (300,150)
Hartmann6 (300. 6)
-50000
-10
3-
-12
1100000
-14
2-
1500000
-16
200000
-18
-250000
-20
-300000
0-
200
400
0
200
400
0
200
400
Figure 7: Runtime performance comparison with BNN- BO methods and SBO with single length parameter. Matém kernel was used conduct SBO. For clarity. VBO was excluded comparison. The results were averaged from 10 runs. The shaded region depicts the standard deviation.
initial dataset for each task. The layer width and depth were selected from (1, 2, 3} X (32, 64, 128, 256). The activation function for BNN Laplace is ReLU and for BNN HMC is Tanh, which is the default choice of each method. For BNN Laplace, the training was conducted by maximizing marginal likelihood with ADAM. with learning rate of 1E-01, and weight decay of 1E-03. For HMC, the BNN is pre-trained by maximizing log likelihood with ADAM. Then followed by HMC sampling procedure. We used their adaptive step size schedule for the Leap frog, and the number of steps set 100.
GP Training in Standard BO. For efficiency, we trained GP via point estimation. The optimizer was chosen from L-BFGS, Adam, and RMSProp. For Adam and RMSProp, we set the initial learning rate to 0.1 and 0.01, respectively. The maximum number of epochs was set to 1,500 for both methods. We maximize the log marginal likelihood of the GP regression model. The positiveness of the length-scale and noise variance parameters ensured via SoftPlus transform. For L-BFGS, we found that maximizing the log marginal likelihood often incurs numerical issues. To achieve numerical stability, employed the prior Un form(0.001. 30) over each length-scale, a diffused Gamma prior over the noise variance Gamma 1I.1, 0.05), and the prior Gamma(2.0, 0.15) over the amplitude. The maximum number of iterations was set to 15K and tolerance level 2E-9.
B.4 Computing Environment
We conducted the experimental investigations on a large computer cluster equipped with Intel Cascade Lake Platinum 8268 chips.
17

---

Published as a conference paper at ICLR 2025
Method
Structural Assumption
ALEBO, HESBO
Low-dim embedding
SaasBO
Partial variable dependency
RDUCB, Tree UCB
Additive function decomposition
SBO. TURBO, BNN
None
Table 4: BO methods and their structural assumption about the target function.
Additional Results
C.1 BNN-based BO and Standard BO with a Single Length-Scale
Due to the space limit, we did not report the results of the BNN-based BO in the main paper. Here, we supplement the comparison with the aforementioned two BNN-based BO approaches, BNN Laplace and BNN HMC, as specified in Section B.3 In addition, we also tested SBO with single length- scale. We used the Matérn kernel. The comparison was performed across all the synthetic and real-worlh benchmarks except Humanoid-SSandiin
( (1003). The results are reported in Fig. 7 Across all the benchmarks, the performance of the BNN based methods is worse than the median of GP based methods. In cases, such as SVM, Hartmann6 (300, 6) and Rosenbrock, the performance of the BNN based methods is close to or overlapping with the performance of the worst GP baselines. This might be due to that under low data regimes, it is more challenging for neural network to capture the landscape of the target function and quantify the uncertainty, especially for high-dimensional targets.
Interestingly, even using a single length-scatt scale parameter, the standard BO, namely SBO-Single as long as it does not encounter gradient vanishing - still consistently shows superior performance, which often elosely matches ARD Matérn, except for Ackley(300, 150), SBO-Single struggled in finding right locations to query, leading to poor performance.
C.2 Different Choice of in UCB
We examined how the exploration parameter À in UCB influences the BO performance in high- dimensional optimization. To this end, used the Maérn kernel and varied À from (0.5, 1.0, 1.5, 2.0}. We show the runtime optimization performance in Figure g can be seen that, in most benchmarks (Mopta08, Rover, NAS201, DNA, etc.), the performance with different X's is close. However, on SVM and Ackley, the performance of X = 0.5 is typically worse than the bigger choice of X. This might be due to that those problems more challenging (e.g., Ackley has many local minima), a larger exploration level (A) is needed to identify better optimization outcomes.
C.3 Alternative Acquisition Functions
Finally, we evaluated three other acquisition functions for nigh-dimensional optimization: Expected Improvement (EI), log-EI and Thompson sampling (TS). At each step, EI computes the expected improvement upon the best sample so far,
El(x) = (p(x) - f(x")) [(()= f(x*) Vo(x)
(u(x) - f(x*) Vu(x)
+ u(x
(16)
where (-) and ((-) are the CDF and PDF of the standard Gaussian distribution, x" is the best sample that gives the largest function value among all the queries up to the current step. p(x) and /o(x) are the posterior mean and standard deviation of fxx given current training data. We used BOTorch to maximize EI through L-BFGS. The most recent work (Ament et al. 2024) introduced log-El to overcome the numerical challenges in optimizing the original HI acquisition function. Accordingly, we also tested log-El using the BoToorc implementation. To conduct TS at each step, we used Sobol Sequence (Sobol" 1967) to sample 3K inputs the domain, and jointly sampled the 3K function values from the GP posterior. We then selected the input with the largest function value as the next query. We denote this method by TS-3K. To explore more input candidates, we also employed the Lancos and conjugate gradient descent to approximate the posterior covariance matrix, and sampled 10K input candidates from Sobol Sequence. We denote this method as TS-10K We used BOTorch implementation of both methods. We used the Matérn kernel.
18

---

Published as a conference paper at ICLR 2025
A=1.5
A=1.0
Best Baseline
Worst Baseline
A=
2=2.0
Median Baseline
Mopta08 (124)
Rover (60)
NA5201 (30)
DNA (180)
-220
2.5
0.74
0.0
-0.30
-760
-80
-0.32
0.70
300
-75
-20
100.0
0.68
-0.34
340
100
200
300
12.5
100
200
300
100
200
300
-0.36
100
200
300
SVM (388)
Rosenbrock (100,100)
Ackley (150,150)
StybTang (200,200)
-7.5
-0.14
50000
-10.0
25000
1000000
12255
550000
-0.16
15.0
75000
2000000
-0.18
-17.5
1000000
250000
-0.20
-300000
-20.0
122000
100
200
300
0
200
400
0
200
400
0
200
400
Rosenbrock (300,100)
Ackley (300,150)
Hartmann6 (300. 6)
-50000
10.00
3-
1000000
12.5
150000
200000
17.5
1
250000
20.00
-300000
0-
200
400
0
200
400
0
200
400
Figure 8: Runtime optimization performance with different exploration level in UCB. In the main paper and Figure and SBO was compared with other methods using 1.5.
The ultimate and runtime optimization performance for using EI, TS-3K and TS-10K, are shown in Fig. and IO, respectively. As we can see, in most cases, such as Mopta08 (128), SVM (388), DNA (180) and Rosenbrock (100, 100), both El and TS perform worse than UCB, and also worse than the best baseline. On Ac5ley(150, 150) and Ackley(300, 150), the performance of El and TS even among the worst. The failure of EI might arise from the gradient vanishing during the optimization, as pointed by the recent work (Ament et al.) 2024). The failure of the TS might be due to the high-dimension of the input space. A few thousand candidates might not be sufficient cover the whole space; as result, the selected queries are inferior in indicating the best direction to optimize the target function. The use of log-El performed excellently, achieving results comparable to those obtained with UCB.
D Length-scale Initialization Strategies Across Different Gaussian Process Libraries
In this section, have investigated the initialization strategy of a (non-exhaustive) list of GP training libraries. Most libraries employ a small constant initialization while a few others created the initialization based on the distances between the training data points.
D.1 Constant Initialization
Gpytorct initializes the length-scale as Sof/Plux0000==0.699.
http:://ggithbb..cmm///mmel http:////o..mmlll1su-ooo//tttoo
19

---

Published as a conference paper at ICLR 2025
UCB
LOG-EI TS-3K
TS-10K
Median Baseline
"*/- Best Baseline
Worst Baseline
Mopta08 (124)
Rover (60)
NAS201 (30)
DNA (180)
0.79
-1.0
-0.0
Pin car
0.2
-00
-0.0
SVM (388)
Rosenbrock (100,100)
Ackley (150,150)
StybTang (200,200)
1.13
-11
-0800
00/0
111
410
-83
Rosenbrock 1300,100)
Ackley (300,150)
Hartmann6 (300, 6)
Figure 9: Final performance with different acquisition functions.
GPy initializes the length-scale 1.
Scikit-Learf initializes the length-scall scale as
GPMI initializes the length-scale as exp(0)
GPstuf initializes the length-scale 1.
GPJax initializes the length-scale 1.
Spearmin initializes the length-scale as
GPllow initializes the length-scale 1.
D.2
Data- Date-Dependent Initialization
initializes the length-scale randomly based on the range of the input features. The initial length scale for each feature is sampled from a uniform distribution, with lower bound of 1 x 10-11 10 and an upper bound set to twice the maximum difference of the variable values in the training dataset.
DiceKriging 10
learn -earn
GPJax
20

---

Published as a conference paper at ICLR 2025
UCB
LOG-EI
TS-10K
Median Baseline
TS-3K
Best Baseline
Worst Baseline
Mopta08 (124)
Rover (60)
NA5201 (30)
DNA (180)
-220
2.5
0.74
0.0
-0.30
-760
-200
-0.32
0.70
300
-75
-0.34
120
110.0
0.68
340
100
200
300
12.5
100
200
300
100
200
300
-0.36
100
200
300
SVM (388)
Rosenbrock (100,100)
Ackley (150,150)
StybTang (200,200)
-7.5
-0.14
50000
-10.0
25000
-0.16
1000000
12.5
550000
15.0
75000
2000000
-77.5
2550000
-300000
-20.0
1120000
100
200
300
0
200
400
0
200
400
0
200
400
Rosenbrock (300,100)
Ackley (300,150)
Hartmann6 (300. 6)
-50000
-10
3-
-12
1100000
-14
150000
-16
200000
-18
-250000
-20
300000
0-
200
400
0
200
400
0
200
400
Figure 10: Runtime performance with different acquisition functions.
hetGr initializes the length-scale to the 50th percentile of the pairwise distances between the training data points.
GPVecchia initializes the length-scale to be proportional to the mean of the pairwise distances between the training data points.
E
Limitation
Our current work has identified gradient vanishing as critical failure mode of SBO high- dimensional settings. We have proposed an effective, robust initialization method to mitigate this issue. However, we have not yet analyzed whether there is an upper limit to the input dimension for which SBO can function properly, what that limit might be. To answer these questions, will continue our investigation in future work, both theoretically and empirically. In addition, we empiri- cally found that the choice of optimizer can also impact GP training performance in high-dimensional settings, For instance, L-BFGS tends to be less stable than ADAM, while RMSProp can outperform ADAM in certain cases. In future work, we will continue exploring initialization strategies to enhance the robustness of GP training across different optimizers.
hitrass///ittbbb.cc//rrr///eee
21


Reviews:
[{"review_id": "3RMPu305Qu", "reviewer": "Reviewer_M8nK", "summary": "This paper considers the problem of high-dimensional Bayesian Optimization with Gaussian Processes, in particular the selection of a kernel. It argues that the poor performance of the standard SE kernel in high dimensions stems from the vanishing gradient problem, and provides theoretical and experimental evidence to this effect. It then proposes a method to initialize the parameters of the SE kernel (setting the length-scale parameters to $O(\\sqrt{d})$) to tackle this problem. Experimental results on synthetic and real datasets show that this initialization brings up the performance of the SE kernel to be comparable to the state-of-the-art.", "questions": "Have you done experiments on problems with even higher dimensions, e.g. in the thousands?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- The paper is well-written, clear, and straightforward to understand.\n- Hyperparameter optimization is a problem with wide-ranging applications.\n- As far as I know, the consideration of the gradient vanishing problem for Bayesian Optimization with Gaussian Processes is novel. It is quite interesting, and indicates that it might be a pertinent direction of analysis for other applications due to the popularity of gradient descent.\n- The paper provides both experimental and learning theoretic evidence of the gradient vanishing problem as the problem dimension increases.\n- The proposed fix brings up the performance of standard BO with GPs and classic SE kernel to be comparable to the state-of-the art, including on NAS201.", "weaknesses": "- The dimensionality of the problems in the experiments are fairly high, but some problems in neural architecture search have even higher dimensionality.\n- The proposed method does not outperform the baselines. However, I don't see this as a major weakness since I feel the main message of the paper is not the proposed method, but the importance of considering the gradient vanishing problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "78SeXWXDjh", "reviewer": "Reviewer_pG21", "summary": "The authors investigate vanishing gradient problems associated with the squared exponential and Matern kernel functions, finding them both to exhibit gradient vanishing at common default initializations. They propose a dimension-dependent initialization and advocate for the Matern kernel, and demonstrate on numerical experiments that this can actually be more important than all the bells and whistles we have been putting on GPs for the last decade.", "questions": "Why not use a Matern kernel with the proposed initialization in your numerical experiments?\n\nUnder the assumption of uniform sampling, does any of the procedures of DiceKriging, hetGP or GPVecchia satisfy your robust initialization?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The aim of this work, to challenge a basic, un-sexy component of the BO framework as the lengthscale initialization/covariance kernel and to attempt to pin down the influence on model performance is great. \n\nAt first, I found it difficult to fully take in the implications. But if my quick refreshing/skim of them is to be trusted, the articles proposing TuRBO, Tree UCB, RDUCB, HESBO, ALEBO, SaasBO do not have any \"vanilla BO\" baseline: they compare only to other sophisticated methods. This article may have revealed that the emperor has no clothes: this thread of literature may indeed have been missing the basic Matern kernel as a suitable approach to high dimensional BO. If correct, this finding will have a significant impact on the future of BO research in high dimension. It should probably also prompt some soul-searching among us who have been doing work in this area.", "weaknesses": "---\n\nInsufficient replicates in the numerical experiments. \n\nCan a study with 10 replicates really be said to be a \"comprehensive evaluation\"? Particularly when there is such a large variance in the results: maybe half of the method-problem combinations seem to span the entire plot. Given how influential the conclusions of this article would be if correct, it is really disappointing that they are based on estimates of performance with such high standard error. Running (significantly) more experiments would solidify the important findings of this article. Really, there should probably be closer to 50 replicates or more.\n\n\n---\nCriticism of the squared exponential and advocacy of the Matern kernel is not new. \n\nStein, M. L. in Interpolation of Spatial Data makes this argument in 1999, and that is referenced in the covariance chapter of Rasmussen and Williams. \n\nIt is probably worth citing at least these documents, and trying to locate other discussion of this topic.\n\n---\n\nWhile there are indeed (surprisingly) many software platforms which use a constant lengthscale initialization (including GLflow, GPJax, Spearmint, GPStuff, GPML, scikit-learn and GPy, based on my quick skim of their source code), this is not universally the case. \n\nThere are also packages that do indeed use the idea of involving the expected distances between points when choosing the lengthscales.\nThese packages actually do so empirically instead of assuming a uniform input distribution, which is surely violated in, say, Bayesian Optimization settings. The ones I found based on a quick search were:\ni) DiceKriging uses multiple random initializations between 1e-10 and twice the range of the distances.\nii) hetGP initializes the lengthscale based on the quantiles of the observed distances.\niii) GPVecchia seems to use 1/4 times the mean distance between training points.\n\nIt is probably worth mentioning these and situating them relative to your proposal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mDDBWkPyaA", "reviewer": "Reviewer_RmNw", "summary": "This paper thoroughly investigates the performance of standard GP in high-dimensional BO. Specifically, the paper focuses on the use of SE and Matern kernels in standard GP, along with analysis on their gradient values during GP training. The authors point out the failure of using SE kernel is mainly due to the vanishing gradient problem and propose a new initialization strategy for kernel length-scale that mitigates the problem and improves performance. Comprehensive experiments are carried out in comparing performance of the proposed standard BO (SBO) against other state-of-the-art BO algorithms on various optimization tasks.", "questions": "1.  In equation (1) and (2), both ARD SE kernel and ARD Matern kernel contain the amplitude $a$. Since in line 405 a Gamma prior is placed on this parameter, I assume the authors also include $a$ as the trainable parameter in GP training. How did the authors initialize $a$ in GP training? In equation (6), is the amplitude $a$ considered as a constant during theoretical analysis? In line 209-212, why is $a$ not presented in the gradient expression of the Matern kernel?\n\n2.  In Figure 2, the relative length-scale differences for both GP-Matern($l_0=0.693$) and GP-SE($l_0=\\sqrt{d}$) are flat lines. Since $l_0$ is initialized at the same value independent of the number of iteration, does this mean the value of length-scale is all the same after training even for different GP (since the GP at each iteration is modeled on different data sets)?\n\n3.  Consider cases: Ackley(150, 150), Mopta08(124), Rover(60), NAS201(30), DNA(180) in Figure 3 and 4 (5 out of the 11 benchmarks), it seems adding the robust initialization to SE kernel worsens its performance. How would the author justify this observation supports their statement in line 345-346 (\"the performance across all the benchmarks is dramatically improved\")?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The authors claim the vanishing gradient in the SE kernel training is the main cause of BO failure in high-dimensional cases and provide detailed theoretical analysis on the gradient values of SE kernel. The theoretical works are easy to follow and provide insights on how likely the vanishing gradient problem can occur when initial values of length-scale are fixed in high-dimensional setting. The proposed initialization strategy is also supported by a probabilistic bound.\n\n2. The authors perform adequate experiments to validate their statement. Both high-dimensional synthetic functions and real-world applications are considered, and the chosen BO algorithms are representative and classic in the field of high-dimensional BO. In addition to theoretical analysis, the authors also provide several experiments results to confirm both the existence of vanishing gradient problem and its relationship with initial length-scale and dimensionalities.", "weaknesses": "1. This paper mainly focuses on evaluating existing method. Vanishing gradient problem has been widely discussed in gradient-based training. I think the novelty of this submission only comes from the proposed initialization method, but the authors spend only limited paragraph highlighting the performance of this proposed method. In the real-world problems, the proposed SBO-SE(RI) method even worsens the performance of the original SBO-SE in most cases (namely every case except for the SVM(388) according to Figure 4).\n\n2.  I'm not convinced by the experiments design in Section 4, where the authors try to further validate their theoretical results by varying $l_0$ and $d$. (1) The evaluation metric used to confirm the occurrence of vanishing gradient lacks generality. To my understanding neither the $L_2$ difference nor the average gradient norm at the first training step reveals the presence of vanishing gradient throughout the training process. The number of training steps is not specified by the authors but it is safe to say there are more than 1 step of kernel parameter training. Therefore, the numerical difference in $L_2$ norm of length-scale before and after training only proves $l$ was updated in some steps but vanishing gradient could still have occurred in other steps. The latter gradient norm metric is only considered for the first step of training without a generalization to the whole training process. Vanishing gradient could still occur in the following steps. Step-wise gradient/length-scale evaluation would be more convincing. (2) The values of length-scale initialization ($l_0$) is not robustly tested. In line 199-200, the authors state fixed $l_0$ will lead to vanishing gradient and in the later numerical verification (line 239-240) several fixed values of $l_0$ are tested. However, given the minimal chosen value of $d$ is $50$, $\\sqrt{d} > 7$ which is much larger than other choices always. For better generality, I think some larger fixed value of $l_0$ (e.g. $5, 15, 25$) should also be considered as a stronger support for the adaptive $l_0 = c\\sqrt{d}$ setting. \n\n3. The readability of this submission needs to improve. Starting from line 233, I keep being instructed to jump between sections to get details and results of the experiments. For example, in line 236 I need to go to section 6.1 to get definitions of the benchmark notations. Similarly in line 264 I'm referred to Figure 3 and 4 which are at the very end of the paper. It would be much easier to read if the authors could gather all the numerical experiments into Section 6 instead of alternating between theoretical presentation and experiments evaluation (e.g. move line 233-324 to the start of section 6).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IALtLB2qI7", "reviewer": "Reviewer_VLW7", "summary": "The common assumption is that Gaussian Process-based Bayesian Optimization inherently scales poorly with the dimensionality (D) of the search space, such that the search is virtually uninformed for D surpassing 15-20. This paper investigates the validity of this assumption and the reason for the typically observed poor performance on search spaces of high dimensionality. They show that, using the popular SE-kernel, the gradient of the log likelihood with respect to the lengthscales vanishes (with high probability) in the case of high dimensionality and typical initial lengthscales, preventing effective gradient-based training of the lengthscales. While Matérn kernels also experience vanishing gradients, the problem occurs at much larger D compared to SE, which aligns with their experimental performance results. Furthermore, the authors show that initialising the lengthscales to \\sqrt{D} is sufficient to overcome the vanishing gradient issue and achieving performance that matches or exceeds that of more advanced Bayesian Optimization methods.", "questions": "I do not have any questions. Please see the suggestions on how the work could be improved under \"Weaknesses\".", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The overall story of the paper is easy to follow, the investigated topic is important, and the proposed approach is sound. The theoretical analysis (Section 4) of the vanishing gradients is interesting and, as far as I am aware, presents a novel perspective on the curse of dimensionality in this context.", "weaknesses": "The central thesis of this paper is that a standard \"no bells and whistles\" implementation of Gaussian Process (GP)-based Bayesian Optimization (BO) can perform competitively with state-of-the-art methods for high-dimensional BO, requiring only a straightforward adjustment to common practices. However, this thesis closely mirrors that of the existing work (1) \"Vanilla Bayesian Optimization Performs Great in High Dimensions\" (ICML 2024), which presents a similar argument.\n\nIn the paper under review, the authors propose scaling the lengthscales proportionally to \\sqrt{D}, and they achieve this by initialising the lengthscales with \\sqrt{D}. The existing work (1) similarly advocates for scaling the lengthscales in proportion to \\sqrt{D}, as discussed in Section 5.1 of their paper. They suggest implementing this adjustment by setting the hyperparameters in the prior on the lengthscale appropriately, noting that using a maximum a posteriori (MAP) estimate for a hyperprior is a straightforward and widely-used approach.\n\nWhile the reviewed paper does cite the existing work, its claim that the proposed method is \"even simpler\" (see Section 5 \"Alternative Method\") is not entirely convincing. Both papers recommend favouring lengthscales near \\sqrt{D}; the difference lies in how they implement this: the reviewed paper achieves it solely through initialisation, whereas the earlier work suggests adjusting the prior on the lengthscales. As it stands, the reviewed paper presents the main thesis as though it were novel (at least in the abstract, introduction and conclusion), despite its overlap with the prior work.\n\nNevertheless, there are valuable contributions in the reviewed paper, such as the analysis of vanishing gradients, which offers an insightful perspective on the consequences of using inappropriately small (but typically used) lengthscales. To improve the work, the authors could clearly position it as an expansion upon the thesis established in the earlier paper, emphasising the novel aspects of their analysis and findings.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "q9hc5TGJx4", "reviewer": "Reviewer_kZCF", "summary": "This paper challenges the widely held notion that Bayesian optimization, using standard kernels, does not scale beyond problems with a, say, single-digit number of independent input variables $d$. The central hypothesis is that standard-kernel GP models _can_ indeed model empirically relevant high-dimensional optimization problems, but are hampered by a problem of hyperparameter tuning. Simply put, if the length-scales $\\ell_0$ of isotropic kernels are initialized of order 1, independent of $d$, then the implied covariance between a fixed number of function evaluations drops exponentially with $d$, and so do the gradients of the marginal likelihood, quickly reaching a regime in which hyperparameter optimization becomes numerically impossible. \n\nInstead, the authors propose to initialize kernel length scales as $\\ell_0 \\sim \\sqrt{d}$, and show that doing so effectively counters the exponential decay in correlation between datapoints, and the associated gradients of the marginal likelihood. They argue further that this surprisingly simple change actually suffices to make standard-kernel GP Bayesian Optimization competitive with methods specifically designed to address high-dimensional global optimization problems.", "questions": "For the moment I have marked the paper as marginally above threshold, but I would be quite willing to raise my score if the point I raised under \"weaknesses\" is addressed convincingly in the discussion. To re-state it as a question: Are you claiming that practically relevant optimization problems do not have an exponentially growing (in $d$) number of local optima? Or perhaps less explicitly, that practically relevant problems have some other property that means they can be globally optimized with a sample complexity that grows sub-exponentially with $d$? Which part of your work supports this claim? (For example, would you say that the problems listed in Figure 4 are a good surrogate for \"practically relevant\" problems? How so? \n\nAs a counter-example, have you run `SBO-SE (RI)` on a genuine sample of a Matérn $\\ell_0=1$ kernel? (Such samples _do_ have an exponentially growing number of local extrema). What happens in such a case?  \n\nThanks in advance for your clarifying reply!\n\n\n---\n\nUpdate post discussion:\n\nThe authors have addressed my questions. One could always discuss further, and I do see some open questions left here. However, I am convinced this paper provides an empirically highly useful technique, reasonably motivated and supported by theory.\n\nThe presence of a related paper does not seem problematic to me, since the related work is amply cited, and the present paper provides a further simplification that users will find highly usable.\n\nOverall, the results presented here will find direct, valuable use with a significant subset of the ICLR community. I thus recommend accepting this paper.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This work makes a nice, compact analytical point about about an algorithmic detail, and provides a host of empirical evidence to support the hypothesis. The theoretical analysis appears sound, although the overarching argument can perhaps be questioned (more below). Prior work is extensively reviewed, and the text itself is an enjoyable read.", "weaknesses": "Overall I think this is a fine paper. I am not entirely sure I follow the whole argument, though. Here is why:\n\nFrom my understanding, the argument was never quite that standard kernels are particularly bad as modelling high-dimensional functions, but that general high-dimensional functions are (in the worst case) just exponentially hard to optimize, period. One can \"hide\" a global optimum among an exponentially growing number of local extrema in, e.g., a sample from a $\\ell_0=1$ SE kernel GP (indeed, the analysis provided in this paper indirectly shows that this is true). So, while I completely agree with the authors that the gradients of the marginal likelihood become exponentially suppressed with growing $d$ if the kernel is initialized with $\\ell_0=1$, and that this can be alleviated by initializing $\\ell_0\\sim \\sqrt[d}$, a second point also needs to be true for the paper's argument to work: **Practically relevant** objective functions must have some property that constrains their shape such that global optima remain identifiable in feasible sample complexity, *and* that property must be at least approximately captured by the isotropic kernels when initialized in the way proposed in the paper. I am not entirely convinced that the paper provides sufficient evidence in this regard, but I realise that this is a debatable point. Hence I'm trying to phrase a question below.\n\nOn a minor note, the authors admit quite clearly that the recent work by Hvarfner et al. achieves pretty much the same effect as their method, but just uses an ever so slightly more involved way to initialize $\\ell_0$ by drawing from a LogNormal distribution. Personally, I don't think this is a major issue, since the authors deal with it quite openly. I guess one could argue that this limits the novelty of the work somewhat, though.\n\n\n### mini issues, easy to fix: \n* Section 2, line 091: It's not technically true that BO aims to find new points that are \"ideally closer to the optimum\". They should just provide helpful information about the location and/or value of the optimium. \n* Proposition 4.3 starts with a typo (\"Given $\\forall$\")\n* Eq. 6: the parentheses around the argument of $\\exp$ are too small. Use `\\bigg` or `\\left`", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)