Understanding High-Dimensional Bayesian Optimization
Authors: Leonard Papenmeier
Matthias Poloczek
Luigi Nardi
Venue: ICML 2025 poster (2025)
Keywords: Bayesian optimization
global optimization
Gaussian process
high-dimensional
Primary area: optimization->zeroorder_and_blackbox_optimization

Abstract: Recent work reported that simple Bayesian optimization (BO) methods perform well for high-dimensional real-world tasks, seemingly contradicting prior work and tribal knowledge. This paper investigates why. We identify underlying challenges that arise in high-dimensional BO and explain why recent methods succeed. Our empirical analysis shows that vanishing gradients caused by Gaussian process (GP) initialization schemes play a major role in the failures of high-dimensional Bayesian optimization (HDBO) and that methods that promote local search behaviors are better suited for the task. We find that maximum likelihood estimation (MLE) of GP length scales suffices for state-of-the-art performance. Based on this, we propose a simple variant of MLE called MSR that leverages these findings to achieve state-of-the-art performance on a comprehensive set of real-world applications. We present targeted experiments to illustrate and confirm our findings.

Understanding High-Dimensional Bayesian Optimization
Leonard Papenmeier'
Matthias Poloczek? Luigi Nardi 13
Abstract
assumptions on the objective function that lower the data demand, such as additivity (Duvenaud al., 2011; Kan- dasamy et al., 2015; Hoang et al., 2018; Han et al., 2021; Ziomek& & Ammar, 2023; Bardou et al., 2024) or the exis- tence of a low-dimensional active subspace (Wang et al., 2016; Nayebi et al., 2019; Letham et al., 2020; Papenmeier et al., 2022). Without such assumptions, was widely be- lieved that BO with Gaussian process (GP) surrogate is limited to approximately 20 dimensions for common evalua- tion budgets (Frazier, 2018; Moriconi et al., 2020). Recently, Hvarfner et al. (2024) and Xu & Zhe (2024) reported that simple BO methods perform well on high-dimensional real- world benchmarks, often surpassing the performance of more sophisticated algorithms.
Recent work reported that simple Bayesian op- timization (BO) methods perform well for high- dimensional real-world tasks, seemingly contra- dicting prior work and tribal knowledge. This pa- per investigates why. We identify underlying chal- lenges that arise in high-dimensional BO and ex- plain why recent methods succeed. Our empirical analysis shows that vanishing gradients caused by Gaussian process (GP) initialization schemes play major role in the failures of high-ding -ighndimensional Bayesian optimization (HDBO) and that meth- ods that promote local search behaviors are better suited for the task. We find that maximum likeli- hood estimation (MLE) of GP length scales suf- fices for state-of-the-an art performance. Based on this, we propose a simple variant of MLE called MSR that leverages these findings achieve state- of-the-art performance on comprehensive set of real-world applications. We present targeted ex- periments to illustrate and confirm our findings.
Due its many impactful applications, high-dimensional Bayesian optimization (HDBO) has seen active research in recent years (Binois & Wycoff, 2022; Papenmeier et al., 2023). While the boundaries have been pushed signifi- cantly, the causes of performance gains have not always been thoroughly identified. For example, in the case of the BODi (Deshwal et al., 2023) and COMBO (Oh et al., 2019) algorithms, later work found that the methods benefited from specific benchmark structures prevalent in their eval- uation (Papenmeier et al., 2023). Similarly, an evaluation of Nayebi et al. (2019) showed that the performance of some prior methods is sensitive to the location of the optimum in the search space.
1. Introduction
Bayesian optimization (BO) has found wide-spread adop- tion for optimizing expensive-to-evalue tions that appear in aerospace engineering (Lukaczyk al., 2014; Lam et al., 2018), drug discovery (Negoescu et al., 201 1), robotics (Lizotte et al., 2007; Calandra et al., 2016; Rai et al., 2018; Mayr et al., 2022) or finance (Baudis & Pošík, 2014).
black-box func-
This paper is motivated by the recent call for more scrutiny and exploratory research (Herrmann et al., 2024). By ex- haustive experimentation, we first identify underlying chal- lenges arising in high dimensions and then examine state- offthh-art HDBO methods to understand how they mitigate these obstacles. Equipped with these insights, we propose a simpler approach that uses maximum likelihood estima- tion (MLE) of the GP length scales called MLE Scaled with RAASP (MSR). We demonstrate that MSR is sufficient for state-of-thee-rt HDBO performance without the need for specifying prior belief on length scales as in maximum a- posteriori estimation (MAP). We note that practitioners usually do not possess such priors and instead rely on empir- ical performances on benchmarks. In particular, we change the initialization of length scales to avoid vanishing gradi- ents of the GP likelihood function that easily occur in high- dimensional spaces but, so far, have been overlooked for BO.
While BO has proven reliable in dinedimnsioona settings, high-dimensional spaces are challenging due to the curse of dimensionality (COD) that demands exponentially more data points to maintain the same precision with increas- ing problem dimensionality. Several approaches have ex- tended BO to high-dimensional spaces under additional
'Department of Computer Science, Lund University, Lund, Sweden " AAmazon (This research does not relate to Matthias' work ) DBtune. Correspondence to: Leonard Papenmeier <lconard, pameamaepapenmmiireeee. cs...h..s............
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Understanding High-Dimensional Bayesian Optimization
Furthermore, we provide empirical evidence suggesting that good BO performance on extremely high-dimensional prob- lems (on the order of 1000 dimensions) is due to local search behavior and not to well-fir fit surrogate model. In summary, we make the following contributions.
to the input space.
1. We identify underlying challenges that arise in HDBO and explain why recent methods succeed. We show that vanishing gradients and local search behaviors are important in HDBO.
High-dimensional BO in the input space. In the liter- ature, HDBO operating directly in the high-dimensional search space X' is often considered infeasible due to the COD. Numerous approaches have been proposed, often leveraging assumptions made on the objective function f such as additivity (Kandasamy et al., 2015; Gardner et al., 2017; Wang et al., 2018; Mutny & Krause, 2018; Ziomek & Ammar, 2023) or axis-alignment (Eriksson & Jankowiak, 2021; Hellsten et al., 2023; Song et al., 2022), which sim- plify the problem and improve sample efficiency if they are met. Other methods identify regions in the search space relevant for the optimization, for example, using trust re- gions (TRs) (Regis, 2016; Pedrielli & Ng. 2016; Eriksson et al., 2019) or by partitioning the space (Wang et al., 2020),
2. We find that MLE of GP length scales suffices for state- of-the-art performance. We propose a simple variant of MLE called MSR.
3. We evaluate MSR on a comprehensive set of real-world applications and a series of targeted experiments that illus- trate and confirm our findings.
2. Problem Statement and Related Work
We aim to find xx E arg min [()), where s: X' - IR is an unknown, only point-wise observable, and expensive- to-evaluate black-box function and X' = [0, 1]d the d- dimensional search space, sometimes called "input space". Bayesian optimization (BO) is a popular approach to op- timize problems with the above characteristics. We give a summary of BO and Gaussian processes (GPs) in Ap- pendix A and restrict the discussion to high-dimensional BO.
Recently, multiple works re-evaluated basic BO setups for high-dimensional problems, presenting state-of-the-art per- formance on various high-dimensional benchmarks with only small changes basic BO strategies. Hvarfner et al. (2024) use a dimensionality-scale. log- normal length scale hyperprior that shifts the mode and mean of the log-normal distribution by a factor of vd. designed to counteract the increased distance between randomly sampled points. To optimize the AF, they change BoTorch's (Balandat et al., 2020) default strategy of performing Boltzmann sampling on a of quasi-randomly generated points by sampling over both a set of quasi-randomly generated points and set of points that are generated by perturbing the 5% best- performing points. By perturbing 20 dimensions on average, this strategy creates candidates closer to the incumbent ob- servations and enforces a more exploitative behavior (Regis & Shoemaker, 2013; Regis, 2016; Eriksson et al., 2019). The effect of the sampling strategy was recently revisited by Rashidi et al. (2024). They argue that TuRBO's random axis-aligned subspace perturbations (RAASPs) are crucial to performance on high- -ittemmimniionlll benchmarks and, moti- vated by this observation, derive the cylindrical Thompson sampling (TS) strategy that maintains locality but drops the requirement of axis alignment. Independently of Hvarfner et al. (2024), Xu & Zhe (2024) reported that "standard GPs can be excellent for HDBO" which they show empirically on several high-dimensional benchmarks They use a uniform (4110 33 3,30) length scale hyperprior, which, in their exper- iments, performs superior to BoTorch's Gamma 1'(3, 6) length scale hyperprior.
Extending the scope of BO to nigh-dimensional problems was, for long time, considered "as one of the holy grails" (Wang et al., 2016) or "one of the most important goals" (Nayebi et al., 2019) of the field. Several contribu- tions extended the scope of BO to specific high- -ighndimensional problems. However, for the longest time, no fully seal- able method has been found to extend in arbitrarily high- dimensional spaces without making additional assumptions about the problem structure. The root problem of extending Bayesian optimization (BO) to high dimensions is the curse of dimensionality (COD) (Binois& & Wycoff, 2022) that not only requires exponentially many more data points to model f with the same precision but also complicates the fitting of the GP hyperparameters and the maximization of the acquisition function (AF). The growing demand for train- ing samples stems from increasing point distances in high dimensions, where the average distance in a dddimensional hypereube is vd (Köppen, 2000).
This paper focuses on high- sigh--immnsional Bayesian optimiza- tion (HDBO) operating directly in the search space. Other methods for HDBO include linear (Nayebi et al., 2019; Wang et al., 2016; Letham et al., 2020; Papenmeier et al., 2022: 2023) or non-linear (Tripp et al., 2020; Moriconi et al., 2020; Maus et al., 2022; Bouhlel et al., 2018; Chen al.,
We discuss an earlier preprint (https://arx.l org/ abs /2402. 02746v3). In a later version, presented at ICLR 2025, the authors concurrently to our work developed an ini- tialization strategy similar to the one presented in Section 4.
2020) embeddings to map from low-dimensional subspace
2

---

Understanding Hiih-Dimennione Bayesian Optimization
GP sample
the initial value for the length scale hyperparameter with which the gradient-based optimizer starts when maximiz- ing the MLE (y-axis). The objective function is sampled from a GP with 6/2-Matérn kemel and & = 0.5, i.e., a GP prior sample. We provide additional implementation details in Appendix B. The dashed line shows the default initial length scale of 2 used in GPyTorch. We consider gradi- ents smaller than the machine precision for single floating point numbers 'vanished', The reason is that even after 500 gradient updates, the length scale would change at most by R 6 x 10-5 from the value with which the gradient-base -based optimizer was initialized.
1.0
100
10
length scale
0.8
10
0.6
10-6
0.4
10-8
f default
0.2
threshold
10-10
10-12
102 dimensions
10a
Figure 1. Maximum MLE gradient magnitude for the 50 first gra- dient steps initialized with different initial length scales (y-axis) and problem dimensionalities (zr-axis). With short initial length scales, the gradients vanish even for low dimensions.
Methods to Mitigate Vanishing Gradients. One strategy to counteract the vanishing gradients is to replace MIE with maximum a-posteriori estimation (MAP) by choosing hyperprior on the length scales that prefers long length scales:
3. Facets of the Curse of Dimensionality
This section discusses how the curse of dimensionality im- pacts -igh-dimensional Bayesian optimization (HDBO) and techniques mitigate these challenges.
MAP arg max log fogplyy) HogplylX,00 () + log
(2)
0
evidence
prior
3.1. Vanishing Gradients Model Fitting
MAP maximizes the unnormalized log posterior, which is the sum of the MLL. and the log prior. We sometimes use the terms 'MLL' and 'unnormalized log-posterior' inter- changeably if what meant is clear from the context. The gradient of the recently popularized dimensionality-sehi log- log-normal hyperprior of (Hvarfner et al., 2024) directs the optimizer toward the mode of the hyperprior log p(0) that corresponds to long length scales.
Bayesian optimization uses a probabilistic surrogate model of f' guide the optimization. Gaussian processes (GPs) are the most popular surrogates due to their analytical tractabil- ity. They allow for different likelihoods, mean, and covari- ance functions, each often exposing several hyperparame- ters, including the function variance of. the noise variance oh and the d model length scales e that need to be fitted to the task at hand. In the absence of prior information about the objective function f. maximum likelihood estima- tion (MLE) is commonly used to fit the model hyperparame- ters by maximizing the GP marginal log-likelihood (MLL):
If the gradient-based based optimizer of the MLE reaches suffi- ciently long length scales, the MLL log p(y) X, 0) no longer vanishes. Fig. 2 shows the length scales of GP condi- tioned on random observations of realization drawn from isotropic 0000ddimeniional GP prior with length scale f = 0.5 and 5/2-Matémn kernel when fitting with MLE and MAP. We initialize the length scales with In 2 z 0.69, draw and 10 observations uniformly at random, and maximize the MLL for 500 iterations with MLE and MAP, using a log- normal hyperprior. We average over 10 random restarts. When conditioning on only 1 observation, the MAP opti- mization starts the initialization In 2 and converges the hyperprior mode. For 10 observations, the length scales first move toward the mode but then converge to a different point, which trades off the attraction of the prior mode and the ground truth length scale. With MLE, the length scale does not change because of the vanishing gradients issue.
0% MLE arg max log pl 0),
Here, X are points in the search space X, y are the associ- ated function values, and 0 is the vector of GP hyperparam- eters. See Appendix A for more information.
The MLL is usually maximized using a multi-start gradient descent (GD) approach. A crucial component of fitting GP is choosing starting points for the MLE hyperparame- ters. In this section, we show that an ill-suited length scale initialization scheme can cause the gradient of the MLL function with respect to the GP length scales to vanish for high-dimensional problems. Thus, the length scales remain at the numerical values that they have been initialized and will not be fitted to the objective function.
To "ensure meaningful correlation" in increasingly high- dimensional space, Hvarfner et al. (2024) scale log-normal hyperprior for the GP length scales with the problem's is di- mensionality. Xu & Zhe (2024) pursue a different approach of initializing the gradient- based optimizer of the length scales with large values. Specifically, they posit uniform u110-3, ((10-3,3)) 3,)0) hyperprior to the length scales and initialize
Fig. shows the severity of the vanishing-gradien phe- nomenon. We plot the maximum magnitude (element-wise) of the length scale gradient across 50 gradient updates of an isotropic GP with a 5/2-Matérn kernel as a function of the input dimensionality of the objective function (z-axis) and
3

---

Understanding High-Dimensional Bayesian Optimization
GP sample (d= 1000)
GP sample
0.5
0.40
6
0.35
length scale
length scale
0.4
0.30
0.25
4
0.3
0.20
0.2
0.15
2
0.10
0.05
0.00
0
200
400
25
50
75
100
# gradient updates
dimensions
MAP (1 obs.)
MLE (10 obs.)
MAP (10 obs.)
prior mode
Figure 3. Average distances between the initial and the final didates of logk1 for various model length scales and dimension- alities without RAASP sampling. Values in the gray region are numerically zero. In high dimensions, the gradient of the AF van- ishes, causing no movement the gradient- optimizer.
Figure 2. Gradient-based optimization of a GP length scale with MAP and MLE. When conditioning the GP on only one random observation, MAP converges to the prior mode. We plot + one standard error, which is too small to be visible.
optimized with gradient- based approaches. Thus, these 'flat" areas of the AF also lead to vanishing gradients, but this time of the AF. This affects the selection of the next sample in the BO procedure.
the gradient- based optimization of the automated relevance determination (ARD) kernel with samples from the hy- perprior. Although their method does not scale the hyper- prior mode with the problem's 's dimensionality, the expected length scale of x 15 is sufficiently long to avoid vanishing gradients for the problems studied by the authors.
The LogEI (Ament et al., 2024) AF provides numerically more stable implementation of EI, mitigating the problem of vanishing gradients but not solving it, as we demonstrate next. Fig. 3 shows the average distances the gradient-based optimizer travels for LogEI. The surrogate is a GP with observations drawn uniformly at random. The unknown objective function is realization of the same GP. We run BoTorch's 's multi-start GD acquisition function opti- mizer, initialized with 512 random samples and 5 random restarts, and measure the distances between the starting and end points of the gradient-based optimization of the AF, Average distances increase with vd; thus the plot shows average distances between the starting and endpoints of the gradient- based optimization that are normalized by d 1 Gray regions correspond to numerically zero average dis- tance, indicating vanishing gradients. Vanishing gradients of the AF remain problem, even when using LogEI I and operating in moderate dimensions.
While these methods mitigate the issue of vanishing gradi- ents, increasing length scales by any constant factor does not always solve it, as Fig. indicates. Whether scaling the hyperprior of the length scale with the problem's 's dimen- sionality achieves good fit of the surrogate model depends on the properties of f. The success of the two methods de- scribed above is related to the effect of the increased length scales on mitigating the problem of vanishing gradients. If the underlying function varies quickly, the GP needs to use a short length scale to model f. In such cases, the GP cannot model the function globally, as shown in Appendix C.4.
3.2. Vanishing Gradients of the Acquisition Function Several popular acquisition functions for Bayesian optimiza- tion (BO), such as upper confidence bound (UCB) (Srini- et al., 2010), expected improvement (EI) (Jones al., 1998), and probability of improvement (PI) (Jones, 2001), rely solely on the GP posterior, exhibiting only small variation when the posterior itself changes only moder- ately. In nigh-dimensional spaces, the expected distance between two points sampled uniformly at random increases with vd (Köppen, 2000). Thus, there are typically large "unexplored regions' in the search space where the algorithm has not sampled. Suppose a commonly used GP with con- stant prior mean function and stationary kernel. The GP posterior corresponds to the prior in those regions unless the kernel has sufficiently long length scales. Therefore, the GP posterior and the acquisition surface are flat in those vast unexplored regions. Acquisition functions (AFs) are usually
Methods to Mitigate Vanishing Gradients of the AF. One technique for handling vanishing gradients in the AF optimization is locality. TuRBO (Eriksson et al., 2019), for example, uses trust regions (TRs) to constrain the optimiza- tion to subregion of the search space. Even if the GP surrogate is uninformative, TuRBO performs local search and can optimize high-dimensional problems even if data is scarce. TuRBO also uses random exis-aligned subspace perturbation (R (RAASP) sampling (Rashidi et al., 2024; Regis & Shoemaker, 2013), i.c., with probability of min (1, 20). it replaces the value of each dimension of the incumbent solution with a value drawn uniformly at random within the TR bounds. This process repeated multiple times to
4

---

Understanding High-Dimensinn Bayesian Optimization
GP sample
GP sample
RAASP an
RAASP off
0.32
0.28
100
length scale
0.4
0.24
0.75
75
0.20
0.16
0.50
0.2
0.12
0.08
0.04
0.25
0.00
10
20
30
10
20
30
50
100
50
100
Iterationn
Iteration
dimensions
dimensions
random
1-0.28
f.0.05
/-0.50
Figure Left: Average distances between the initial and final candidates of LogÉI I with RAASP sampling. The vanishing gradi- ent issue decreases, Right: Fraction of multi-start GD candidates originating from the RAASP samples when evaluating LogEl random samples. In high dimensions, RAASP samples are increas- ingly more likely to get picked, even for longer length scales.
Figure S. OTSD for BoTorch's AF maximizer operating on a 100-dimensional space and GPs of various length scales with and without RAASP sampling. The behavior of short-length-scan GPs. reverts local search (f 0.05 in the left panel) with RAASP sampling and to local search without RAASP sampling (e = 0.05 and 0.28 in the right panel). Shaded areas show the standard error of the mean obtained by 10 random repetitions. In the right panel, the blue line masks the black line.
create several candidates evaluated on a realization of the GP posterior, a process known as Thompson sampling. The point of maximum value is then chosen for evaluation in the next iteration of the BO loop. This design choice further enforces locality as new candidates only differ on average from the incumbent in 20 dimensions for d 2 20.
by the lack of gray color. This indicates reduction of the vanishing- -erriigegg issues. We attribute this candidates being created close to the incumbent observations where the AF is less likely to be flat.
Differing slightly from TuRBO's approach, RAASP sam- pling has been implemented in BoTorch's AF maximiza- tion and can optionally be enabled with the parameter sample best. BoTorch augments set of globally sampled candidates with the RAASP samples, re- sulting in twice as many initial candidates. It perturbs the best-performing points by replacing the value of each di- mension with a probability of min (1, 2 by a value drawn from a Gaussian distribution, truncated to stay within the bounds of the search space. BoTorch then chooses the points of maximum initial acquisition value start the GD optimization of the AF.
In general, when the length scales of the GP are short and the dimensionality is high, BO shows a local-search-like be- havior with RAASP sampling and a random search behavior without it. We demonstrate this using the observation travel- ing salesperson distance (OTSD) (Papenmeier et al., 2025), which quantifies algorithm's 's exploration by finding the shortest path connecting all observations made up to cer- tain iteration. The OTSD curve of an algorithm A consis- tently lying above the curve of algorithm B indicates that al- gorithm A is more explorative as its observations are spread more evenly across X. The OTSD is always monotonic in- creasing. Fig. 5 shows the OTSDs for 100-dimensional GPs with different length scales, each initialized with 10 random samples in the design of experiments (DOE) phase and sub- sequently optimized with LogEI and RAASP sampling for 20 iterations. Unless the model length scale is sufficiently long for the AF gradient not to vanish (as for = 0.5 in the right panel of Fig. 5), the AF maximizer picks one of the initial random candidates without further optimizing This is supported by the trajectories for the BO phase (iteration 9) following the trajectory of the DOE phase (iteration 9) for f = 0.05 and & 0.28 in the right panel of Fig. 5. We generally recommend the RAASP sampling method as it improves BO by automatically reverting to local search when encountering flat AFs.
With increasing dimensionality or descending length scale, the starting points for the multi-start GD routine chosen by the AF maximizer are increasingly more likely to originate from the RAASP samples. Fig. 4 (right panel) illustrates this. Here, we draw realizations from GPs, initialized with different dimensionalities (z-axis) and length scales (y-axis). For each realization, pling and plot the percentage of candidates of maximum ac- quisition value originating from the RAASP samples across all candidates. A higher percentage indicates a more 'local' sampling. We further average across five random restarts. The percentage of RAASP candidates with maximum ac- quisition value increases with the input dimensionality and decreases with the length scale. At the same time, those candidates stay close to the initial candidates. This is shown in the bottom right of Fig. 4 (left panel), which shows the average distance traveled by candidate points of the AF maximizer when using RAASP sampling. With RAASP sampling, candidates travel a positive distance, visualized
maximize the AF with RAASP sam-
3.3. Bias- Variance Trade-off for fitting the GP GP models are commonly fitted by maximizing the MLL, either using unbiased MLE estimation or using MAP, which places hyperprior on one or several GP hyperparameters.
5

---

Understanding High-Diiennional Bayesian Optimization
GP sample (d 10)
GP sample (d 50)
GP sample (s= 5)
GP sample (s 50)
0.00
10'
penalty
-0.25
length scale
MAP
0.50
100
0.0
data fit
-1
-2.5
100
200
500
1000
2000
100
200
500
1000
2000
-2
-5.0
# samples
# samples
MLE
ground truth prior mode (MAP)
MLL
MAP
MLE
MAP
MLE
MAP
Figure 6. Average length scales (y-axis) obtained by MLE (blue) and MAP (orange) for different numbers of randomly sampled observations (z-axis) for 10- and for a 50-dimensional GP prior sample. The obtained length scales differ substantially for the higher dimensional function if few points have been observed.
length scale
length scale
d 10
= 100
prior mode
ground truth
Figure 7. The MLL surface (bottom), the penalty (top row), and data fit terms (center) for various length scales and numbers observations. Fewer observations lead to more erratic changes in the data fit term, leading higher variance the length scale estimates unless a prior gives additional shape to the surface.
MLE exhibits a higher variance and sensitivity to noise, particularly when fitting model in high-dimensional spaces with scarce data. MAP, on the other hand, has lower variance in the length-scale estimates but comes at the cost of bias unless accurate prior information is available. The MLL is given by
are shown in the first two rows, respectively; the constant term - log 2n is omitted from the figure. We account for the different number of samples between the left and right figures by scaling the penalty, data fit, MLE, and MAP terms with s We show the surface for 5 (left) and 50 (right) data points. As the length scales increase, the entries of the kernel matrix increase. The determinant KK decays more quickly, and the penalty term decreases, adding a more distinct global trend to the likelihood surface. In the limit, & + oo, the kernel ma- trix becomes matrix of ones, and the determinant becomes zero. In low dimensions, the data fit term decreases for long length scales, but the decreasing penalty compensates for this, resulting in relatively flat MLL surface. The fast de- cay of the data fit term increases the hignallto-noiss ratio, making it easier for the optimizer to converge in 10 than 100 dimensions. This can be seen by comparing the green and brown MLL curves in Fig. for = 50 samples. For more observations, MLL becomes smoother in all dimensions, as indicated by the left vs. right panel in Fig. 7.
-
1 logp(y|X,0) - 2 (K(X, x)+ + 1)
data
(3)
n log 2n
log (og] + of
complexity penalty
The first and second terms are often called data fit and complexity penalty (Williams & Rasmussen, 2006). For more details, see Appendix A. This section explores the bias- variance trade-off between these two popular approaches for GP model fitting.
MLE. Fig. shows the length scale obtained when using (blue) or MAP (orange) to fit a GP surrogate model with an 5/2-ARD-Mare --t-- to a realization drawn from an isotropic GP prior with length scale e and noise term 10-$ We examine 10 and 50-dimensional GP and repeat each ex- periment 50-times. As before, the gradient-based optimizer starts with an initial length scale of 10 We observe that the length scales estimated by MLE vary significantly less for the 10-dimensional function than for the 50-dimensional one. As increase the number of observations which the GP surrogate is conditioned, the variance of the esti- mated length scales decreases.
MAP. MAP allows for incorporating prior beliefs about reasonable values for hyperparameters. However, practi- tioners often do not possess such prior information and hence resort to hyperpriors that reportedly perform well in benchmarks. Karvonen & Oates (2023) criticized this an 'arbitrary" determination of hyperparameters.
The orange distributions in Fig. show the average length scales obtained by MAP with a Gamma(3, 6) prior, which has been the default in BoTorch before version 12. O. We use this prior as it has a substantial mass around its mode 1/3.
The dotted curves in the bottom row of Fig. 7 show the likelihood surface for MLE as specified in Eq. (3). The penalty term 1/ log [K] and the data fit term
6

---

Understanding High- Mighh-Dimensiona Bayesian Optimization
M8 = 124)
Lasso DNA (d 180)
0.40
Table I. Comparison of the simple BO methods used for the em-
best value
300
pirical evaluation.
0.35
Method
length scale scaling
RAASP sampling
250
0.30
MSR
/initiil value)
500
1000
0
500
1000
MLE (scaled) MLE (E = In 2)
/ /initial value)
X
Ant (d 888)
Humanoid (d 6392)
2000
x
X
best value
200
DSP
(prior)
1000
400
600
500
1000
200
400
and we refer to this new BO method as 'MLE scaled'. The second method in the evaluation is the same BO method, but now using the default value & In 2 of GPyTorch as initial length scale in MLE. This value is not scaled with d. Based on analysis of the impact of random axis- aligned subspace perturbation (RAASP) sampling when op- timizing the acquisition function (AF), we combine 'MLE scaled' with RAASP sampling and call the method 'MLE Scaled with RAASP' (MSR). We detail the RAASP sam- pling and the AF maximization in Appendix B. The next method, DSP, is the new default in BoTorch (Balandat et al., 2020), which uses maximum a- a-posteriori estima- tion (MAP) estimate of length scales and initializes the optimization with the mode of the dimensionality-scui length scale prior (DSP) (Hvarfner et al., 2024). Table I summarizes the methods with basic BO setups we use for the empirical evaluation.
SVM (d 388)
Rover (d = 60)
10
best value
0.2
5
0.1
500
1000
500
1000
iteration
iteration
(ours)
SAASBO
(scaled)
TURBO (RAASP)
MLE (t 2n2)
DSP
Bounce
MLE (2- ln2, RAASP)
Figure 8. BO with the 'scaled' initialization of MLE performs com-
parably to the state-of-the-art in HDBO.
simplifying our analysis compared to wider priors, which reduce the difference between MLE and MAP. Compared to MLE. the MAP estimates vary less but exhibit significant bias. This is pronounced for the 50-dimensional GP sample, where the MAP estimates for the length scales revert to the prior mode for 100, 200, 500, and 1000 initial samples. The solid lines in the lower row of Fig. 7 show the surface for the MAP estimation, using the same GP sample as for MLE. The log prior term adds additional curvature, resulting in length scale estimates of lower variance. This is particu- larly noticeable for little data (left column of Fig. 7) and consistent with Fig. 6. With more data, MLE and MAP be- come increasingly more similar, with MAP's log posterior decreasing faster for longer length scales due to the Gamma prior.
We also compare against SAASBO (Eriksson & Jankowiak, 2021), TuRBO (Eriksson et al., 2019), and Bounce (Pa- penmeier et al., 2023). SAASBO has a large computational runtime. Hence, we run it only for 500 iterations and termi- nate runs that exceed 72 hours. That is why SAASBO has fewer evaluations for Ant and Humanoid.
Our benchmarks are the 124-dimensional soft-constrained version of the Mcpta 08 benchmark (Jones, 2008) in- troduced by Eriksson & Jankowiak (2021), the 180- dimensional Lasso-DNA (Šehić et al., 2022), the 388- dimensional SVM (Eriksson & Jankowiak, 2021), the 60- dimensional Rover (Eriksson et al., 2019), and two 888- and 6392-dimensional Mujoco benchmarks used by Hvarfner et al. (2024). The first four benchmarks are noise-free, while the others exhibit observational noise.
4. Discussion
Experimental Sctup and Benchmarks. We propose simple initialization for the gradient- based optimizer used for fitting the length scales of the Gaussian process (GP) surrogate via MLE and evaluate its performance for BO tasks. In what follows, we suppose a BO algorithm with a 5/2-ARD-Matte --ttrn kernel and LogEI (Ament et al., 2024). To address the issue of vanishing gradients at the start of the MLE optimization, we choose the initial length scale as 0. and scale with vd to account for the increasing distances of the randomly sampled design of experiments (DOE) points. vi Thus, the initial length scale used in the optimization
MLE Works Well for HDBO. Fig. 8 shows the perfor- mance of the BO methods on the four real-world applica- tions. Each plot gives the average objective value of the best solution found so far in terms of the number of itera- tions. We show the ranking of the methods according to the final performance in Table 2 in Appendix C. The confidence bands indicate the standard error of the mean. MSR achieves competitive performance across all bench- marks, matching the SOTA DSP. Notably, MSR outper- forms DSP on 124d Mopt. 0 8 and 888d Ant, and performs
T0
7

---

Understanding High-Dimensionn Bayesian Optimization
M8pta08 (d= 124)
Lasso DNA (d= = 180)
Mopta08 (d 124)
Lasso DNA (d = 180)
length scale
4000
10a--
0151
2000
10
500
1000
500
1000
500
1000
500
1000
Ant (d 888)
Humanoid 6392)
Ant (d 888)
Humanoid (d 6392)
length scale
10a
OTSD
10000
10a
500
1000
200
400
500
0
100
200
SVM (d 388)
Rover (d 60)
SVM 388)
Rover (d = 60)
5000
10
015
10a
500 Iteration
1000
500 Iteration
1000
500
1000
0
500 iteration
1000
iteration
MSR (ours)
DSP
MSR (ours)
DSP
MLE (scaled)
random
(scaled)
MLE (1 In2, RAASP)
MLE (2 (n)
MLE (I=
Figure 10. DSP exhibits the least exploration. MLE with fixed initial length scales performs like random search on Ant and Humanoid.
Figure 9. Average length scales of MSR and the other methods. RAASP sampling gives more consistent length scale estimates.
slightly worse than DSP other benchmarks. Bounce (Pa- penmeier et al., 2023) is consistently outperformed by MSR, MLE, and DSP but surpasses SAASBO (Friksson & Jankowiak, 2021), especially on Ant. Although the con- stant length scale initialization (f = In 2) without RAASP achieves satisfactory results on power-dimensionll bench- marks such as 124d 08 and 180d Lasso-DM it fails on higher-dimensional benchmarks like 888d Ant and 6392d Humanoid i d. We attribute this breakdown to vanish- ing gradients as shown in Fig. 13 in Appendix C.2.
the lower-dimensional benchmarks, being in line with our analysis of the bias-variance trade-off in Sec. 3.3. At the beginning of its execution, the BO algorithm that uses MLE with scaled initial length scales ('MLE (scaled)") uses longer length scales than all other methods. The resulting esti- mates vary significantly for the high- nigh-dimenssonal Ant and Humanoid d problems, supporting our analysis in Sec. 3.3 where we study the comparatively high variance of MLE compared MAP.
The length scales obtained by MSR lie between the values of DSP and of 'MLE without RAASP sampling' (green dots, 'MLE (scaled)') for most benchmarks. An exception is Ant where MSR sometimes results in shorter length scales than DSP. Overall, the RAASP sampling, which is the only difference between MSR and 'MLE (scaled)', obtains more consistent length scale estimates.
Fig. compares the mean length scales per BO iteration, averaged over dimensions and 15 repetitions. For the 124- dimensional Mopt aC8 and 180-dimensional Lasso-DNA applications, 'MLE (f = In 2)' learns length scales similar to the other methods. However, this constant initialization strategy fails to make progress from the initial value for the more high-dimensional problems in the bottom row of Fig. 9. We attribute this to vanishing gradients of the marginal log- likelihood (MLL) discussed in Sec. 3 and highlighted in Fig. 13 in Appendix C.2.
RAASP Promotes Locality. Fig. 10 compares the amount of exploration that the algorithms perform through the lens of the observation traveling salesperson distance (OTSD) metrie; see Section 3.2 for details on OTSD. We observe that DSP (blue curves) is the most exploitative method on all benchmarks, being in line with the fact that, after MLE with constant length scale initialization ('MLE (e In 2)'), DSP has the shortest length scales on most benchmarks. The fact that 'MLE (E = In 2)' is the most explorative method, coinciding with the 'random search' line in Fig. 10, de- spite having the shortest length scales explained by our
RAASP Reduces Variance. Fig. 9 shows a surprising behavior for DSP. As one would anticipate, the estimated length scales are typically close to the mode of the hyper- prior at the start of the optimization. However, they then converge to an even higher value on all benchmarks but the dige-dima Humanoid benchmark. Furthermore, the deviation from the prior mode more pronounced for
8

---

Understanding High-Dimensional Bayesian Optimization
analysis in Sec. 3.2: the random initial points for the AF optimization not further optimized if the gradient of the AF vanishes, because the method does not employ RAASP sampling. We observe that 'MLE (é- = 2)' does not learn longer length scales during the optimization for Ant and Sumanoid, as indicated by the horizontal brown line in in Fig. 9. Thus, does not recover later. For Mopt and Lasso-DNA, which have a lower input dimension, the ef- fect is less pronounced because 'MLE (E- = In 2)' sometimes learns longer length scales that avoid vanishing gradients of the AF. MLE with scaled initial length scale values (green curves in Fig. 10) is the second-most explorative method. However, this is not due to vanishing gradients but caused by overly long length scales, shown as green dots in Fig. 9. MSR (red curves) is more explorative than DSP and more ex- ploitative than 'MLE (scaled)', which does not use RAASP sampling. This is consistent with the shorter length scales of MSR (red dots in Fig. 9), confirming that MSR not only yields more consistent length scale estimates but also acts more local than its RAASP- SamsS--nmmiin--iiii free equivalent.
this perspective.
Our work focuses on GP surrogate models but we will ex- plore in how far our findings can be extended to other sur- rogate models, such as random forests or Bayesian neural networks.
Finally, we will explore how our findings help improve the performance of established techniques for HDBO by combining MSR with trust regions (TRs) (Eriksson et al., 2019), adaptive subspace embeddings (Papenmeier et al., 2022), or additive structures (Duvenaud et al., 2011).
Impact Statement
This paper presents work that aims to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
A Akknowledgements
This project was partly supported by the Wallenberg AI, Au- tonomous Systems, and Software program (WASP) funded by the Knut and Alice Wallenberg Foundation. The compu- tations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725
Notes on Popular HDBO Benchmarks. In Fig. 9, all methods converge to similarly long length scales on the M8 and Lasso-DNi benchmarks. This is likely at- tributable a specific property of these popular benchmarks, as we discuss in Appendix D. In short, these benchmarks have simple structure that benefits models with long length scales, posing the risk of algorithms being 'overfitted" to these benchmarks.
References
Ament, S., Daulton, S., Eriksson, D., Balandat, M.. and Bakshy, E. Unexpected improvements to expected im- provement for bayesian optimization. Advances in Neural Information Processing Systems, 36, 2024.
5. Conclusions and Future Work
Our analysis reveals underlying challenges in high- dimensional Bayesian optimization (HDBO) while offer- ing practical insights for improving HDBO methods. We demonstrate that common approaches for fitting Gaussian processes (GPs) cause vanishing gradients in high dimen- sions. We propose an initialization for maximum likelihood estimation (MLE) that achieves state-of-the-ar performance on established HDBO benchmarks without requiring the as- sumptions of maximum a-posteriori estimation. Finally, provide empirical evidence that combining MLE with random axis-aligned subspace perturbation (RAASP) sam- pling reduces the variance of the length scale estimates and yields values closer to the ones learned by DSP, providing a fresh argument for the inclusion of RAASP sampling in HDBO.
Balandat, M., Karrer, B., Jiang, D., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. BoTorch: A frame- work for efficient Monte-Carbh Bayesian optimization. Advances in neural information processing systems, 33: 21524-21538, 2020. URL https: : / /gittutt./...// com/ pytorch/ wwnorrch/botee http///ttttttttttteeeeeeeeeeeeee 12. 0. Last access: Jan 16, 2025.
Bardou, A., Thiran, P., and Begin, T. Relaxing the additivity constraints in decentralized no-regret high-dimensional bayesian optimization. In The Twelfth International Con- ference on Learning Representations, 2024.
In future work, we will continue to carefully vet popular benchmarks and propose novel, challenging benchmarks that preserve the traits of real-world applications. Further- more, this work emphasizes the importance of numerical stability for the performance of Bayesian optimization (BO) in high dimensions. Thus, we propose approaching the de- velopment of models and acquisition functions (AFs) from
Baudis, P. and Pošík, P. Online Black-Box Algorithm Port- folios for Continuous Optimization. In Parallel Problem Solving from Nature - PPSN XIII, pp. 40-49, Cham, 2014. Springer International Publishing.
Binois, M. and Wycoff, N. A Survey on High -ighe-imnnsional Gaussian Process Modeling with Application to Bayesian

---

Understanding High-Dimensional Bayesian Optimization
Optimization. ACM Trans. Evol, Learn. Optim., 2(2), aug 2022. doi: 10.1145/354561.
Machine Learning. In Forty-first International Confer- ence Machine Learning, 2024.
Bouhlel, M. A., Bartoli, N., Regis, R. G., Otsmane, A., and Morlier, J. Efficient global optimization for high- dimensional constrained problems by using the Kriging models combined with the partial least squares method. Engineering Optimization, 50(12):2038-205. 2018.
Hoang, N., Hoang, Q. M., Ouyang, R., and Low, K. H. Decentralized high-dimensional bayesian optimization with factor graphs. In Proceedings of the AAAI Confer- ence Artificial Intelligence, volume 32, 2018.
Hvarfner, C.,, Hellsten, E. O., and Nardi, L. Vanilla Bayesian optimization performs great in high dimensions. In Salakhutdinov, R., Kolter, Z., Heller, K., Weller, A., Oliver, N., Scarlett, J., and Berkenkamp, F. (eds.), Pro- ceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learn- ing Research, pp. 20793-2088 PMLR, 21-27 Jul 2024.
Calandra, R., Seyfarth, A., Peters, J., and Deisenroth, M. Bayesian optimization for learning gaits under uncer- tainty. Annals of Mathematics and Artificial Intelligence, 76(1):5-23, 2016.
Chen, Zhu, G., Yuan, C., and Huang. Y. Semi- Semi-supervised Embedding Learning for High-dimensional Bayesian Op- timization. arXiv preprint arXiv:2005. 14601, 2020.
Jones, D. R. A taxonomy of global optimization methods based on response surfaces. Journal of global optimiza- tion, 21:345-383 -33,, 2001.
Deshwal, A., Ament, S., Balandat, M., Bakshy, E., Doppa, R., and Eriksson, D. Bayesian optimization over high- dimensional combinatorial spaces via dictionary-based embeddings. In International Conference on Artificial Intelligence and Statistics, pp. 7021-7039. PMLR, 2023.
Jones, D. R. Large-Scale Multi- mulii-Disciplinary Mass Optimiza- tion in the Auto Industry. In MOPTA 2008 Conference (20 August 2008), 2008.
Duvenaud, D. K. Nickisch, H., and Rasmussen, C. Addi- tive gaussian processes. Advances in neural information processing systems, 24, 2011.
Jones, D. R., Schonlau, M., and Welch, W. J. Efficient global optimization of expensive black-bor functions. Journal of Global optimization, 13:455-492, 1998.
Eriksson, D. and Jankowiak, M. High- diig--dimensional Bayesian optimization with sparse axis-aligned subspaces. In Un- certainty in Artificial Intelligence, pp. 493-503. PMLR, 2021.
Kandasamy, K., Schneider, J., and Póezos, B. High di- mensional Bayesian optimisation and bandits via additive models. In International conference on machine learning. pp. 195-304. PMLR, 2015.
Eriksson, D., Pearce, M., Gardner, J.. Turner, R. D., and Poloczek, M. Scalable global optimization via local Bayesian optimization. Advances in neural information processing systems, 32, 2019.
Karvonen, T. and Oates, C. Maximum likelihood estima- tion in Gaussian process regression is ill-posed. Journal of Machine Learning Research, 24(120):1-47, 2023.
Köppen, M. The curse of dimensionality. In Sth online world conference on soft computing in industrial applica- tions (WSC5), volume 1, pp. 4-8, 2000.
Frazier, I. A tutorial on Bayesian optimization. arXiv preprint htt.:/807.02//// 2018.
Gardner, J.. Guo, C., Weinberger, K., Garnett, R., and Grosse, R. Diseovering and exploiting additive structure for Bayesian optimization. In International Conference on Artificial Intelligence and Statistics, pp. 1311-1319, 2017.
Lam, R., Poloczek, M., Frazier, P., and Willcox, K. E. Advances in Bayesian optimization with applications in aerospace engineering. In 2018 AIAA Non-Deterministic Approaches Conference, pp. 1656, 2018
Han, E., Arora, I., and Scarlett, J. -igh-ddmensional Bayesian optimization via tree-structured additive mod- els. Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pp. 7630-7638, 2021.
Letham, B., Calandra, R., Rai, A., and Bakshy, E. Re-examining linear embeddings for high-dimensional Bayesian optimization. Advances in neural information processing systems, 33:1 33:1546-1555 1558, 2020.
Hellsten, F. O., Hvarfner, C., Papenmeier, I., and Nardi, L. digh--iimensional Bayesian Optimization with Group Testing. arXiv preprint arXiv:2310.03311. 2023.
Lizotte, D. Wang, T., Bowling, M. H.. Schuurmans, D., et al. Automatic Gait Optimization With Gaussian Pro- cess Regression. In IJCAI, volume 7, pp. 1444949,995 2007.
Herrmann, M., Lange, F. J. D., Eggensperger, K., Casal- icchio, G., Wever, M., Feurer, M., Rügamer, D., Hüllermeier, E., Boulesteix, A.-L., and Bischl, B. Po- sition: Why We Must Rethink Empirical Research in
Lukaczyk, T. W., Constantine, P., Palacios, F., and Alonso, J. J. Active subspaces for shape optimization. In 10th AIAA multidisciplinary design optimization conference, pp. 1171, 2014.
10

---

Understanding High-Dimensional Bayesian Optimization
Maus, N., Jones, H. T., Moore, J., Kusner, M., Bradshaw, J., and Gardner, J. R. Local Latent Space Bayesian Opti- mization over Structured Inputs. In Advances in Neural Information Processing Systems, 2022.
Rai, A., Antonova, R., Song, S., Martin, w., Geyer, H., and Atkeson, C. Bayesian optimization using domain knowledge on the ATRIAS biped. In 2018 IEEE Interna- tional Conference on Robotics and Automation (ICRA), pp. 1771-1778. IEEE, 2018.
Mayr, M., Ahmad, F., Chatzilygerouddss, K. I., Nardi, L., and Krügerr, V. Skill-based Multi-objective Reinforcement Learning of Industrial Robot Tasks with Planning and Knowledge Integration. CoRR, abs/2203. 10033, 2022.
Rashidi, B., Johnstonbaugh, K.. and Gao, C. Cylindri- cal Thompson Sampling for High-Dimensional Bayesian Optimization. In International Conference on Artificial Intelligence and Statistics, pp. 3502-3510. PMLR, 2024.
Mockus, The Bayesian approach to global optimization. In System Modeling and Optimization: Proceedings of the 10th IFIP Conference New York City, USA. August 31 31-September 4, 1981, pp. 473-48) Springer, 2005.
Regis, R. G. Trust regions in Kriging-basee optimization with expected improvement. Engineering optimization, 48(6): 1037-105,, 2016.
Moriconi, R., Deisenroth, M. P., and Sesh Kumar, K. High-dimensional Bayesian optimization using low- dimensional feature spaces. Machine Learning, 109: 1925-1943, 2020.
Regis, R. G. and Shoemaker, C. A. Combining radial ba- sis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization. Engineering Optimization, 45(5):529 555, 2013.
Mutny, M. and Krause, A. Efficient high dimensional Bayesian optimization with additivity and quadrature Fourier features. Advances in Neural Information Pro- cessing Systems, 31, 2018.
Šehic, K., Gramfort, A., Salmon, J., and Nardi, L. Las- sobench: A high- figh-dimensional hyperparameter optimiza- tion benchmark suite for lasso. In International Confer- ence on Automated Machine Learning, pp. 2-1. PMLR, 2022.
Nayebi, A., Munteanu, A., and Poloczek, M. A framework for Bayesian optimization in embedded subspaces. In International Conference Machine Learning, pp. 4752- 4761. PMLR, 2019.
Song, L., Xue, K., Huang, X., and Qian, C. Monte carlo tree search based variable selection for high dimensional bayesian optimization. Advances in Neural Information Processing Systems, 15:28488-28500 2022.
Negoescu, D. M., Frazier, P. I., and Powell, W. B. The Anoiledg-GGradeet
Knowledge- Gradient Algorithm for Sequencing Experi- ments in Drug Discovery. INFORMS Journal on Com- puting, (3(3))346-3633 2011.
Srinivas, N., Krause, A., Kakade, S., and Seeger, M. Gaus- sian process optimization in the bandit setting: no regret and experimental design. In Proceedings of the 27th In- ternational Conference on International Conference on Machine Learning, ICML' 10, pp. 1015-1022, Madison, WI, USA, 2010. Omnipress. ISBN 978 6055899777
Oh, C., Tomozak, J., Gavves, E., and Welling, M. Combi- natorial bayesian optimization using the graph cartesian product. Advances in Neural Information Processing Systems, 32, 2019.
Tripp. A., Daxberger, E., and Hernandez-Lobat J. M. Sample-Efficient Optimization in the Latent Space of Deep Generative Models via Weighted Retraining. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), Advances in Neural Information Pro- cessing Systems (NeurIPS), volume 33, pp. 11259-11272. Curran Associates, Inc., 2020.
Papenmeier, L., Nardi, L., and Poloczek, M. Increasing the scope as you learn: Adaptive bayesian optimization in nested subspaces. Advances in Neural Information Processing Systems, 35:1 1586-11601, 2022.
Papenmeier, L., Nardi, L., and Poloczek, M. Bounce: Reli- able high-dimensional bayesian optimization for combi- natorial and mixed spaces. In Thirty-seventl Conference on Neural Information Processing Systems, 2023.
Wang, L., Fonseca, R., and Tian, Y. Learning Search Space Partition for Black-box Optimization using Monte Carlo Tree Search. Advances in Neural Information Processing Systems, 33:1951 1-95522, 2020.
Papenmeier, L., Cheng, N., Becker, S., and Nardi, L. Explor- ing exploration in bayesian optimization, arXiv preprint DrXiv:2502.08208, 2025.
Wang. Z., Hutter, F., Zoghi, M., Matheson, D., and De Feitas, N. Bayesian optimization in a billion dimen- sions via random embeddings. Journal of Artificial Intel- ligence Research, 55:361-387. 2016.
Pedrielli, G. and Ng. S. H. G-STAR: A new kriging-based trust region method for global optimization. In 2016 Winter Simulation Conference (WSC), pp. 803-8144 IEEE, 2016.
Wang. Z., Gehring, C., Kohli, P., and Jegelka, S. Batched large-scale Bayesian optimization in high-dimensional

---

Understanding High-Dimensional Bayesian Optimization
spaces. In International Conference on Artificial Intelli- gence and Statistics, pp- 745-754. PMLR, 2018.
Williams, C. and Rasmussen, C. E. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Wolpert, D. and Macready, W. No free lunch theo- rems for optimization. IEEE transactions on evolutionary computation, 1(1):67-82, 1997.
Xu, Z. and Zhe, S. Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization. arXiv preprint 47X11:2400.002466:3, 2024.
Ziomek, J. K. and Ammar, H. B. Are random decomposi- tions all we need in high dimensional Bayesian optimisa- tion? In International Conference Machine Learning, pp. 43347-43368. PMLR, 2023.
12

---

Understanding High-Dimensional Bayesian Optimization
A. A Review of Gaussian Processes and Bayesian Optimization A.1. Gaussian Processes
Gaussian processes (GPs) model a distribution over functions, i.e., assume that f is drawn from GP: / (P(m())) k(x, x') + o (===1)) where m and & are the mean and covariance function of the GP, respectively (Williams Rasmussen, 2006). Common kernel functions include the radial basis function (RBF) kernel:
kRBF( x') = of exp (-²) )
(4)
or the 9/zautomated relevance determination (ARD)-Matém) kernel:
Mau (z, x') o (1+1++/
+ 5r ( V5r
[1+
exp (-v5r)
(5)
3
with Ei= (x)=)()) Here, e is d-dimensional vector of component-wise length scales. Thus, the kernel's number of e
hyperparameters (HPs) in Eq. (5) is d+L 1.
Given some training data = {(x11 v1)...., (ÆN, uN)]: X = [xI y= = (81.... UN) )T, the function values y. of a set of query points X. is normally distributed as
X, ~N(K(XX, + , K(X., X)KXXX (X)
(6)
(7)
Let 0 - (0 o7. 07 2} be the set of GP hyperparameters. The GP surrogate is then typically fitted by maximizing the marginal log-like log-likelhhood w.r.t. 0, also known as maximum likelihood estimation (MLE), i.c.
0" E arg max log p(y) X,
(8)
1 - (X 0 y (X, + -, " log 2n
log ply|X,00) =
With a gradient-based approach, this is done by maximizing Eq. (9), which is usually multi-modal and difficult to optimize. The gradient of the marginal log-likelihood w.r.t. 0g is given by
a
1 3 (K o logply X, 0)
0K 00,
1 tr (KK o y (K + +Ko 1
0K
(10)
00,
2
00,
where 28,
0K is the symmetric Jacobian matrix of partial derivatives w.r.t. 0:.
One often endows the GP hyperparameters with hyperpriors and seeks the mode of the posterior distribution, known as maximum a-posteriori estimation (MAP):
0" E arg log ((y) log p(0).
(11)
A.2. Bayesian Optimization
Bayesian optimization (BO) is an iterative approach, alternating between fitting the model and choosing query points. Query points are found by maximizing an acquisition function (AF), e.g., expected improvement (EI) (Mockus, 2005). The BI AF measures how much observing point is expected improve over the best function value observed thus far. It is defined as
EI(x) R(())) [[((x)
v*1+] = (a)( --) -x) (N))))())))))
(12)
with Z = [=x)=v" T and o being the standard normal cumulative distribution function (CDF) and probability density function (PDF), and AN and o being the posterior mean and posterior variance at æ, i.c.,
oxx)
(x(x)= K(x, (I, + +0,1)
(13)
(x) (k(<x, + + -
(14)
13

---

Understanding High--Dimensioal Bayesian Optimization
2000
100
1000
1500
d 10000
Count
1000
500
10 10
10"
10"
EI
Figure 11. Distribution of El values for GPs various dimensionalities. When conditioning the amount of data points and maintaining the length scale as the dimensionality grows, the distribution of El values becomes more peaked.
As discussed by (Ament et al., 2024), El often suffers from vanishing gradients, which only worsens in high-dimensionall spaces due to the plethora of flat regions. This is shown in Fig. 11. Here, we condition GP on 100 random points in [0, 1] for which we obtain function values by drawing from GP prior with an isotropic RBF kernel with E = 10. We evaluate the AF on 2000 points drawn from a scrambled Sobol sequence and plot the histograms for various dimensionalities. As the dimensionality grows, there are more equal El values, indicating flat regions the AF and possible problems with vanishing gradients of the EI function. (Ament al., 2024) propose LogEI, numerically more stable version EI that solves many of the numerical issues with EI.
B. A Additional Implementation Details B.1. Implementation of RAASP
We perturb the top 5% observations using normal distribution with o = 10-3, truncated within X. For d 2 20, we also generate samples with only a subset of dimensions perturbed, cach with 20
d probability. Of 4m random samples, 2m are global samples from a scrambled Sobol sequence, m are local samples around the top 5%, perturbing all dimensions, and m are local samples around the top 5%, perturbing 20 dimensions on average d 2 20, or all dimensions if d < 20. We call the 2m local samples the RAASP samples. The starting points for the gradient-basee optimization of the AF are drawn from the 4m overall samples using Boltzmann sampling (Ament al., 2024).
B.2. Optimization of the Acquisition Function
We use the LogEI AF (Ament al., 2024) for its numerical stability. is maximized by evaluating on 512 scrambled Sobol points, then selecting five starting points via Boltzmann sampling for gradient-based optimization using L-BFGS- B, with up 2000 iterations. The budget of 2000 iterations is rarely exhausted, the optimizer typically converges much earlier (see Fig. 12). This aligns with previous work claiming that El attains maxima close to good observations (Ament et al., 2024; Hvarfner et al., 2024) and our observation that the starting points for the gradient- based AF maximizer predominantly originate from the RAASP samples.
14

---

Understanding Highh-Dmennional Bayesian Optimization
SVM
Humanoid
Mopta08
Lasso DNA
0.008
0.002
0.000
1000
2000
1000
2000
1000
2000
1000
2000
1000
2000
AF gradient
# AF gradient
#F gradient
#AF gradient
AF gradient
w/ RAASP
w/o RAASP
Figure Number of gradient updates for the AF optimization for MSR, and with and without RAASP sampling. RAASP sampling reduces the number of gradient updates.
C. Additional Experiments
C.1. Ranking of Optimization Algorithms
MSR
DSP
Bounce
MLE (scaled)
MLE (& In 2)
Mopta08 (d = 124)
1
2
5
Lasso-DNA (d = 180)
4
1
5
Ant (d = 888) Humanoid (d 6392)
2
4
3
2
1
Table 2. Ranking for the different optimizers on the benchmark problems according to their final performance. SAASBO is excluded from the comparison as was not run entirety optimization; Bounce ran into memory issues on Humanoic and, therefore, does not have rank on this benchmark.
C.2. MLE Gradients for Real- World Experiments
Complementing our analysis Sec. 4, Fig. 13 shows the average absolute gradients of the different MLE methods, including our proposed MSR method. The constant length scale initialization ('MLE (L = In 2))) is the only method consistently exhibiting vanishing gradients on the -i8-dmmensionll Ant and the 6392-dimensional Humanoid problems as depicted by the solid orange lines for those benchmarks. On Mopta 08 and Lasso DNA, all methods have non-vanishing gradients. Furthermore, both MLE methods scaling the initial length scale do not suffer from vanishing gradients Ant and Humanoid.
15

---

Understanding High- Bigeh-Dimnsiooal Bayesian Optimization
Mopta08 (d 124)
Lasso DNA (d 180)
abs avg grad
10
10a
10-2
10-1
500
1000
500
1000
Ant (d 888)
Humanoid (d 6392)
102
abs avg grad
10
10-15
10
10 10-32
0
500
1000
100
200
SVM (d = 388)
Rover (d 60)
abs avg grad
1035
101
1029
10
103
10-3
0
500
1000
0
500
1000
iteration
iteration
MSR (ours)
MLE (/ In 2)
MLE (scaled)
Figure 13.. Mean absolute value of the gradients for the different MLE methods, including the proposed MSR. The constant length scale initialization exhibits vanishing gradients for the high-dimensional Ant and Humanoie d problems.
16

---

Understanding High-Dimensional Bayesian Optimization
C.3. High-dimeensionll Benchmark Functions
By the no-free-lunch theorem (Wolpert Macready, 1997), the relative performance of an optimization algorithm depends on the properties of the problems operates on. Here, we show that for several benchmark problems, state-of-the-art algorithm strictly dominates the other methods.
Table 3 shows the relative performances of CMA-ES, DSP, and TuRBO after 1000 optimization steps on the 100-dimensional versions of the Levy, Schwe fel, 1, and Griewank benchmarks. We evaluate Levy in the bounds [-10, 10), Schwefel in the bounds [-500, 500], and Griewank the bounds [-600, 600] by scaling from the unit hypercube in which the GP operates to the respective bounds before evaluating the function.
To better understand the reason for the performance differences, we study the observation traveling salesperson distance (OTSD) for the different functions, shown in Fig. 14. Plots of the 2-dimensional versions of all three benchmarks shown in Fig. 17 and the performance and OTSD plot of 100-dimensional Levy figure can be found in Fig. 15, CMA ES shows the lowest level of exploration and has the lowest OTSD on the Schwefel function, where it outperforms the two other algorithms. On Gri ewank (see Fig. 16), DSP has the highest average OTSD performs best while the less explorative CMA-ES shows the worst performance. For Levy, both TuRBO and DSP are relatively explorative and outperform CMA ES by a considerable margin. We conclude that more explorative algorithms are advantageous on the benchmarks with a clear global trend like Griewank, which resembles a paraboloid, and Levy, which has a parabolic shape along the æ1 dimension (see Fig. 17). In contrast, the Schwe benchmark is more "stationary" in that point's function value depends less that point's s absolute position in the space. Noteworthy, stationarity assumed by GP models with stationary covariance function, which, by far, are the most common covariance functions for high- -igh-dimensional Bayesian optimization (HDBO). A more local approach such as CMA-Bi ES is beneficial on this highly multi-modal benchmark.
Levy100
10a
1000
best value
OTSD
10a
500
500
1000
iteration
CMA-ES
DSP
TuRBO
best value
OTSD
Figure 15. OTSD (solid lines) and performance curves (dashed lines) of the 100-dimensional Levy function
Figs. 15 and 16 show the OTSD and performance plots for the 100-dimensional Levy and Gr Griewank functions. Fig. 17 shows the two-dimensional versions of the Levy, Gri trink, and Schwefel benchmark functions.
Benchmark
Rank
Rank 2
Rank
Levy100
TuRBO
DSP
CMA-ES
Schwefel100
CMA-ES
DSP
TuRBO
Griewank100
DSP
TuRBO
CMA-ES
Table Relative performances of CMA-ES, DSP, and TuRBO after optimizing for 1000 itcrations averaged over 10 repetitions. No algorithm performs best for all benchmarks.
17

---

Understanding High-Dimennional Bayesian Optimization
Schwefel100
1500
x 10
best value
1000
3 x 10
500
2 x 10
0
500
1000
iteration
*
CMA-ES
.....
DSP
TuRBO
best value
OTSD
Figure 14, OTSD (solid lines) and performance curves (dashed lines) of the 100-dimensional Schwe fel function
Griewank100
best value
10a
400
10a
015D
200
10a
500
1000
iteration
CMA-ES
DSP
TuRBO
best value
OTSD
Figure 16. OTSD (solid lines) and performance curves (dashed lines) of the 100-dimensional Gri ewank function
Levy
Griewank
Schwefel
150
1500
60
100
1000
40
20
50
500
10
500
500
250
-10
-500
-500
250
0
5
0
0
5
500
-500
250
-500
Figure 17. The two-dimensional versions of the Levy, Griewank, and Schwefel benchmark functions used above.
18

---

Understanding High-Dimennional Bayesian Optimization
C.4. Hard Optimization Problems We reiterate that the curse of dimensionality (COD) remains reality and exists even for low-dimensional problems. Fig. 18 shows 100 evaluations made by El on 2-dimensional GP prior sample the benchmark function. There is model mismatch; the length scales of the surrogate model are set to the correct value (L 0.025). However, El operates locally and fails to find the global optimum (marked by red cross).
GP sample, 0=0.025
XXX
x
f(x)
GP mean
Figure 18. LogEI run on a two-dimensional GP prior sample for 100 evaluations. The right panel shows the posterior mean at the end of the optimization, highly multimodal benchmarks. El reverts to local search behavior and does not obtain global optimum (red cross).
19

---

Understanding High-Dimensional Bayesian Optimization
D. Popular Benchmarks Seem Simpler Than Expected
This section examines two popular high-dimensiona benchmarks, the 180d Lasso-DNA and 124d Mopta08. In what follows, we will demonstrate that many input variables seem to have little influence on the objective value. These empirical findings suggest that these benchmarks are not truly as high-dimensional as their nominal number of input variables might suggest, potentially misleading the perceived difficulty of these benchmarks and confounding the assessment of what state-of-the-art algorithms will be able to accomplish in practice.
We begin by collecting the top 10% best points identified by the SOTA algorithm using dimensionaliyy-salee scalee length scale prior (DSP) (Hvarfner et al., 2024) for each benchmark. For each dimension z; separately, we then count how often these points lie on the boundary of the search space. If over half of the points place 2; on the boundary, we label 2; as secondary, and as dominant otherwise. To ensure consistency across runs, we perform 15 repetitions and consider those dimensions that have been identified as dominant in eight or more of the repetitions.
Table 4 reports the number of dominant and secondary dimensions. In both benchmarks, large fraction of the dimensions are classified as secondary, meaning that the best solutions obtained by DSP method set the corresponding input variables to value near the boundary of the search space. Note that our finding is directionally aligned with the 43 active dimensions that Šehić et al. (2022) reported for Lasso-DNA using different method to estimate the number of relevant dimensions.
We further investigate the impact of each set of dimensions by replacing, at each iteration, either the dominant dimensions (faominant or the secondary dimensions (fecomiary) with randomly chosen binary values. The rows faominan and fsecondury in Table 4 indicate the corresponding average objective values (with standard errors in parentheses). Replacing secondary dimensions impairs the result only slightly, whereas randomizing the dominant dimensions yields markedly larger performance drop. A two-sided Wilcoxon test shows that all differences are statistically significant at < 0.001.
Next, we examine how Gaussian process (GP) models account for this structure during HDBO. As illustrated in Fig. 19, the estimated length scales of secondary dimensions are typically large, which is consistent with the observation that they influence the (predicted) objective value less. Intuitively, if a dimension's 's optimal value often lies at the boundary, the GP can assign a very large length scale, learning the trend toward the boundary with few evaluations and leaving more of the budget to explore the truly dominant dimensions. In effect, the BO loop focuses on those dimensions that genuinely drive performance.
Below, we demonstrate that for these problems, BO will set most of the input variables to exactly zero or one -that is, to the boundary of the search space- regardless of whether the GP is fit by MLE or MAP. Thus, we conclude that the task effectively has a lower dimensionality than its nominal number of input variables. While this property was already recognized for LasSo-DNA, our results demonstrate for the first time that also applies to Mopta 88. Furthermore, it is interesting to note that simple GP surrogate fitted via MLE or MAP can leverage this structure, a behavior that had not been previously observed.
Above, we showed that 1) dimensions of the best points observed by MLE and MAP that predominantly lie on the border have significantly longer length scales than the "dominant" dimensions, and 2) most of the dimensions have marginal impact. MLE shows a similar behavior and is omitted for brevity. We further omit the 3s8-dimensional SVM benchmark as is known to have a dimeeineniional effective subspace (Eriksson & Jankowiak, 2021).
To complement our analysis, we show that MLE and MAP assign dominant dimensions mainly values at the boundary. This
Lasso -NA (d=180)
Mopta08 8 (d=124)
dominant
69.33
30.80
# secondary
110.67
93.20
dominant
0.315 (+:4 x 10-4)
328.824 (+ (=:1))
secondary
0.445 (#:4 x 10-3)
430.474 164)
Frand
0.410 (++3.3 x 10-2)
403.428 (+48.976)
Table The number of dimensions identified dominant and secondary and the average function values obtained when replacing the dominant (fi) or secondary (fz) parameters with uniformly random values. The average values and standard deviations for uniformly random points are shown as fane- Replacing secondary parameters harms performance considcrably less than replacing dominant parameters.
20

---

Mopta08 Understanding Hiihh-Dimeniionll Bayesian Optimization
Mopta08 (d= 124)
Lasso-DNA (d = 180)
length scale
20
93 20 second. 30.80 dom.
110.67 second 69.33 dom.
10
500
1000
500 iteration
1000
iteration
dominant dims
secondary dims
Figure 19. The mean average length scales of "dominant" and "secondary" dimensions for the Mopta08 (left) and Lasso-DNA (right) benchmarks for DSP.
Mopta08
Lasso DNA
frac. params. border
0.8
frac. paramss. at border
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
0
250
500
750
1000
250
500
750
1000
Iteration
Iteration
Figure 20. Fraction of dimensions set to value the border (0 1) by DSP. The shaded area shows the standard error of the mean across 15 repetitions.
Lasso-DNA
0.8
0.6
frac. params. border
frac. params. at border
0.6
0.5
0.4
0.4
0.3
0.2
0.2
0.1
0.0
0.0
200
400
600
800
1000
0
200
400
600
800
1000
Iteration
Iteration
Figure 21. Fraction of dimensions set value the border (0 or 1) MLE method. The shaded area shows the standard of the mean across 15 repetitions.
indicates that the GP model actually makes use of the specific characteristics of these benchmarks.
Figs. 20 and 21 further show that BO consistently evaluates a large share of the parameters at the border. Fig. 20 shows this for DSP by (Hvarfner et al., 2024) whereas Fig. 21 uses MLE as proposed Section 4. The general trend is that, during the course of the optimization, increasingly many parameters are evaluated at the border, which consistent with Fig. 19.
21

---

Understanding High-Dimensional Bayesian Optimization
We thus argue that two HDBO benchmarks do not fully capture the complexity of HDBO because of this simple structure. While this is to be expected for Lasso-DNA and SVM, this property has not been discussed for the soft-constrained Mopta08 benchmark introduced in (Eriksson et al., 2019) to the best of our knowledge.
22


Reviews:
[{"review_id": "SMGRFpEtUc", "reviewer": "Reviewer_xFof", "summary": "The paper provides an investigation on the problems of BO in high dimensions (the curse of dimensionality). First is the problem of vanishing gradients when modelling GP, which makes some GP hyperparameters, such as lengthscales cannot be updated. Based on recent work of Hvarfner et al., 2024, the authors find that increasing the lengthscale prior such that it scales with problem’s dimensionality can help mitigate the problem. The second mentioned problem is the vanishing gradients when optimizing AF. The author mentions 2 solutions: numerically stable LogEI acquisition function (Ament et al., 2024) and RAASP sampling strategy from TuRBO (Eriksson et al., 2019). Third, the authors provide some insights on the trade-off between bias and variance when fitting GP, particularly a comparison between MLE and MAP. While MLE has higher variance, MAP exhibits a bias towards the prior, which can be devastating in case of inaccurate prior assumption. Finally, the authors combine these findings into a proposed method MSR (MLE Scaled with RAASP) and compare them against several variants and other recent SOTAs.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful for the reviewer's thoughtful feedback and recognition of our work's strengths, particularly the systematic summary of high-dimensional Bayesian optimization (HDBO) issues and the effectiveness of our proposed methodology. In response to the constructive comments, we have addressed the concerns and further strengthened our contributions.\n\n**Structure of the Lasso-DNA and Mopta08 benchmarks and inclusion of additional benchmarks**\n\nWe recognize the reviewer's concern about the structural risks associated with Lasso-DNA and Mopta08. Since these benchmarks are nearly always included in related work, we included them to demonstrate the effectiveness of our method in scenarios with specific structural properties against other methods. In the spirit of ensuring a broader evaluation, we added the SVM benchmark [1] (https://ibb.co/jkx3MXN1) and will add another benchmark for the camera-ready version at the latest.\n\n**Add baselines that use the RAASP sampling strategy, such as TuRBO, and compare against recent SOTA methods like RDUCB.**\n\nAs demonstrated in [2], the DSP method we compare to already significantly outperforms RDUCB. Given this established performance gap, we have focused our comparisons on DSP. However, we acknowledge RDUCB and other milestones in HDBO in the related work section. To provide additional data points, as suggested by the reviewer, we have now added TuRBO (https://ibb.co/jkx3MXN1). TuRBO is competitive on the Mopta08 and Lasso-DNA benchmarks and has good anytime performance on Humanoid for the first half of the optimization but stagnates in the second half. \n\n**Include a baseline that uses RAASP without scaled length scales to isolate the effect of RAASP.**\n\nWe have included TuRBO as a baseline that uses RAASP without scaled length scales (https://ibb.co/jkx3MXN1). We further compare MLE without scaled length scales to RAASP (https://ibb.co/6J7g8XRn) and will add this result to the appendix.\n\n**MSR is simply a combination of previous works.**\n\nWhile MSR builds on existing techniques, its practical significance lies in effectively combining them to address the challenges of HDBO and provide a straightforward yet powerful approach that practitioners can readily adopt.\n\n**Revise Section 4 to identify and describe each method mentioned clearly.**\n\nWe will revise Section 4 to improve clarity, ensuring each method is identified and described, making the experimental setup and benchmarks easier to follow.\n\n**Potential for incorporating MSR with other HDBO techniques.**\n\nIntegrating MSR with methods such as TRs, adaptive subspace embeddings, and additive structures is a relevant future work direction, which could further enhance the performance and applicability of MSR in various HDBO scenarios. We will add this idea in the future work section.\n\n[1] Eriksson, David, and Martin Jankowiak. \"High-dimensional Bayesian optimization with sparse axis-aligned subspaces.\" Uncertainty in Artificial Intelligence. PMLR, 2021.\n\n[2] Hvarfner, Carl, Erik Orm Hellsten, and Luigi Nardi. \"Vanilla Bayesian Optimization Performs Great in High Dimensions.\" International Conference on Machine Learning.", "claims_and_evidence": "All of the three claims are supported by clear experimental results and also in line with previous paper findings.", "methods_and_evaluation_criteria": "•\tThere are only four benchmark problems used in the experimental results, while common BO papers have at least 6.\n\n•\tTo make it worse, the authors criticize (in Section 4 – Notes on Common HDBO Benchmarks) that LassoDNA and Mopta08 problems have structures that pose the risk for HDBO algorithms but then fail to use other benchmarks to support the proposed methods.", "theoretical_claims": "There is no theoretical claim.", "experimental_designs_or_analyses": "The baseline SOTAs, including SAASBO, Bounce and DSP, are reasonable comparisons. However, I think it would be better if the author could include:\n\n•\tA baseline that also use RAASP sampling strategy. For example, TuRBO (Eriksson et al., 2019) or MSR variant that does not use scaled lengthscale. Note that Bounce (different from its previous version BAxUS, does not use RAASP anymore).\n\n•\tRDUCB (Ziomek & Ammar, 2023) are included, a recent SOTA leveraging additivity assumption of f. This is to confirm the performance against other HDBO approaches.", "relation_to_broader_scientific_literature": "The three key findings are in line with previous results. The authors did quite well to summarize them and illustrate in a comprehensive way.", "essential_references_not_discussed": "To my knowledge, all of the essential related works are discussed.", "other_strengths_and_weaknesses": "Strengths:\n\n-\tThe discussed issues of HDBO are worth investigating and are systematically summarized.\n-\tThe proposed methodology is simple yet effective, even when compared to recent more sophisticated methods.\n\nWeaknesses:\n\n-\tThe authors did not provide any new significant suggestions for any of the three problems. The proposed method is simply a combination of previous works.\n-\tThe experimental results are not significantly strong due to the lack of some baselines and inadequate benchmark problems.", "other_comments_or_suggestions": "Minor:\n\nThe writing in Section 4 – Experimental Setup and Benchmarks are hard to follow. The authors mention “the second method” and “the fifth method”, but do not mention the other methods.", "questions_for_authors": "1.\tWhat do the authors think about incorporating MSR with other techniques from HDBO baselines? As MSR improve the GP modelling and AF optimization, it can be incorporated with several heuristics such as trust region (TuRBO), adaptive subspace embedding (BAxUS/Bounce), additivity structure (RDUCB). These heuristics for HDBO eventually must deal with the GP fitting and AF optimizing problems, so improving them by MSR techniques can be a great help.", "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": "cGTfXDM8Hm", "reviewer": "Reviewer_9KgT", "summary": "The paper investigates the perfomance of BO for high-dimensional problems, building upon a recent focus in the literature on working out how to make standard BO frameworks work better in high-dimensions, rather than proposing entirely new frameworks. Lots of empirical evidence is provided to justify the simple solution they propose, i.e. clever lengthscale initialisation + local acquisition optimisation. Along the way, lots of intersting/practically important observations are made.\n\n## update after rebuttal\n\nThanks for the rebuttal. I hope to see all these extra details make it into the final version. I have raised my score.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the positive feedback regarding the clarity of our writing and the potential impact of our simple yet effective solution. We implemented enhancements to clarify our methodology and bolster the impact of our findings.\n\n**Implementation of RAASP**\n\nWe will improve the description of the RAASP implementation in Section 3.2. We will further clarify that MSR uses BoTorch’s RAASP implementation. Since we cannot upload a revised manuscript, we will briefly outline the details here.\n\nWe perturb the top 5% observations using a normal distribution with $\\sigma=0.001$, truncated within $\\mathcal{X}$. For $D\\geq 20$, we also generate samples with only a subset of dimensions perturbed, each with a $\\frac{20}{D}$ probability. Of $4m$ random samples, $2m$ are global samples from a scrambled Sobol sequence, $m$ are local samples around the top 5%, perturbing all dimensions, and $m$ are local samples around the top 5%, perturbing 20 dimensions on average if $D\\geq 20$, or all dimensions if $D<20$. We call the $2m$ local samples the RAASP samples. The starting points for the gradient-based optimization of the acquisition function are drawn from the $4m$ overall samples using Boltzmann sampling [1].\n\n**Details on the AF and its optimization**\n\nWe will expand the methodology section to detail the AF and optimization process without RAASP to improve reproducibility and clarity. We briefly outline the details below.\n\nWe use the LogEI AF [1] for its numerical stability. It is maximized by evaluating it on 512 scrambled Sobol points, then selecting five starting points via Boltzmann sampling for gradient-based optimization using L-BFGS-B, with up to 2000 iterations.\n\n**Magnitude of the contribution**\n\nThis work sheds new light on key challenges and fundamental issues in HDBO, like vanishing gradients in GP model fitting and AFs. These are critical for understanding the limitations of popular algorithms and were not part of prior work. This revision of prior concepts will prompt broader comfort for BO practitioners to use methods such as DSP and MSR.  We expect that our identification of a more straightforward, high-performing implementation - supported by uncovering fundamental principles and efficiencies - enables researchers to investigate novel algorithmic variations and applications. MSR introduces a simple initialization method for GP length scales using MLE that scales with the problem dimensionality, eliminating the need for prior beliefs (Section 4). We further identify the importance of RAASP for avoiding vanishing gradients and reducing variance in length scale estimates (Sections 3.2 and 4). We also show that the MLE estimates have high variance while MAP estimates are biased, relating it to the variance-bias tradeoff (Section 3.3).\n\n**Comparisons with other HDBO algorithms**\n\nWe now compare our method with TuRBO (https://ibb.co/jkx3MXN1), which is competitive on Mopta08 and Lasso-DNA and performs well on Humanoid for the first half of the optimization but stagnates in the second half.\n\n**Discussion of AF optimization budgets**\n\nWe interpret this comment as 1) gradient descent budget, 2) RAASP budget, and 3) multi-start initialization for the AF. We will clarify this in the manuscript by adding these explanations.\n\n1. We capped the number of gradient descent steps to 2000, the BoTorch default budget. We will add to the manuscript that this budget is rarely depleted, as shown in this histogram plot https://ibb.co/bjQQvskj\n2. When augmenting the candidate set of starting points for the AF maximization with RAASP samples, we see fewer gradient descent steps required to optimize the AF. This aligns with previous work claiming that EI attains maxima close to good observations [1,2] and our observation that the starting points for the gradient-based AF maximizer predominantly originate from the RAASP samples. We will add this explanation to the manuscript and pose the question of how much optimization is sufficient in the future work section.\n3. Evaluating the AF on more initial points before selecting starting points can reduce the risk of vanishing gradients. This figure (https://ibb.co/XkL4yCpn) illustrates this by plotting the fraction of points that moved a positive distance after maximizing the AF against the number of initial points. We maximized LogEI on a 50-dimensional GP prior sample ($\\ell=0.2$, isotropic kernel) and repeated the experiment 50 times. With $2^{11}$ random samples, only 2% of final candidates moved a non-zero distance. With $2^{17}$ samples, this increased to 50%. Using RAASP sampling with just 512 raw samples, none of the starting points exhibited vanishing gradients across all 50 repetitions.\n\nThese discussions will corroborate our arguments, and we will add them to the manuscript.\n\n[1] Ament et al. \"Unexpected improvements to expected improvement for bayesian optimization.\" NeurIPS 2023\n\n[2] Hvarfner et al. \"Vanilla Bayesian Optimization Performs Great in High Dimensions.\" ICML 2024", "claims_and_evidence": "The paper claims that the real issue (in contrast to that proposed in recent works) when fitting Gps in high-dimension is to do with vanishing gradients in the MLL. Significant evidence is provided to support this claim, including targeted analysis as well as full BO pipeline runs.", "methods_and_evaluation_criteria": "For what is largely empirical work, the evaluations are somewhat limited with regards to ablations. For example, it would be powerful to see results for differnet acqusition function optimisation budgets (and indeed some actualy details on how this is perfomed....).", "theoretical_claims": "N/A", "experimental_designs_or_analyses": "I think a fair few key pieces of information are missing (or need better signposting) in order to claim reproducibility:\n1. exact details for RASSP implementation (see below, i.e. Gaussian proposal params, and budget e.t.c)\n2. details for the acquisition function used in the experiments\n3. details on how the acquisition function is optimised when not doing RASSP", "relation_to_broader_scientific_literature": "This work builds upon recent work selling the use of sutably modified traditional BO algorithms for high-dimensional BO. This strand of work has oppurtunities for high-impact, especially as the proposed algorithm is so simple/easy to apply in exisiting frameworks/pipelines.\n\nI do worry a little bit about the utility of this work above the recent papers in this field (e.g. the Vanilla BO paper referenced throughout). I agree that the method proposed here seems to be better motivated + more generally applicable, but in absolute terms, the contribution is not huge.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "The paper is very well-written and was a pleasure to read. I feel that I learned a couple of things along the journey, identifying many of the strange behaviours that I have observed empirically when applying BO in practice.", "other_comments_or_suggestions": "It would be nice for a full discusson of RASSP (i.e. what exact params it has for its Gaussian proposal e.t.c), rather than just referencing Botorch (apologies If i missed this).\n\nI know the aim of the paper is to investigate how standard BO algorithms can be adapted to high-dim BO, but it would be a stronger sell if you included some results with more of the popular high-dim BO algorithms, like TURBO e.t.c", "questions_for_authors": "Improved discussion/experimentation around acquisition optimisation budgets and how they interplay with the proposed method would encourage me to raise my score.", "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": "PEnhlKcynr", "reviewer": "Reviewer_MivE", "summary": "This paper considers High-dimensional Bayesian Optimization and thoroughly studies the use of vanilla Gaussian Processes (GP) for that matter. A recent line of work illustrated how simple tweaks can lead to considerable improvements in optimization outcomes, with GPs outperforming complex state-of-the-art baselines. Such tweaks rely on proper lengthscale prior scaling and clever initialization strategies for Marginal log likelihood optimization. The paper thoroughly analyzes the effect that these modifications have on lengthscale optimization, finding that they effectively induce different kinds of search behavior in high dimensions, e.g., local or random. Based on these findings, a simple variant is proposed, and various benchmarks assessing its performances are provided.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their positive feedback, particularly for recognizing the thorough analysis of Gaussian Processes in high-dimensional Bayesian Optimization (HDBO), as well as the clarity and practical insights provided in the paper. We appreciate the reviewer's thoughtful comments and considered each point to enhance the clarity and impact of our manuscript, as detailed below.\n\n**Focus on Gaussian Processes**\n\nWe will clarify in the discussion section that our analysis is specifically tailored to GPs and does not extend to other surrogates. This focus allows us to provide a deep dive into the behavior of GP length scales in high-dimensional spaces. However, studying other surrogates is an important avenue for future work. The vanishing gradients we observe are primarily related to popular stationary kernels for which the covariance of two points vanishes quickly with their distance. Other surrogates, such as random forests or Bayesian neural networks, use different techniques to relate input and output similarity and may show different susceptibility to vanishing gradients.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "- The high-dimensional benchmarks are restricted to vanilla Gaussian Processes and do not involve other surrogates like Bayesian Neural Networks, nor baselines using complex nonlinear embeddings like VAEs. Given the scope of the paper, which focuses on explaining how GP lengthscales behave throughout optimization in high-dimensional spaces, I have no issues with that.\n\n- The high-dimensional datasets employed for evaluation are standard. However, the authors highlight that the high-dimensional benchmarks employed by the community, which they also use in their evaluation, might not be \"that high-dimensional\", as detailed in Appendix C. They argue that \"these benchmarks have a simplifying structure that benefits models with long length scales, posing the risk of algorithms being ‘overfitted’ to these benchmarks.\"", "theoretical_claims": "There is no theoretical claim.", "experimental_designs_or_analyses": "I do not see any particular issue with the experiments and their subsequent analysis.", "relation_to_broader_scientific_literature": "This paper is part of a recent trend, consisting of reverting back to simple methods for high-dimensional black-box optimization. The work it builds on is cited, and the field of high-dimensional BO as a whole is well-described in the related work.", "essential_references_not_discussed": "Not that I am aware of.", "other_strengths_and_weaknesses": "The paper is well-written and well-organized. The technical novelty is light, but that is not the goal here. Instead, the paper does a great job at dissecting the complex behavior of GP lengthscale optimization in high-dimensional BO. As a practitioner myself, I appreciate the insights delivered by this analysis, and I believe the community will highly benefit from them.", "other_comments_or_suggestions": "Fig 5, typo: The behavior of short-length-scale GPs\nreverts to local search (ℓ = 0.05 in the left panel) with RAASP\nsampling and to **local search** without RAASP sampling (ℓ = 0.05\nand ℓ = 0.28 in the right panel).", "questions_for_authors": "I do not have any questions.", "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": "yNQBUAGYFV", "reviewer": "Reviewer_bU1B", "summary": "This paper focuses on High-Dimensional Bayesian Optimization (HDBO). More precisely, it seeks to explain the good performance of vanilla BO for some high-dimensional tasks that has been observed by recent works. To do so, the paper empirically studies the vanishing gradients that appear when fitting a high-dimensional GP or when maximizing an acquisition function and derives multiple insights about the behavior of BO methods in high-dimensional input spaces.\n\nThe paper also proposes ways to mitigate the negative effects of vanishing gradients, such as Random Axis-Aligned Subspace Perturbations (RAASP) sampling or scaled Maximum Likelihood Estimation (MLE), and illustrates their benefits on empirical experiments.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer's positive feedback regarding the insights and empirical observations presented in our paper. We are glad to hear that our work is seen as valuable to the HDBO community, well-written, and enjoyable. \n\n**General claims in the title and abstract**\n\nWe will tone down the \"identifying fundamental challenges\" part of the text and edit the abstract and introduction to communicate the empirical focus of our work to set the right expectations for the reader. \n\n**Sublinear regret bounds**\n\nMost theoretical analyses of Bayesian Optimization, including those that derive sublinear regret bounds, assume that the acquisition function can be maximized exactly [1-4]. RAASP, as a technique to avoid vanishing gradients of the acquisition function, is orthogonal to this. Arguably, it was a long-standing belief that maximizing the acquisition function is simple enough for proofs not to focus on this aspect [5]. Our analysis shows that the error introduced during the acquisition function maximization deserves more attention, but we leave the derivation of sublinear regret bounds under RAASP for future work. \n\n**References**\n\nWe appreciate the suggestion to include more recent references on additivity assumptions. We will incorporate the recommended papers to contextualize our work better.\n\n**Typos and minor issues**\n\nWe will correct the typos and clarify the interchangeable use of 'acquisition function' and 'AF.' We will also fix the undefined reference and ensure consistency throughout the paper.\n\n\n[1] Srinivas, Niranjan, et al. \"Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design.\" Proceedings of the 27th International Conference on Machine Learning.\n\n[2] Wu, Kaiwen, et al. \"The behavior and convergence of local bayesian optimization.\" Advances in neural information processing systems 36 (2023)\n\n[3] Scarlett, Jonathan. \"Tight regret bounds for Bayesian optimization in one dimension.\" International Conference on Machine Learning.\n\n[4] Wang, Ziyu, and Nando de Freitas. \"Theoretical analysis of Bayesian optimisation with unknown Gaussian process hyper-parameters.\" arXiv preprint arXiv:1406.7758 (2014).\n\n[5] Bull, Adam D. \"Convergence rates of efficient global optimization algorithms.\" Journal of Machine Learning Research 12.10 (2011).", "claims_and_evidence": "This paper derives intriguing insights about HDBO and an interesting explanation about why recent HDBO methods perform well in high-dimensional settings.\n\nHowever, I believe that \"understanding HDBO\" or \"identifying fundamental challenges\" are very strong and general claims, however this paper relies exclusively on empirical observations to support them. The results clearly point in an interesting direction, but more theoretical research is needed for, e.g., derive an upper bound on the magnitude of the gradients for the MLE and the acquisition function regardless of the objective function or quantify the impact of RAASP sampling on the magnitude of the gradients / on the regret of the algorithm.", "methods_and_evaluation_criteria": "The paper identifies vanishing gradients as the main reason why vanilla BO does not perform well in high-dimensional input spaces and proposes a method that seems to mitigate this issue.\n\nThe method is compared against common strategies to learn the model hyperparameters and optimize high-dimensional objective functions, on high-dimensional problems and using multiple evaluation criteria.\n\nOverall, this makes a lot of sense for the problem at hand.", "theoretical_claims": "Excluding the (overly) general claims made in the title and the abstract of the paper, this work does not make any theoretical claim.", "experimental_designs_or_analyses": "The paper provides a lot of empirical observations and visualizations to support their conclusions on HDBO. All of them are insightful and sound.", "relation_to_broader_scientific_literature": "As I mentioned before, the paper provides very relevant insights on HDBO. I think this exploratory research is valuable to the HDBO community and should spark theoretical research efforts to better study the phenomena highlighted in this paper.", "essential_references_not_discussed": "I would have liked to see more recent papers on the additivity assumptions (e.g., [1, 2]), especially [2] which is the most recent and general work on the subject.\n\nHowever, these are only useful in the contextualization of the paper. As far as I know, the paper properly discusses prior work.\n\n**References**\n\n[1] Hoang, T. N., Hoang, Q. M., Ouyang, R., & Low, K. H. (2018, April). Decentralized high-dimensional Bayesian optimization with factor graphs. In Proceedings of the AAAI Conference on Artificial Intelligence (Vol. 32, No. 1).\n\n[2] Bardou, A., Thiran, P., & Begin, T. (2024). Relaxing the Additivity Constraints in Decentralized No-Regret High-Dimensional Bayesian Optimization. In The Twelfth International Conference on Learning Representations.", "other_strengths_and_weaknesses": "I would like to point out that the paper is well-written, I actually enjoyed reading it.\n\nI have no other major strength / weakness to report.", "other_comments_or_suggestions": "As mentioned before, I would suggest to rephrase some claims to point out the empirical / exploratory nature of the research conducted in this paper.\n\nAdditionnally, here is a list of typos:\n* in the abstract, \"investigates the 'why'.'\",\n* \"acquisition function\" and its acronym \"AF\" are used interchangeably throughout the paper,\n* on page 5, \"more spread more evenly\",\n* Undefined reference to section on page 16.", "questions_for_authors": "Here is a question to spark the discussion with the authors:\n\n1) Have you thought about the implications of vanishing gradients (and the use of mitigating techniques like RAASP sampling) on the asymptotic performance of HDBO algorithms? Do you think that it is still possible to derive a sublinear regret bound under RAASP?", "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:

Meta-review: This paper investigates the problems of Bayesian optimization in high-dimensional space, discussing vanishing gradient and the length-scale prior.  It also proposes mitigations for the problems.

Although there is no theoretical contribution, and the experiments are not comprehensive, all reviewers found that the paper is well-written, well-organized, and provides many impactful insights, which facilitate new research in the ML community.  Minor concerns have been addressed in the rebuttal, and expected to be reflected to the camera-ready version.
Decision: Accept (poster)