High-dimensional Bayesian Optimization via Covariance Matrix Adaptation Strategy
Authors: Lam Ngo
Huong Ha
Jeffrey Chan
Vu Nguyen
Hongyu Zhang
Venue: Accepted by TMLR ()

Abstract: Bayesian Optimization (BO) is an effective method for finding the global optimum of expensive black-box functions. However, it is well known that applying BO to high-dimensional optimization problems is challenging. To address this issue, a promising solution is to use a local search strategy that partitions the search domain into local regions with high likelihood of containing the global optimum, and then use BO to optimize the objective function within these regions. In this paper, we propose a novel technique for defining the local regions using the Covariance Matrix Adaptation (CMA) strategy. Specifically, we use CMA to learn a search distribution that can estimate the probabilities of data points being the global optimum of the objective function. Based on this search distribution, we then define the local regions consisting of data points with high probabilities of being the global optimum. Our approach serves as a meta-algorithm as it can incorporate existing black-box BO optimizers, such as BO, TuRBO, and BAxUS, to find the global optimum of the objective function within our derived local regions. We evaluate our proposed method on various benchmark synthetic and real-world problems. The results demonstrate that our method outperforms existing state-of-the-art techniques.

Published in Transactions on Machine Learning Research (01/2024)
High-dimensions Bayesian Optimization via Covariance Matrix Adaptation Strategy
Lam Ngo
63962378 wwwt:eeeeeetttt..eeeee...
RMIT University, Australia
Huong Ha
bhong.hharrmit.com wew.:mmtt..ddddooo
RMIT University, Australia
Jeffrey Chan
jeffrey.chon&rmit.cc
RMIT University, Australia
Vu Nguyen
eutngn&@amaaznn.com
Amazon, Australia
Hongyu Zhang
hychang@oqu.edd edu.co
Chongging University, China
Reviewed on OpenReview: https: // openreview. net/ forum? id= eTger?gPuU
Abstract
Bayesian Optimization (BO) is an effective method for finding the global optimum of expen- sive black-box functions. However, it is well known that applying BO to high-dimensional optimization problems is challenging. To address this issue, a promising solution to use a local search strategy that partitions the search domain into local regions with high likeli- hood of containing the global optimum, and then use BO optimize the objective function within these regions. In this paper, we propose novel technique for defining the local regions using the Covariance Matrix Adaptation (CMA) strategy. Specifically, to learn search distribution that can estimate the probabilities of data points being the global optimun of the objective function. Based on this search distribution, we then define the local regions consisting of data points with high probabilitics of being the global op- limum. Our approach serves as mela-alyorithm as it can incorporate existing black-box Bо optimizers, such as BO, TuRBD (Eriksson al., 2019), and BAxUS (Papenmeier al., 2022), to find the global optimum of the objective function within our derived local regions. We evaluate proposed method various benchmark synthetic and real-world problems. The results demonstrate that our method outperforms existing state-of-the-art techniques.
use CMA
1 Introduction
Optimizing expensive black-box functions is important task that has various applications in machine learning, data science, and operational research. Bayesian Optimization (BO) (Jones et al., 1998; Brochu et al., 2010; Shahriari et al., 2016; Binois & Wycoff, 2022; Garnett, 2023) is a powerful approach to tackle this challenging problem in an efficient manner. It has been successfully applied in a wide range of applications, including but not limited to hyperparameter tuning of machine learning models (Snoek et al., 2012; Turner et al., 2020), neural architecture scarch (Jenatton al., 2017: Kandasamy et al., 2018b), material design (Heráández-obatte al., 2017), robotics (Calandra et al., 2016), and reinforcement learning (Brochu al., 2010; Harker-Holler et al., 2022).

---

Published in Transactions on Machine Learning Research (01/2024)
BO operates in an iterative fashion by repeatedly training surrogate model based the observed data and using an acquisition function to suggest promising data points for evaluation (Garnett, 2023). This method is inspired by Bayes' theorem, which aims to improve the prior belief about the objective function by incorporating observed data to obtain a posterior with better information. In this way, BO selects the next data points by considering previous information and maximizes the knowledge gained about the objective function with new observations, making it sample-cfficient in finding the objective function's 's global optimum.
Despite being a powerful optimization method, BO still suffers from various problems, including the curse of dimensionality issue, ie., it often performs poorly when applied to high-dimensional problems (Rana et al., 2017; Eriksson et al., 2019; Binois & Wyeoff, 2022; Papenmeier et al., 2022). One reason for this is that, as the search domain grows in dimensionality, more local optima may appear, making it difficult for the algorithm to find the global optimum. Additionally, a larger search space also implies more regions with high uncertainty, which can potentially cause the acquisition function to overemphasize exploration and fail to exploit potential regions within a fixed budget. Morcover, in high-dimensional spaces, the objective function is typically heterogeneous, making it difficult to fit a global surrogate model across the entire domain.
There have been various works attempting to make BO work well for high-dimensional optimization problems. One of the most promising approaches, which have demonstrated significant success, is the use of local search strategy that partitions the search domain into promising local regions where the optimizat ion process is performed within (Munos, 2011; Wang et al., 2014; Eriksson et al., 2019; Wang et al., 2020; Wan et al., 2021). These works, however, have certain limitations. For instance, the works in Munos (2011); Wang et al. (2014); Eriksson et al. (2019) employ a search space partition technique with fixed parameters that maybe difficult to optimally specify in advance, and thus, may not provide adequate flexibility for various problems (Wang et al., 2020). The work in Wang et al. (2020) learns promising local regions by partitioning the search space into non-linear boundary regions using an unsupervised classification algorithm, but there is no guarantee that these local regions can be learned accurately with a limited amount of training data.
In this paper, we follow the aforementioned local search approach to tackle the high-dinensional optimization problem. We propose to use the Covariance Malriz ddpptation (CMA) strategy to systemat ically define the local regions be used within this local search approach. CMA is a technique developed in the Evolutionary Algorithm literature (Hansen & Ostermeier, 2001), aiming to estimate the probability distribution of data points in the search domain being the global optimum of the objective function (i.e., search distribution). Typically, CMA is combined with Evolutionary Strategy (ES) techniques like CMA-ES (Hansen & Ostermeier, 2001). These techniques leverage the search distribution derived from CMA to guide the search for the global optimum toward promising regions in the search domain, i.e., regions that highly likely contain the global optimum. It has been shown that CMA-based ES techniques, such as CMA-ES, perform very well in finding the global optima of high-dimensional optimization problems, demonstrating CMA's effectiveness in identi- fying promising regions that likely contain the global optima of high-dimensional optimization problems. A drawback of these techniques is that they generally require a large number of function evaluations (Loshchilov & Hutter, 2016; Nomura et al., 2021).
Inspired by the effectiveness of CMA in working with higl-dimensional optimization problems, we propose to incorporate CMA into BO methods by using it to define the local regions in the BO local search approach. In particular, we define the local regions as the regions with the highest probabilitics of containing the global optimum based on CMA's 's search distribution. Subsequently, we can use an existing BO optimizer, e.g., BO, TuRBO (Friksson et al., 2019), BAXUS (Papenmeier et al., 2022), within these local regions to find the global optimum of the objective function. By leveraging information from CMA's search distribution, BO hothods can focus the search within the local regions that have high likelihood of containing the global optimum of the optimization problem. CMA-based BO methods are therefore expected to work well with high-dimensional optimization problems, due to CMA's effectiveness, whilst preserving data-efficiency, a property lacking in CMA-based ES techniques. We derive the CMA-BO, CMA-TURBO and CMA-BAXUS algorithms corresponding to the cases when we incorporate CMA with the optimizers BO, TuRBO, BAXUS, respectively. Our experimental results on various synthetic and real-world benchmark problems confirm that our proposed approach helps existing BO methods to work better for high-dimensional optimization problems whilst being data-efficient.
In summary, our contributions are as follows:
2

---

Published in Transactions on Machine Learning Research (01/2024)
Proposing novel meta-olgoritlmm using the local search approach and the CMA strategy to enhance the performance of existing BO methods for high-dimensional optimization problems;
Deriving the CMA-based BO algorithus corresponding to the cases when incorporate the pro- posed meta-algorithm with the state-of-the-ar BO methods (e.g., BO, TuRBO, BAXUS);
Conducting comprehensive evaluation on various high-dimensional synthetic and real-world bench- mark problems and demonstrating that our proposed CMA-based mela-algorithm outperforms ex- isting state-of-the-art methods.
The implementation of our method is available at httpppr//ggttto.cco///////////ag
2
Background
In this section, we present the fundamental background BO. Then we revisit two state-of-the-ar BO methods (TuRBO, BAxUS) that we will incorporate into our CMA-basod meta-algorithm.
2.1 Bayesian Optimization
Bayesian optimization (BO) is powerful optimization method to find the global optimum of an expensive black-box objective function by sequential queries (Jones et al., 1998; Brochu et al., 2010; Shahriari et al., 2016; Husain et al., 2023; Garnelt, 2023). Let us consider the minimization problem: given an unknown objective function f X + where c IRd is a compact space, the goal of BO is to find a global optimmm xx of the objective function f:
x² E'g m fxx).
(1)
BO solves an optimization problem in an iterative manner. First, the objective function f is approximated by a surrogate model trained with the current observed dataset Do = {æ,, Mo with yi = f(;)) + 6 and e; ~ ((0, o²) being the corrupted noise. Then an acquisition function a : X + is constructed from the surrogate model to assign scores all data points in the domain X based on their potential to improve the optimization process. The next data point to be evaluated, denoted as Fnext; is selected as the maximizer of the acquisition function. Subscquently, the objective function is evaluated Enext and the new observed data (Xnext, Whext); with Inext = ((xnex + 8 and ε (O o²), is added to the current observed dataset. The process is conducted iteratively until a pre-defined budget is exhausted, and then BO returns the best value found from the observed datasct as an estimate of the global optimum x².
There are different choices for the surrogate models to be used in BO, including Gaussian Process (GP) (Ras- mussen & Williams, 2006), Tree-structured Parzen Estimator (TPE) (Bergstra et al., 2011), and neural networks (Springenberg et al., 2016). In work, we focus the GP surrogate model, which is one of the most popular surrogate models in BO. GP is probabilistic model that can provide both scalar pre- dictions for the objective function values and their associated uncertainty. A GP is completely specified by its mean function μ(x) and covariance function (kernel) k(æ, x') (Rasmussen & Williams, 2006). While the mean function indicates the most probable values for the objective function values, the covariance function captures the properties of the objective function, e.g. its smoothness.
There is also variety of common acquisition functions. Examples include Expected Improvement (EI) (Mockus et al., 1978), Probability of Improvement (PI) (Kushner, 1964), Upper Confidence Bound (UCB) (Srinivas et al., 2010), Thompson Sampling (TS) (Thompson, 1933), and Knowledge Gradient (KG) (Frazier et al., 2009). While each of these acquisition functions has its own strengths and weak- nesses, they all aim to balance between exploration and exploitation. Exploration encourages the algorithm to look for promising values in highly uncertain locations, whilst exploitation favors refining the knowledge around the currently optimal locations. In this work, we use Thompson Sampling (TS) (Thompson, 1933) as the acquisition function following Eriksson et al. (2019); Papenmeier et al. (2022). In the next section, will describe in detail the TS acquisiiion function which will be later used in our method.
3

---

Published in Transactions on Machine Learning Research (01/2024)
2.2 The Thompson Sampling Acquisition Function
The Thompson Sampling (TS) acquisition function (Thompson, 1933), commonly used in BO research (Kan- dasamy et al., 2018a; Eriksson et al., 2019; Papenmeier et al., 2022), follows a stochastic policy. The main idea is to sample a random rcalization (i.c., sample path) of the objective function from its posterior distri- bution, then optimize this sample path to find the next data point to be evaluated. Specifically, at iteration 1, given the observed dataset. Dr, the TS acquisition function a a= (x;; D.) can be defined as,
a (x; D,) (a) where (xr) ~ P.))
(2)
and pt)) denotes the posterior distributton of the objective function ʃ given the observed data D,. GP is used as the surrogate model for the objective function, then (x) (P(u(z)) k(z, DD)) For minimization problem as defined in Eq. (1), the TS process is then to minimize the acquisition function TS (x;
(s:; D:) to find the next data point Fe+1 to be evaluated,
F+ = arg mingex a (x;; D.).
The TS acquisition function selects the data point to be evaluated by drawing a random realization of the objective function f from its posterior distribution, therefore, it encourages exploitation of regions with higher probabilities of being optimal, while allowing for exploration of other regions, owing to its random nature. It thus satisfies the desirable property of balancing between exploitation and exploration, which is a requirement for any acquisition function BO.
2.3 TurBO
TuRBO (Eriksson et al., 2019) a state-of-the-art BÔ method, which proposes to use the local search strategy to solve the high-dimensional optimization problem. In TuRBO, the global search space partitioned into smaller domains, called trust. regions (TR) (Yuan, 2000), within which the surrogale model is believed to accurately model the objective function. TuRBO uses GP the surrogate model and trains the GP using all the previous observed data then selects the next observed data point by optimizing the acquisition function locally within the TR. This local strategy makes TuRBO more efficient than standard BO methods in high-dimensional problems, as it only requires modeling of local surrogate models, abandoning the global surrogate model used in standard BO methods. The local surrogate models of TuRBO do not suffer from the heterogeneity of the objective function, as they only need to accurately caplure the objective function within the TR. Moreover, as the TR reduces the regions with large uncertainty, TuRBO mitigates the expllxxpoon issue commonly encountered by standard BO methods when solving high-dimensional optimization problems. The detailed description of TuRBO can be found in Appendix Section A.
2.4 BAXUS
BAXUS (Papenmeier et al., 2022) is state-of-tho-art BO method that tackles the high-dimensional optimiza- tion problem using subspace embedding approach (Wang et al., 2016; Nayebi et al., 2019; Letham et al., 2020). The main idea is to assume the existence of low-dimensional subspace (active subspace) Z c Rd (de < d), function : Z + R and a projection matrix T E Rd xd such that Væ E XTT((T)) = f(x). This property enables the optimization process to be conducted in the active subspace, which has a lower dimension compared to the original high-dimensional space. In practice, the effective dimensionality de is generally unknown, therefore, existing approaches (Wang al., 2016; Nayebi et al., 2019; Letham et al., 2020) randomly choose subspace V C Rd (target space) to project the original space into, and perform optimization within this chosen subspace. The target dimension dy must be chosen such that the probability of the target space containing the global optimum is high. In practico, choosing appropriate value of the larget dimension dy is challenging as a small value of dy does not guarantee that the target space contains the global optimum whilst a large value of dy could be subject to the curse dimensionality issue. BAXUS proposes an adaptive strategy to gradually increase the target dimension during the optimization process, guarantecing a higher probability, compared to existing approaches, that its embedding contains the global optimum of the objective function. The detailed description of BAXUS is in Appendix Section A.2.
4

---

Published in Transactions on Machine Learning Research (01/2024)
3 Related Work
Various research works have been conducted to address the challenge of applying BO to high-dimensional optimization problems. One approach is to exploit the additive structure of the objective functions, then construct and combine a large number of Gaussian Processes (GPs) to approxinate the objective fune- tion (Kandasamy et al., 2015; Gardner et al., 2017). Some works suggest replacing GPs with surrogate models that might scale better with high-dimensional data such as random forcsts (Hutter et al., 2011), deep neural networks (Snoek et al., 2015), or Bayesian neural networks (Springenberg el al., 2016). Oh et al. (2018) propose BOCK whose main idea is to employ cylindrical transformation to transform the geometry of the search space, thus mitigating the over-exploration issue of BO in high-dimensional optimization prob- lems. Other methods, such as the work by Garnett et al. (2014), REMBO (Wang et al., 2016), HeSBO (Nayebi et al., 2019), ALEBO (Letham et al., 2020), and BAxUS (Papenmeier et al., 2022) propose mapping the orig- inal high-dimensional space into a low-dimensional space and then conducting optimization within this low-dimensional space. More recently, Song et al. (2022) propose MCTS-VS, a meta-algorithm that employs Monte Carlo tree search (MCTS) to construct a low-dimensional subspace and then use a BO optimizer (e.g.; BO, TuRBO) to optimize the objective function within this subspace.
Another approach that has recently attracted lot of attention is the use of a local search strategy that parti- tions the search domain into promising local regions where the optimization can be performed within (Munos, 2011; Wang et al., 2014; Eriksson et al., 2019; Wang et al., 2020; Turner et al., 2020; Fröhlich et al., 2021; Wan et al., 2021; 2022). Notable methods in this direction include TuRBO (Eriksson et al., 2019), whose main idea is to construct local regions as hyper-rectangles centered around the best-found values, and dynami- cally expand or shrink these regions based on the function values. Wan et al. (2021) extend this idea for the categorical and mixod scarch space scttings. Another method proposod in Wang ct al. (2020), namely LA-MCTS, is meta-level algorithm that uses an unsupervised K-mean algorithm to classify the search space into good and bad local regions, within which a BO optimizer such as BO TuRBO can be employed.
Other related works, including the rescarch by Müller et al. (2021); Nguyen et al. (2022), also employ a local search strategy. However, different. compared to the local search approach of partitioning the search space into local regions, this approach aims to incorporate gradient information into the BO process to enhance its effectiveness. In particular, the work in Müller et al. (2021) develops a probabilistic model that can incorporate gradient information, and selects data points for evaluation as those that maximize the gradient. Building on this work, Nguyen et al. (2022) propose a. new acquisition function designed to maximize the probability of gradient descent, thereby enabling BO to rapidly converge toward the local optimumm region. The local stralegies in these hod are different from the space partitioning mechanism used in TuRBO and LA-MCTS These methods aim to perform local optimization using the current solution and seek to descend the objective function values via the gradient information computed from objective function values in nearby regions. These methods complement our proposed CMA-based meta-algorithm, as they can be used as optimizers within our approach.
Evolutionary Algorithms (EAs) represent a widely-used family of algorithms for optimizing high-dimensional black-box functions. Among these, CMA-ES (Hansen & Ostermeier, 2001) is well-known for its impressive performance in finding global optimum of the objective function. In this paper, we propose novel meta- algorithm that leverages the CMA technique of CMA-ES to define the local regions. There are several EA methods that also combine the CMA technique with the GP surrogate model to enhance the optimization performance. GPOP (Buche et al., 2005) is an EA method that uses CMA-ES to optimize a merit. function defined by the predicted mcan and standard deviation function of a trained GP. DTS-CMAES (Bajer et al., 2019), another EA method which has been shown to outperform GPOP, constructs a doubly-trained GP inside CMA-ES to select the data points for forming the mean vector and covariance matrix of CMA. Closely related to our approach, BADS (Acerbi & Ma, 2017) is a global optimization method that adopts the mesh adaptive direct search framework (Audet & Dennis, 2006) which uses CMA-ES as the search oracle and a GP-based method to find the global optimum of the objective function.
In our experiments, we compare the CMA-basod BO methods (crcated by combining our proposed CMA- based meta-algorithm with the BO optimizers BO, TuRBO, BAxUS) with a comprehensive list of related meth- ods: the standard BO method, TuRBO, BAXUS, LA-MCTS, MCTS-VS, CMA-ES, DTS-CMAES and BADS.
5

---

Published in Transactions on Machine Learning Research (01/2024)
4 High-dimensional Bayesian Optimization via Covariance Matrix Adaptation
In this section, we first discuss the key idcas of the CMA strategy (Section 4.1), then we propose the CMA- based meta-algorithh (Section 4.2). Finally, we derive the CMA-based BO algorithms (CMA-BO, CMA-TuRBO, CMA-BAxUS) corresponding to the cases when we integrate the CMA-based meta-algorithm with the BO optimizer B0, TuRBO, and BAxUS (Sections 4.2.1, 4.2.2, and 4.2.3).
4.1 The Covariance Matrix Adaptation Strategy
The CMA strategy was initially developed in the Evolutionary Algorithm literature, particularly in the CMA-ES method (Hansen & Ostermeier, 2001). Its main idea is based on stochastic search, which maintains search distribution, p(æ:), to estimate the probabilities of data points in the search domain being the global optimum of the objective function. First, a search distribution (x) is initialized. Then a population of data is sampled from (z)), and the search distribution is updated based the objective function values of these data points. This process is conducted iteratively until the algorithm converges (Hansen & Ostermeier, 2001; Abdolmaleki et al., 2017). CMA provides well-established formulas (details below) for the updates of the search distribution, enabling the stochastic search algorithm to eventually allocate the highest probability to the global optimum. It has been shown that C2MAbased ES techniques, such as CMA-ES, perform very well in finding the global optima of high-dimensional optimization problems (Toshchilov Hutter, 2016; Eriksson et al., 2019; Letham al., 2020; Nomura al., 2021).
In practice, the most popular choice for the search distribution used in CMA is the multivariate normal distribution (Hansen & Ostermeier, 2001; Hansen, 2016). Consequently, the focus of CMA is on updating the two principal moments of this distribution: the mean vector m and the covariance matrix Σ. In CMA, the covariance matrix is normally decomposed as Σ = o²C where a > is the overall standard deviation (step size) of the search distribution, thus the goal of CMA is then to update m, a, and C. At iteration t, given A observed data points {a]²) .3
(t) ) sampled from the search distribution / (o)t-))) ((-1) of
the previous iteration t - 1, the mean vector m covariance matrix ((O), and step size g( of the search distribution at iteration & are updated as follows (Hansen & Ostermeier, 2001; Hansen, 2016),
m (t) m (t-1) + Cm Σ a
Wi (2 m !
( (1 - -1 -c,c (1--)) cc.)
Cp Exm
wi (0t1)) m (((-)) (af2 m m(--))) (OT, (4)
+ (o(t-1)/3
i=1
((() = (t1))
where
((t) EL (m) (o) denotes the evolution path which quantifies the overall movement of the search distribution, i.e., the movement of the mean vector across iterations,
Tyi denotes the i-th best candidate out of X data points (ax 1 based on their noisy function values, ie., their corresponding noisy function values satisfy: B.d S v < < v
E X is a hypeparameter denoling the number of data points selected to update the search distri- bution; by default, is usually set as [A/2].
f denotes the weight coefficients associated with the data points (a such that w1 Z W2 2 w/ > 0 > Wp+1 2 ... 2 wx Σ 1, and Σ w, R
Cm, C1 and e₁ denote the learning rates at which the search distribution changes, where larger rates result in faster change and smaller rates reduce the adaptation rate of search distribution,
(t)) denotes a modification rule for the step size, depending on the evolution path (t),
Detailed information for the suggested settings of these hyperparameters can be found in Appendix Section A.3. From Eq. (4), it can be seen that in CMA, the mean vector m!!) is updated based the previous

---

Published in Transactions on Machine Learning Research (01/2024)
Data sampled BO optimizer
CMA Update
Ellipsoid mean
Global optimum
(1) Intialization
(2) Data Sampling
(3) Local Region Update
(n) Convergence
Figure 1: Illustration of the proposed CMA-based meta-algorithm. In step (1), a hyper-ellipsoid local region is initialized. In step (2), a BO optimizer (c.g., BO, TuRBO, BAxUS) is used in this local region collect a population of candidates (data points be evaluated). In step (3); the local region is updated using the CMA technique. The process is conducted iteratively until the evaluation budget is depleted.
mean vector (t-1) and the highest-ranking observed data points (--2)=- The covariance matrix CO) is updated based the previous covariance matrix C(t-1), all the observed data points {- ( 1 and the evolution path p of the search distribution from previous ilerations. The step size o! is updated based on 8, which depends on the overall movement of the scarch distribution (the evolution path ())) If the evolution path is short, e.g., when the vectors Am = m (f-1) in consecutive iterations cancel each other out, the step size a is decreased, as the search distribution is likely to start converging toward a solution. On the contrary, when the evolution path is long, e.g., the vectors Am are in the same direction, the step size is increased, as the search distribution is likely be far away from the true one.
Theoretically, it has been shown that the CMA strategy can be interpreted as natural gradient learning method that updates the parameters (mean, covariance matrix, step size) of the search distribution p(æ) to minimize the expected function value B[f(z)] under this distribution (Akimoto et al., 2010; Nomura et al., 2021). With the updates in Eq. (4), CMA tends to maximize the probability of generating successful data points (a ) (c.g., data points with lower function values for minimization problem) in the subsequent iterations (Hansen & Auger, 2011). Empirically, as discussed at the beginning of this section, CMA-based ES techniques like CMA-ES perform very well in finding the global optimum of high-dimensional optimization problems (Loshchilov & Hutter, 2016; Nomura et al., 2021), demonstrating the effectiveness of CMA in deriving search distributions that can estimate the probabilities of data points being the global optimum of the objective function.
4.2 The CMA-based Meta-algoritha
In this section, we present our proposed CMA-based meta-algorithm. We first discuss the overall process of this meta-algorithm, then we derive three CMA-based BO methods where integrate the proposed meta-algorithm with the three state-of-the-art BO optimizers (BO, TuRBO, BAXUS)
Overall Process. We ralustrate the CMA-based meta-algorithm in Fig. and the pseudocode Algo- rithm 1. First, an initial search distribution (m ((O)) is set (line 6), and local region S(O) is computed based on this search distribution and the chosen BO optimizer bo_apt (line 9). Then the BO optimizer bo_ vo_opt is used within the local region S(0) to suggest data points Dx = {æ,} to be evaluated (lines 11-14). The search distribution is then updated based on these À observed data points Dx (line 15). The process is conducted iteratively until the evaluation budget depleted, and the algorithm terminates. Note that restart strategy is also included to restart the algorithm when the current optimization process is stuck at local minimum (line 17).
Local Region Formulation. We first define the base local region S, for the CMA-basod meta-algorithm. Depending on the employed BO optimizer, we will further derive the local region S corresponding that particular BO optimizer in the later sections (Sections 4.2.1, 4.2.2, and 4.2.3). As discussed in Section

---

Published in Transactions on Machine Learning Research (01/2024)
Algorithm 1 The CMA-based meta-algorithm.
1: Input: Objective function f()), search domain [1, u]] maximum number of function evaluations N;
number of initial points no, BO optimizer bo_opt Output: The optimum x²
3: Set L 0, T 4 [(N- nn)/Aj, global dataset D -- 0, local dataset Ω < 0, population size À
4: while t S T do
5: Sample no initial data points Do
b Latin hypercube
Set the initial search distribution , (a(0)) based on Do
Sec. A.5
7: Update global dataset D + D U Do- Reset local dataset 12 4- Do, restart + False
8: while t T and not restart do Compute the local region ((t) from N(ml)) (o(0) depending on bo bo_ Eqs. (6))
10:
Initialize a dataset to collect À observed data points Dx + 0 for i=1:A do
11:
12:
Apply BO optimizer bo_opt to propose an observed data {ar yi from a pool of data points Update D2 + DA U fz;, y
13:
14:
end for
15:
Update (m++)) c(++), (((t)) + CMA((m!), "0),,()), Dx) Update DUD, S Q DUD ++t+1
Eq. (4)
16:
17:
Update restart t True if stopping criteria satisfied end while
18:
19:
end while
20:
Return x² = arg mn from global dataset D = {(x;,, 1:)
4.1, the search distribution by CMA can assign higher probabilities to more promising data points; there- fore, can define the local regions as the regions containing data points with high probability values from this search distribution. As CMA's search distribution is multivariate normal distribution, we pro- pose defining the local region as the a-level confidence hyper-ellipsoid centered at the mean vector of this search distribution, containing a percent of the population of data points that follows this search distribu- tion. Specifically, at iteration t, given the multivariate normal search distribution ((m((11), D)), where ((-1)) (t) ((t-1), obtained in the previous iteration, the base hyper-ellipsoid local region sg() can be computed as,
56 = fo = (() S 2 x
where (x) -((--m(-1)) (2(-1)))- ( (xt m(t-1)) is the Mahalanobis distance (Mahalanobis, 1936) from æ to the search distribution (t- ()-))) and xi=a. is the Chi-squared 1 1 critical value with d degree of freedom. In our proposed CMA-based meta-algorithm, we set a to be 99.73% corresponding to the 3-sigma rule that is commonly used in practice. With this setting the selected observed data in each iteration will always fall within the three standard deviations of the mean vector of the search distribution.
Local Optimization. In each iteration 4, given the multivariate normal search distribution (t,, obtained in the previous iteration, we first sample a pool of data points that follow this search distribution and are within the local region S(t), Then, we use the employed BO optimizer, bo bo_opt, to select À data points from this pool of data points. The rationale behind this step that in the CMA strategy, the search distribution (tm((t))), E((-)) provides estimates of the probabilities of data points in the search domain being the global optimum of the objective function. By sampling data points following this scarch distribution, we can have pool of data points having high probabilities of being the global optimum. By using BO to select data points from this pool, we have a higher probability of select- ing the better data points that are close to the global optimum. Finally, after obtaining X observed data points, we update the CMA's scarch distribution following Eq. (4). is worth noting that, in our proposod approach, in each iteration, we sample and evaluate À data points rather than just one as in standard BO methods, iee, in our algorithm, one iteration is equal to À iterations in standard BO methods.
8

---

Published in Transactions on Machine Learning Research (01/2024)
Restart Strategy. A local search strategy is typically biased toward the starting point, and the optimiza- tion process can be trapped in local minima (Friksson et al., 2019). To enable global optimization for the CMA-based meta-algorithm, we use the CMA's restart strategy: a new local search will be initialized when the current one is stuck at a local minimum (Auger & Hansen, 2005; Hansen, 2016). The conditions for a restart CMA normally involve checking if the objective function values are flat for a number of iterations or if some numerical indicators (e.g., condition number) of the search distribution are violated. Further- more, since some BO optimizers (e.g,, TuRBO, BAzUs) have their own restart strategies, we also incorporate these restart strategies when applying our proposed CMA-based meta-algorithm to the corresponding BO optimizers (detailed information in the subsequent sections).
4.2.1 CMA-BO: CMA-based Meta-algorithm with Standard BO
In this section, we describe CMA-BO, the corresponding CMA-based BO method obtained when incorporating the proposed CMA-based meta-algorithm with the standard BO optimizer. Note that from this section, we remove the superscript denoting the iteration index for brevity.
Local Region Formulation. For CMA-BO, we define the local region 8 equal to the base local region S, described in Eq. (5), i.e., the local region is the a-level confidence hyper-ellipsoid of the search distribution ((m, E). The value a is also set at 99.73%, corresponding to the 3-sigma rulc.
Local Optimization. Following the base algorithm described in Section 4.2, in each iteration, we first sample pool of data points that (1) follow the previous search distribution N(m, E), and, (2) are within the local region S. Then we sequentially apply BO with the TS acquisition function select the best X data points from this pool. Note that when training the GP, in Friksson et al. (2019), we also use all the observed data in all the previous iterations. The pseudocode of the local optimization step in CMA-BO is in Appendix Section A.4, Algorithm 2.
Restart Strategy. CMA-BO has the same restart strategy with CMA as described in the base algorithm.
4.2.2 CMA-TURBO, CMA-based Meta-algorithm with TuRBO
We derive CMA-TURBO the CMA-based BO method obtained when incorporating our CMA-based meta- algorithm with the TuRBO optimizer. The challenge here is that TuRBO has its own local region adaptation mechanism to shrink or expand. Furthermore, the shape of the local regions of TuRBO is hyper-roctangle which is very different from our CMA local regions' shape which is hyper-ellipsoid.
Local Region Formulation. For CMA-TURBO, we incorporate TuRBO's local region adaptation mechanism, which is based on the success and failure state of the optimization process, with the local region strategy defined by CMA. Specifically, with the search distribution N(m, Σ), we define the local region SCMA-TURBO as the hyper-ellipsoid with: (1) the center being at the mean vector m, (2) the radii (lengths of the semi-axes of the hyper-cllipsoid) computed based on the covariance matrix Σ and scaled with factor L that is based the success and failure state of the optimization, similar to the local region adaptation mechanism in TuRBO. In particular, L is initially set as 0.8, and after Tsuce consecutive success, L is doubled, while it is halved when the optimization fails to progress after Tiail consecutive times. Therefore, compared to the local regions defined by CMA-BO. the local regions of CHA-TuRBO are scaled based on the historical optimization success record as in TuRBO. With a scale factor L, the local regions of CMA-TuRBO are defined as follows,
SCnA- TURBO = {x ACMA-TuBBO (x) S
2 xi-ad))
(6)
where (us))))))) = V (x - m) TuRuA-TRRO (a m) is the Mahalanobis distance from x to the scaled search distribution N(m, ECHA-TUROO) and ECCA-TARO L.S.. The derivation of this scaled covariance matrix is as follows. By definition, the radii of the hyper-ellipsoid constructed by Σ is r = diag(A (A 1/2), where A is the diagonal matrix whose diagonal elements are the cigenvalues of Σ. Note that, using the cigendecomposition. we have that, Σ- = UAU-1 with U being the eigenvector matrix. Thus, when scaling the radii of this hyper-ellipsoid by T., i.e., TCMA-TURBO = Lr, the covariance matrix Ecu-TuRao becomes U(T2A)U-1 = I.2.. Finally, similar to TuRBO, we also set upper and lower bounds for L, i.e., L cannot exceed a threshold Lmax and when L becomes smaller than a threshold Lmiu, the algorithm restarts.
9

---

Published in Transactions on Machine Learning Research (01/2024)
Local Optimization. After obtaining the scaled search distribution (mm, Ecu-runso) and the local region SCNA-TARBD; the optimization process is conducted similarly as in the base algorithm described in Section 4.2. Specifically, we first sample a pool of data points from the search distribution ((m. ECHA-TURRO), and then apply BO with the TS acquisition function to select data points from this pool. Besides, when training the GP, as with Eriksson et al. (2019), we use all the observed data points so far. The pseudocode of the local optimization step in CMA-TuRBO is in Appendix Section A.4, Algorithm 3.
Restart Strategy. Apart from the restart stralegy of CMA, we also employ the restart strategy of TuRBO, i.e., when L shrinks below a minimum threshold Limi; we terminate the CMA local region SunA-TuRa0 and restart it at a new location randomly.
4.2.3 CMA-BAXUS: CMA-based Meta-algorithir with BAXUS
We present CMA-BAxUS, the method resulted when incorporating our proposed CMA-based meta-algorithm with the BAXUS optimizer. The difficulties in deriving CMA-BAxUS are that BAXUS is operated within a series of search space projections on different dimensionalities and BAxUS also includes the local search idea from TuRBO within its optimization process.
Local Region Formulation. As described in Section 2.4, the core idea of BAxUS is to perform optimization in the target space V of dimension dy rather than in the original search space X' of dimension d (d 2 dy). Therefore, to incorporate BAxUS into our proposed CMA-based mcta-algorithm, we need to compute CMA's local regions in the targel space V. Note that since BAXUS obtains the objective function evaluations in the original search domain X, thus, using the CMA's update formula in Eq. (4), we can only compute the search distribution in X. To compute the CMA's local regions in the target space V, our main goal is to project the search distribution (m Ex) from X' to V, and then use the projected search distribution (m Ev) to constrruct the local regions. When performing function evaluation X, BAXUS uses a sparse embedding matrix Q V + X, so that for any vector U E V, we can compute the corresponding vector æ X as Qu, and thus evaluate the objective function value f(x). Our problem is then to find projection matrix P that maps X' to V. This is basically a linear regression problem, and the solution can be derived as P (QTQ) QT (Golub & Van Loan, 2013). Having defined the linear transformation P: - V. given the search distribution (x Ex), we can compute the projected search distribution Ny Ev) as follows (Tong, 1990),
my = Pmx.
Σy PT.
Furthermore, note that BAXUS employs TuRBO as their optimizer, so when defining the local regions for BAXUS, we also make use of the local region adaplation mechanism in TuRBO, which is to include a scale factor L to scale the local region based on the success and failure state of the optimization process. With this, the local region Sy,cHI can be defined as,
Sy. Sy,crA-BAXIS {v ACNA-BAXUS S x²-a,d)
(8)
CHA-BAXUS (v my) is the Mahalanobis distance from E V to the where ACHA-BAXUS (v) V (v my) Darl
scaled search distribution (m ECHA-BAXUS) and EcMA-BAXUS = L'Ey.
Local Optimization. After defining the local region Sy.ck4-BAxus in the target space V, we perform the optimization process similarly to the base algorithm described in Section 4.2, which is to sample a pool of data points from the search distribution (m,, EcM-)BAus) and use BO with the TS acquisition function to pick data points in the target space V. Note that, to make use of the observed data collected in previous (arget spaces, we employ the same splitting strategy as BAXUS to transform the data obtained in previous target spaces into the current target space, and add them to the observed dataset of the current target space. Finally, after collecting À observed data points in the target space V, we then find the corresponding data points in the original search space X' by using the projection 2 = Qu and evaluate their objective function values f(x). From here, we can update the search distribution of CMA, (x Ex), using Eq. (4), and then recompute the local region Sy,cmA-BAxUs via the projected search distribution (m,, Ev). The pseudocode of the local optimization step in -MXUS in Appendix Section A.4, Algorithm 4.
10

---

Published in Transactions on Machine Learning Research (01/2024)
Restart Strategy. Apart from the restart strategy of CMA, we also employ the restart strategy of BAXUS. Specifically, when the local region of the largest target dimension shrinks smaller than the minimum thresh- old, ie., when dy = d and L < Lmin: we restart the embedding with target dimension dy = d and identity embedding matrix Q = Id-
5
Experiments
5.1
Experimental Setup and Baselines
We compare our proposed CMA-based mela-algorithm against a comprehensive list. of related baselines, including BO, TuRBO (Eriksson al., 2019), BAxUS (Papenmeier et al., 2022), LA-MCTS (Wang al., 2020), MCTS-VS (Song ct al., 2022), CMA-ES (Hansen & Ostermeier, 2001), DTS-CMAER (Bajer et al., 2019), and BADS (Acerbi & Ma, 2017). Since LA-MCTS and MCTS-VS are also meta-algorithms like proposed method, we therefore compare against the corresponding methods obtained when incorporating these meta-algorithms with the BO and TuRBO optimizers, resulting in LAMCTS- MCTSVS-B0, NCTSVS-TuRBO. Note that neither LA-MCTS nor MCTS-VS provide guidance on how to incorporate BAXUS as a BO optimizer, so we are unable include the BAxUS-based methods with these two meta-algorithms. Besides, we are also unable to include LAMCTS-BO as its running time prohibitively slow on the problems used in this paper (each repeat takes approximately 3 days run). Details of the experiment setups are in Appendix Sections A.5 and A.6. The average running time for each method are also reported in Appendix Section A.8.
5.2 Synthetic and Real-world Benchmark Problems
We conduct experiments on eight synthetic and three real-world benchmark problems to evaluate all methods. For synthetic problems, use Levy-100D, Alpine-100D, Rastrigin 100D, Elllipoidd-00DD, Schaffer2-100D Branin2-500D, and two modified versions, Shifted-Levy-1000 and Shifted-Alpine-100 100ee For real-world prob- lems, we use Half-cheetab--102,, MassoDNA-188 180D and Rover-100D, These are the benchmark BO problems that used related works including Wang et al. (2016; 2018); Eriksson et al. (2019); Nguyen et al. (2020); Wang et al. (2020); Eriksson & Poloczek (2021); Papenmeier et al. (2022); Song al. (2022); Nguyen et al. (2022); Ziomek Ammar (2023). Details of these problems are in Appendix Section A.7.
5.3
Experimental Results
5.3.1
Comparison against State-of-the-ard BO Optimizers
In this section, we aim to evaluate whether the proposed CMA-based mela-algorithm enhances the per- formance of the BO optimizers by comparing the performance of the CMA-based BO methods (CMA-BO, CMA-TURBO, CMA CMA-BAxUS) with the corresponding BO optimizers (BO, TuRBO, BAXUS). In Fig. 2, it can be clearly seen that our proposed CMA-based meta-algorithm significantly enhances the performance of the correspond- ing optimizers. CMA-BO and CMA-TuRBO outperform BO and TuRBO by a very high margin across all 11 bench- mark problems. CMA-BAXUS outperforms BAZUS significantly on problems (Levy-1001), Rastrigin-100D, 2chaffer2-100D, Shifted-Allpne-1000 Rover-100D, LassoDNA-180D) and performs similarly 5 problems (Alpine-100D, Ellipsoid-110OD, Shifted--Lev--100D, Branin2-500D, Hall-cheetah-102D), Besides, it is worth noting that, in the shifted functions, both the performance of BAXUS and CMA-BAXUS degrade drastically. This is because the global optimum is not at the center of the search domain and these methods no longer have the advantageous benefit of the sparse embedding technique. However, CMA-BAXUS still outperforms BAXUS Shifted-Alpine-1000 and has similar performance in Shfftad-Levy-100D.
5.3.2 Comparison against State-of-the-arr Meta-algorithms
Here, we aim lo evaluate whether the proposed CMA-based mela-algorithm is better than existing state-of- the-art algoritggsss Specifically, we compare with other methods created by applying existing meta- algorithms to the associated BO optimizers (LAMCTS-TURBO, MCTSVS-BO, aCTSVS-TuRBO) In Fig. 3, compared to existing meta-algorithms (LA-MCTS and MCTS-VS), our proposed CMA-based algo-itgo also outper- forms these taae-of-thheart meta-algorithms significantly. It. can be clearly seen that. CMA-TuREE outperforms
11

---

Published in Transactions on Machine Learning Research (01/2024)
Alpine-100D
Levy-100D
Shitted-Alplne1000
Shifted-Levy-1000
10a
Regret
10a
102
10²
10a
101
10
500
1000 1500
2000
500
1000
1500 2000
0
500
1000
1500
2000
500
1000
1500 2000
Number of evaluations
Number of evaluations
Number evaluations
Number of evaluations
Branin2-5008
Schaffer2-100D
Ellipsoid-100D
Rastrigin-1000
10
Regret
101
10c
10
10-2
10a
10?
250
500
750
1000
250
500
500
1000
1500
2000
500
1000 1500 2000
Number of evaluations
Number of evaluations
Number of evaluations
Number of evaluations
Best value found
Rover-100D
1c3
Haff<chetthh-102D
1e-1
LassoDNA-1800
CMA-BO
2.5
CMA-TURBO
3.2
0.0
CMA-BAXUS
-5
3.0
BO TuRBO
-2.5
250 500
750
1000
500
1000
1500 2000
250
500
750
1000
BAXUS
Number of evaluations
Number of evaluations
Number of evaluations
Figure 2: Comparison between the CMA-based BO methods (CMA-BO, CMA-TURBO, CMA-BAXUS) against the original BO optimizers (BO, TuRBO, BAxUS). Plotting the mean and standard error over 10 repetitions. The CMA-based BO methods outperform their respective BO optimizers in most cases.
Alpine-100D
Levy-100D
Shifted-Alpine-1100
Shifted-Levy-1000
10a
Regret
10?
102
10a
10a
101
10a
500
1000
1500
2000
500
1000
1500
2000
0
500
1000
1500
2000
500
1000
1500 2000
Number of evaluations
Number evaluations
Number evaluations
Number evaluations
Branin2-500D
Schaffer2-100D
Ellipsoid-100D
Rastrigin-100D
10a
101
10
Regret
10
10-3
10-2
10a
10a
10-6
250
500
750
1000
250
500
750
1000
1500 Number evaluations
500
1000
2000
500
1000 1500 2000
Number of evaluations
Number of evaluations
Number of evaluations
Best value found
Rover-100D
1c3
Half-cheetah-102D
1e
LassoDNA-
CMA-BO
2.5
CMA-TURBO
3.2
0.0
CMA-BAXUS
-5
2.5
3.0
LAMCTS-TURBO
MCTSVS-BO
250
500
750
1000
500
1000
1500 2000
250
500
750
1000
MCTSVS-TURBO
Number of evaluations
Number of evaluations
Number of evaluations
Figure Comparison between the CMA-based BO methods (CMA-BO, CMA-TURBO) against existing meta- algorithms (LAMCTS-TURED, MCTSVS-BO, MCTSVS-TURBO). Plotting the mean and standard error over 10 repe- titions. The CMA-based BO methods outperform the other meta-algorithms given the same BO optimizer.
both LAMCTS-TURBO and MCTSVS-TuRBO by a very high margin on all of the problems. Similarly, CMA-BO also outperforms MCTSVS-BO on all of the problems by a very high margin. Note that, as mentioned in Section 5, we are unable include LAMCTS-BO due to its prohibitively slow running time (approximately days per one repeat). Furthermore, these meta-algorithms do not suggest on how to incorporate BAXUS as a BO optimizer, we are also unable to compare them with CMA-BAzUS.
12

---

Published in Transactions on Machine Learning Research (01/2024)
5.3.3
Comparison against other Related Baselines
Alpine-100D
Levy-100D
Shifted Alpnne---pnnee-
Shifted-Levy-1000
10a
Regret
102
102
10
10a
101
10
500
1000 1500
2000
500
1000
1500 2000
1500 Number evaluations
500
1000
2000
500
1000
1500 2000
Number of evaluations
Number of evaluations
Number of evaluations
Branin2-500D
Schaffer2-100D
Ellipsoid-100D
Rastrigin-1000
10
10a
10
Regret
10
10
10-2
10a
10
1
750 Number of evaluations
250
500
1000
250
500
750
1000
1500 Number of evaluations
500
1000
2000
500
1000 1500 2000
Number of evaluations
Number of evaluations
Best value found
Rover-100D
1e3
Half<cheetah-102D
1e
1
LassoDNA-180D
CMA-BO CMA-TURBO
2.5
3.2
0.0
CMA-BAXUS
-5
3.0
CMAES
2.5
DTS-CMAES
250 500
750
1000
500
1000
1500 2000
0
250
500
750
1000
BADS
Number of evaluations
Number of evaluations
Number of evaluations
Figure 4: Comparison between the CMA-based BO methods (CMA-BO, CMA-TURBO, CMA-BAXUS) against the CMA-based ES methods (CMA-ES, DTS-CMAES) and BADS, a global optimization method which combines BO and CMA-ES. Plotting the mean and standard error over 10 repetitions. The CMA-based BO methods outperform CMA-ES, DTS-CMAES and BADS consistently.
We compare the performance of our proposed CMA-based meta-algorithm with other related methods such as the CMA-based ES methods (CMA-ES, DTS-CMAES) and BADS, a global optimization method combining BO and CMA-ES. Fig. demonstrates that the CMA-based BO methods outperform the related CMA-based ES methods in the FA literature such CMA-ES and DTS-CMAES. This improvement could be attributed to the use of BO optimizers to select data points for the CMA strategy, instead of randomly sampling, is the case with these evolutionary algorithms. This approach makes the CMA-based BO methods to be more data-efficient. BADS's good performance on the Ellipsoid-100DD problem is thanks to the directed search mechanism based on the mesh points, making BADS perform well on unimodal functions. However, BADS performance is still similar to CMA-BAXUS within given budget. Apart from that, BADS seem to struggle with the high high-dimensional optimization problems with limited data and perform poorly on most of our problems.
5.4
Analysis of the Effectiveness of the CMA-based Meta-algorithm
5.4.1
The Trajectory of the Local Regions by the CMA-based Meta-algorithm
We conduct a study to understand the trajectories of the local regions defined by the CMA-based meta- algorithm in various 2D problems. Note that we only plot the local regions for two derived CMA-basod BO methods, CMA-BO and CMA-TURBO, as the behavior of CMA-BAXUS in 2D problems is similar to that of CMA-TURBO. In Fig. 5, we show the local trajectories for the problem Shifted-Alpine-2D. Additional results for all other synthetic problems can be found in Appendix Section A.9.
We can see that at the beginning (Iteration 0), when the prior information is insufficient for BO, the selected data points scatter randomly throughout the scarch domain. In the later iterations, owing to the use of BO or TuRBO combining with the local regions defined by the CMA strategy, the selected data points converge closer to the global optimum.. Note that in these plots, both CMA-BO and CMA-TuRBO start with the same CMA's search distribution, however, the local regions CMA-TURBO are further scaled by a factor of compared to the base local region due to its local region adaptation mechanism. In the example plotted here, the local regions of CMA-BO gradually shrink toward the global optimum whilst the local regions of CMA-TuRBO shrink
13

---

Published in Transactions on Machine Learning Research (01/2024)
lteration 0
Iteratioo
Iteration
Iteration 5
Iteration
Iteration 10
10
10
0
Global minimum by proposed
Iteration
Iteration L-0.1
iteration L-0.0125
Iteration 5 L-0.8
iterationn L-0.05
Iteration L-0,05
Eilipse vector
L-0,8
10-1
reg
Base local gionn
10
10
10
-10
10
-00
0
10
10
10
10
0
10
10
0
10
Figure 5: The local regions' trajectories defined by the proposed CMA-based meta-algorithm when paired with BO (upper row) and TuRBO (lower row) for the Shifted- AlD function. In this case, the local regions of CMA-BO gradually move towards the global minimum the objective function whilst the local regions of CMA-TuRBO quickly converge sub-optimal location, then restart and move loward the global optimum.
much faster, converge to local optimum at Iteration (with 0.0125), then restart and ultimately converge the global optimum.
5.4.2 How the CMA-based Meta-algorithh Approaches the Global Optimum Compared to Baselines
In Section 5.3, we present the performance of our proposed CMA-based BO methods in terms of the best function values found. In this section, we evaluate how close the selected data points by the CMA-based BO methods are the global optimum of the objective function compared to existing baselines. In Fig. 6, plot the Euclidean distance between the selected data points in each iteration and the global optimum of the objective functions for all the methods. Note that we can only evaluate using the synthetic problems as it is not possible to know the global optima of the real-world problems.
AlpD
Levy-100D
Shifted-Alpine-1100
Shifted-levy-1000
Log Distance
4.
4
4
2
3
1000
2000
0
1000
2000
1000
2000
0
1000
2000
Number of evaluations
Number of evaluations
Number of evaluations
Number of evaluations
Ellipsoid-100D
Rastrigin-100D
Schaffer2-100D
Log Distance
2
-5
1000
2000
1000
2000
500
1000
Number evaluations
Number of evaluations
Number Ivvaluatoos
CMA-BO
CMA-TURBO
CMA-BAXUS
MCTSVS-BO
LAMCTS-TURBE
CMA-ES
BO
TuRBO
BAXUS
MCTSVS-TURBO
BADS
DTS-CMAES
Figure 6: The Euclidean distances between selected data points in each iteration and the global optimum. The CMA-based BO methods can guide the search closer to the global optimum than other baselines.
that the data points selected by the CMA-based BO methods come closer to the global optimum than other baselines. Specifically, the data points selected by CMA-BO, CMA-TuRBO, and CMA-BAXUS are closer to the global optimum than those selected by BO, TuRBO, and BAXUS, respectively in all the problems except. Schaffer2-100D) Compared lo other baselines, it is also clear that the CMA-based BO methods can approach closer to the global optimmu all the problems. It's worth noting that for CMA-BAXUS, there are some big jumps in the distance plots which is due to the changes in dimensionality of the target space, similar BAXUS. When the target dimension increases, the search is performed in a higher dimension, it needs more data to find good solution than when in low dimension. However,
From Fig. 6,
14

---

Published in Transactions on Machine Learning Research (01/2024)
BAXUS suffers this change much more than CMA-BAXUS, ie., the jumps of BAxUS are more significant than CMA-BAXUS as CMA-BAXUS benefits from the guidance of the CMA strategy.
5.4.3
How the CMA-based Meta-algorithm Locates Promising Local Regions Compared to Baselines
Alpine-100D
Levy-100D
Shifted-Alpine-1000
Shifted-lev-11000
Log distance
5.0
5.0
5
5
4
2.5
2.5
3
1000
2000
1000
2000
0
1000
2000
1000
2000
Number of evaluations
Number of evaluations
Number of evaluations
Number of evaluations
Ellipsoid-100D
Rastrigin-100D
Schaffer2-100D
Log distance
5
7.5
5.0
0
1000
2000
0
1000
2000
500 Number of evaluations
1000
Number of evaluations
Number of evaluations
CMA-BO
CMA-TURBO
CMA-BAXUS
LAMCTS-TURBO MCTSVS-TURBO
CMA-ES
TURBO
BAXUS
Figure 7: The Euclidean distances between centers of local regions and the global optimumm The CMA-based BO methods can guide the search closer to the global optimum compared to other basclines.
In this section, we investigate the capability of our CMA-based meta-alee algoriito in guiding the search closer to the promising regions that have high probabilities of containing the global optimum. We investigate the movement of the centers of the local regions defined by the methods by plotting their distances to the global optimum. We conduct this study for BO methods that define center for their local regions, i.e., CMA-BO, CMA-TURBO, CMA-BAXUS, TuRBO, BAXUS, LAMCTS-TURBO, MCTSVS-TuRBO, and CMA-ES. For CMA-BO, CMA-TuRBO, OMA-BAXUS, and CMA-ES, we compute the distance between the mcan vectors of the CMA scarch distributions at each iteration and the global optimum. For TuRBO, BAxUS, LAMCTS-TURBO and MCTSVS-TuRBO, we compute the distance between the centers of the hyper-roctangular local regions of these methods (defined by TuRBO) at each iteration and the global optimuum. In Fig. 7, we can see that all the CMA-based BO methods can effectively guide their respective local regions toward the global optimum better than other baselines. CMA-ES, on the other hand, suffers significantly from the curse of dimensionality issue in the high-dimensional sctting, and steers the scarch distribution away from the global optimum. This is likely due to the over- exploration of CMA-ES, which will be discussed in detail in Section 5.4.4. These results further confirm the capability of our proposed CMA-based meta- algoritho in identifying promising local regions.
5.4.4 The Effectiveness of the CMA-based BO Methods versus CMA-ES
The Bias Issue of the CMA Updates in the CMA-based Meta-algorithm. We investigate key difference between the CMA-based BO methods and CMA-ES: the CMA update process in Eq. (4). In CMA-ES, the CMA search distribution is updated using X data points randomly sampled from the CMA multivariate normal search distribution. In our proposed CMA-based meta-algorithm, BO is used to pick A data points from a candidate pool of ne data points sampled from the CMA search distribution. The use of BO might introduce bias updating the mean vector m, covariance matrix C, and step-size a Eq. (4). One possible issue that can arise is that when ne is large, the À selected data points could be located near single point. Even if n. is not large, the CMA search distribution could be concentrated around point, and this could cause premature convergence of the algorithm. However, we argue that in the high-dimensional setting for BO with a limited budget, this bias issuc is not critical. The first reason is that the search space in a high-dimennsinall optimization problem is very large and thus, a standard size of ne (e.g., thousands) is not possible to make the pool of ne data points very dense, and therefore the scenario of selected data points located very close to a single point is rare when the evaluation budget is limited. strategy, so it does not only select data points
The second reason that BO has an exploitation-explraatee with the imatedd function values (exploitation). but also data points with uncertain function values
15

---

Published in Transactions on Machine Learning Research (01/2024)
Alpine-100D
Levy-100D
Shifted-Alpine1100
Shifted-Levy-1000
10a
Regret
102
Regret
Regret
Regret
102
10'
10
10:
100
10a
0
1000
2000
1000
2000
0
1000
2000
1000
2000
Number of evaluations
Number of evaluations
Number of evaluations
Number of evaluations
Branin2-500D
Best value found
Rover-100D
Best value found
1e3
Half-cheetah-1022
Best found
le-1
LassoDNA-180D
Regret
10t
0
2.5
3.25
10-2
0.0
2.5
3.00
500
1000
500
1000
0
1000
2000
500
1000
Number of evaluations
Number of evallations
Number of evaluations
Number of evaluations
CMA-BO
CMA-BO
CMA-TURBO
CMA-TURBO (x3)
CMA-BAXUS
CMA-BAXUS
CMA-BO (x2)
CMA-TURBO (x2)
CMA-BAXUS (x2)
Figure 8: Performance of CMA-BO, CMA-TuRBO and CHA-BAXUS when increasing (double and triple) the mum- ber of sampled data points ne when optimizing acquisition function Overall, the CMA-based BO methods maintain the performance, indicating the methods robustness w.r.t to Reg-
Iteration: 0
Iteration:
Iteration: 22
Iteration:
Iteration: 44
Iteration:
Total eval:
Total eval: 224
Total eval: 411
Total eval: 598
Total eval: 785
Total eval: 972
100
10-1
10-3
-100
-100
0 100
-00
0
100 -10
0
100
-00
0
100 -100 0 100 -100 0
100
Global minimum
CMA-BO
CMA-ES
Figure 9: The 2D plot for Schaffer--100D function projecting the 2 effective dimensions. CMA-ES over- explores while CMA can guide the search to focus on the promising region around the global optimum.
(exploration). In the high-dimensional setting, with a limited evaluation budget, the number of observed data points used to build the GP is even much smaller compared lo the search space size, and this results in a GP with high uncertainty in many areas, making B0 to select data points with some levels randomness.
We conducted some analysis to validate our arguments. First, as discussed in Sections 5.4.2 and 5.4.3, the results from Figs. 6 and 7 show that the selected data points and the identified local regions by our CMA-based BO methods can approach the global optimum faster than other baselines. Furthermore, we also evaluate the robustness of our proposed CMA-based meta-algorithm w.r.t the mumber of sampled data points Ne. We increase to double and triple the value we use in our default setting which is min(100d. 5000), resulting in the values: min(200d, 10000) and min(300d, 15000). In Fig. 8, we plot the results our proposed CMA-based BO methods with different values n.. We can see that the performance of our proposed methods remains similar, demonstrating their robustness to the choice of ney and thus the bias issue mentioned above is not critical to the high-dinensional setting we use in this paper.
The Over-exploration Issue of CMA-ES in the High-dimensional Setting. It is worth noting that, in practice, CMA-ES tends to over-explore the search space due to its random sampling strategy when selecting data points to update the scarch distribution. This behavior can alrcady be seen in Fig. 6 where we can see that the data points sampled from CMA-ES are very lar from the global optimum within our budget, and in Fig. 7 where it can be observed that the centers of the search distributions by CMA-ES diverge significantly from the global optimum within our evaluation budget. We further illustrate this over-oxploration issuc of CMA-ES on the Schaffer2-10n 100D function by displaying the sampled data points of CMA-ES and CMA-BO across the optimization process. We choose ScDD it has 2 effective dimensions and 98 dummy
16

---

Published in Transactions on Machine Learning Research (01/2024)
dimensions, so can project the selected data points to these two effective dimensions and visualize the selected data points. As we have the results of 10 repeats, we present the first one in Fig. 9, and leave the remaining ones in the Appendix Section A.10, In Fig. 9, we can see that as CMA-ES selects data points randomly and in the high-dimensional setting, these data points are scattered everywhere in the search space, causing CMA-ES over-explore. On the other hand, CMA-BO can select more meaningful data points, mitigating the over-exploration issue of CMA-ES in the high-dimensional setting.
6 Conclusion
In this paper, propose novel CMA-based meta-algorithm to address the high-dimensional BO problem by incorporating a local search strategy and the CMA strategy to enhance the performance of existing BO methods. We further derive the CMA-based ВО algorithms for the cases in which our proposed meta- algorithm is applied to some common state-of-the-art BO optimizers such as BO, TuRBO, and BAxUS. Our extensive experimental results demonstrate the effectiveness and efficiency of the proposed CMA- meta- algorithm, which can significantly improve the BO optimizers and outperform other state-of-tho-at meta- algorithms and related methods.
Acknowledgmentt
This rescarch is supported by Australian Research Council Discovery Project DP220103044. The first and second authors (L.N. & H.H.) would like to thank the Google Cloud Research Credits Program for the computing resources on this project.
References
Abbas Abdolmaleki, Bob Price, Nuno Lau, Luis Paulo Reis, and Gerhard Neumann. Deriving and improv- ing CMA-ES with information geometric trust regions. In Proceedings of the Genetic and Evolutionary Computation Conference, pp. 657-664. ACM, 2017. doi:
Luigi Acerbi and Wei Ji Ma. Practical bayesian optimization for model fitting with bayesian adaptive direct search. In Advances Neural Information Processing Systems (NewrIPS), volume 30, 2017.
Youhei Akimoto, Yuichi Nagata, Isao Ono, and Shigenobu Kobayashi. Bidirectional Relation between CMA Evolution Strategics and Natural Evolution Strategies. In Proceedings of the 11th International Conference on Parallel Problem Solving from Nature: Part pp. 154-163, Berlin, Heidelberg, 2010. Springer- Verlag. TSBN 3642158439.
Charles Audet and J. E. Dennis. Mesh adaptive direct search algorithms for constrained optimization. SIAM Journal on Oplimizaliin, 17(1):188 217, 2006.
Anne Auger and Nikolaus Hansen. A restart ema evolution strategy with increasing population size. In 2005 IEFE congress evolutionary computation, volume 2, pp. 1769-1776. IEFE, 2005.
Lukáš Bajer, Zbynök Pitra, Jakub Repický, and Martin Holeña. Gaussian process surrogate models for the ema-es. In Proceedings of the Genelic and Evolulionary Compulation Conference Companion, pp. 17-18. Association for Computing Machinery, 2019. ISBN 9781450367486. doi: 10. 10.1145//311661933326644
James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for tyyper-prameeee opti-
mization. In Advances in Newral Information Processing Systems 24, pp. 2546-255 2254,, 2011.
Mickael Binois and Nathan Wycoff. A survey on high-dimensional gaussian process modeling with application to bayesian optimization. ACM Transactions Evolutionary Learning and Optimization, 2(2):1-26, 2022.
Eric Brochu, Vlad M. Cora, and Nando de Freitas. A Tutorial on Bayesian Optimization of Expensive Cost Functions, with A Application to Active User Modeling and Hierarchical Reinforcement Learning, 2010.
17

---

Published in Transactions on Machine Learning Research (01/2024)
Dirk Buche, Nicol N Schraudolph, and Petros Koumoutsakos. Accelerating evolutionary algorithms with gaussian process fitness function models. IEEE Transactions on Systems, Mon, and Cybernelics, Part C (Applications and Reviews), 35(2):183 194, 2005.
Roberto Calandra, André Scyfarth, Jan Peters, and Marc Peter Deisenroth. Bayesian optimization for learning gaits under uncertainly: An experimental comparison on a dynamic bipedal walker. Annals of Mathematics and Artificial Intelligence, 76:5 23, 2016.
David Eriksson and Matthias Poloczok. Scalable Constrained Bayesian Optimization. In Proceedings of The 24th International Conference on Artificial Intelligence and Stalistics, volume 130, pp. 730 738. PMLR, 2021.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable Global Optimization via Local Bayesian Optimization. In Advances in Neural Information Processing Systems, pp. 5496-5507, 2019.
Peter I. Frazier, Warren B. Powell, and Savas Dayanik. The Graaieng Policy for Correlated Normal Beliefs. INFORMS J. Comput., 21(4):599 613, 2009. doi: 110100//00001000003111..
Lukas P. Fröhlich, Melanie N. Zeilinger, and Edgar D. Klenske. Cautious bayesian optimization for efficient and scalable policy search. In Proceedings of the 3rd Annual Conference on Learning for Dynamics and Control, L4DC 2021, 7-8 June 2021, Virtual Event, Switzerland, volume 144 of Proceedings of Machine Learning Research, pp. 227 227-240, 2021.
Jacob R. Gardner, Chuan Guo, Kilian Q. Weinberger, Roman Garnett, and Roger B. Grosse. Discovering and Exploiting Additive Structure for Bayesian Optimization. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, volume 54, pp. 1311-1319. PMLR, 2017.
Jacob R. Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew Gordon Wilson. GPyTorch: Blackbox Matrix-Matris Gaussian Process Inference with GPU Acceleration. In Advances in Neural Information Processing Systems 31, pp. 7587-7597, 2018.
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.
Roman Garnett, Michael A. Osborne, and Philipp Hennig, Active learning of linear embeddings for gaussian processes. In Proceedings of the 30th Conference on Uncertainty in Artificiol Intelligence, pp. 230-239. AUAI Press, 2014.
Gene H. Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
Nikolaus Hansen. The CMA Evolution Strategy: A Tutorial, 2016.
Nikolaus Hansen and Anne Auger. CMA-ES: Evolution Strategies and Covariance Matrix Adaptation. In Proceedings of the 13th Annual Conference Companion on Genetic and Evolutionary Computation, pp. 991-1010. Association for Computing Machinory, 2011. ISBN 9781450306904. doi: 10.1145/2001858. 2002123.
Nikolaus Hansen and Andreas Ostermeier. Completely Derandomized Self- Adaptation in Evolution Strate- gics. Evolutionary Computation, 152):159-195, 2001. doi: 10.118/100365500150000033
José Miguel Hernándoz-Lobato, James Requeima, Edward O. Pyzer-Knapp, and Alán Aspuru-Guzik. Parallel and Distributed Thompson Sampling for Large-scale Accelerated Exploration of Chemical Space. Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 1470-1479. PML.R 2017.
Hisham Husain, Vu Nguyen, and Anton van den Hengel. Distributionally robust bayesian optimization with p-divergences. In seirth- Conference Neural Information Processing Systems, 2023.
18

---

Published in Transactions on Machine Learning Research (01/2024)
Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown. Sequential Model-Based Optimization for General Algorithm Configuration. In Learning and Intelligent Oplimization 5th International Conference, volume 6683, pp. 507 523. Springer, 2011. doi: (001007/978-3--42-255533 0-01444222556444400 40.
Rodolphe Jenatton, Cedric Archambeau, Javier González, and Matthias Seeger. Bayesian Optimization with Tree-structured Dependencies. In Proceedings of the 34th International Conference on Machine Learning, volume 70, pp. 1655-1664. PMLR, 2017.
Donald R. Jones. A Taxonomy of Global Optimization Methods Based on Response Surfaces. J. Glob. Optim., 21(4):345 383, 2001. doi: 10.1023/A:10122710255533
Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient Global Optimization of Expensive Black-Box Functions. J. Glob. Optim., 13(4):455-45) 492, 1998. doi: 10.003/AA0000330441111.
Kirthevasan Kandasamy, Jefl G. Schneider, and Barnabás Póczos. High Dimensional Bayesian Optimisa- tion and Bandits via Additive Models. In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pp. 295-304. JMLR, 2015.
Kirthevasan Kandasamy, Akshay Krishnamurthy, Jeff Schneider, and Barnabas Poczos. Parallelised bayesian optimisation via thompson sampling. In Proceedings of the Twenty- First International Conference on Artificial Intelligence and Statislics (AISTATS), volune of Proceedings of Machine Leurning Research, pp. 133-1 142, 2018a.
Kirthevasan Kandasamy, Willie Neiswanger, Jeff Schneider, Barnabás P6czos, and Eric P. Xing. Neural Ar- chitecture Search with Bayesian Optimisation and Optimal Transport. In Advances in Neural Information Processing Systems, pp. 2020 2029, 2018b.
H. J. Kushner. A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. Journal of Basic Engineering, 86(1):97-100 1964. ISSN 0021-9223. doi: 10.1115/1. 3653121.
Ben Letham, Roberto Calandra, Akshara Rai, and Eylan Bakshy. Re-examining linear embeddings for high- dimensional bayesian optimization. In Advances in Neural Information Processing Systems (NeurIPS), volume 33, pp. 1546-1558, 2020.
llya Loshchilov and Frank Hutter. CMA-ES for Hyperparameter Optimization of Deep Neural Networks, 2016.
Prasanta Chandra Mahalanobis. On the generalized distance in statistics. Proceedings of the National Institute of Sciences (Calcutta), 2:49-55, 1936.
Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas. The Application of Bayesian Methods for Seeking the Extremum. Towards Global Optimization, 2(117-129):2, 1978.
Sarah Müller, Alexander von Rohr, and Sebastian Trimpe. Local policy search with Bayesian optimization. In Advances in Neural Information Processing Systems 34, pp. 20708-20722 2021.
Rémi Munos. Optimistic Optimization of Deterministic Function without the Knowledge of its Smoothness. In Advances in Neural Information Processing Systems 24, pp. 783-791, 2011.
Amin Nayebi, Alexander Munteanu, and Matthias Poloczek. A Framework for Bayesian Optimization in Em- bedded Subspaces. In Proceedings of the 36th International Conference on Machine Learning, volume 97, pp. 4752-4761. PMLR, 2019.
Quan Nguyen, Kaiwen Wu, Jacob Gardner, and Roman Garnett, Local Bayesian optimization via maximizing probability of descent. In Advances in Neural Information Processing Systems, volume 35, pp. 13200- Curran Associates, Inc., 2022.
Vu Nguyen, Sebastian Schulze, and Michael Osborne. Bayesian optimization for iterative learning. Advances in Neural Information Processing Systems, 33:9361-999 93771, 2020.
19

---

Published in Transactions on Machine Learning Research (01/2024)
Masahiro Nomura, Shuhei Watanabe, Youhci Akimoto, Yoshihiko Ozaki, and Masaki Onishi. Warm starting cma-es for hyperparameter optimization. In Proceedings of the AAAT Conference on Artificial Intelligence, volume 35, pp. 9188 9196, 2021.
ChangYong Oh, Fistratios Gavves, and Max Welling. BOCK: Bayesian Optimization with Cylindrical Kernels. In Proceedings of the 35th International Conference Machine Learning, volume 80, pp. 3865- 3874. PMLR, 2018.
Leonard Papenmeier, Luigi Nardi, and Matthias Poloczek. Increasing the scope as you learn: Adaptive bayesian optimization in nested subspaces. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), Advances in Neural Information Processing Systems (NeurIPS), 2022.
Jack Parker-Holder, Raghu Rajan, Xingyou Song, André Biedenkapp, Yingjie Miao, Theresa Eimer, Baohe Zhang, Vu Nguyen, Roberto Calandra, Alcksandra Faust, Frank Hutter, and Marius Lindauer. Automatod Reinforcement Learning (AutoRL): Survey and Open Problems. J. Artif. Intell. Res., 74:517-568, 2022. doi: 1at.3//aar...138880
Santu Rana, Cheng Li, Sunil Gupta, Vu Nguyen, and Svetha Venkatesh. High dimensional bayesian opti- mization with elastic gaussian process. In International conference on machine learning, pp. 2883 2891. PMLR, 2017.
Carl Edward Rasmussen and Christopher K. 1. Williams. Gaussian processes for machine learning. MIT Press, 2006. ISBN 026218258X.
Kenan Schié, Alexandre Gramfort, Joseph Salmon, and Luigi Nardi. Lassobench: A high-dimensional hy- perparameter optimization benchmark suite for lasso. In International Conference Aulomated Machine Learning, pp. 2-1. PMLR, 2022.
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proc. IEEE, 104(1):148-175, 2016. doi: 10.1109/
JPROC.2015.2442218
Jasper Snock, Hugo Larochelle, and Ryan P Adams. Practical Bayesian Optimization of Machine Learming Algorithms. In Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Md. Mostofa Ali Patwary, Prabhat, and Ryan P. Adams. Scalable Bayesian Optimization Using Deep Neural Networks. In Proceedings of the 32nd International Conference on Machine Learning, volume 37, pp. 2171 2180. JMLR.org, 2015.
Lei Song, Ke Xue, Xiaobin Huang, and Chao Qian. Monte carlo tree search based variable selection for high dimensional bayesian optimization. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (ods.), Advances Neural Information Processing Systems, 2022. URL cttps://operrevie net/ forum?id=SUzPPe
Jost Tobias Springenberg, Aaron Klein, Stefan Falkner, and Frank Hutter. Bayesian Optimization with Robust Bayesian Neural Networks. In Advances in Neural Information Processing Systems 29, pp. 4134- 4142, 2016.
Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias W. Seeger. Gaussian Process Optimiza- tion in the Bandit Setting: No Regret and Experimental Design. In Proceedings of the 27th International Conference on Machine Learning, pp. 1015-1022. Omnipress, 2010.
William R Thompson. On the likelihood that on unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3-4):285-294 1933. ISSN 0006-3444, doi: 10.1093 10.1093/biomett/
Yung Liang Tong. Fundamental properties and sampling distributions of the multivariate normal distribution. Springer, 1990.
20

---

Published in Transactions on Machine Learning Research (01/2024)
Ryan Turner, David Eriksson, Michael McCourt, Juha Kiili, Ecro Laaksonen, Zhen Xu, and Isabelle Guyon. Bayesian Optimization is Superior to Random Search for Machine Learning Hyperparameter Tuning: Analysis of the Black- Box Optimization Challenge 2020. In NewrIPS 2020 Competition and Demonstrution Track, volume 133, pp. 3-26. PMLR, 2020.
Xingchen Wan, Vu Nguyen, Huong Ha, Binxin Ru, Cong Lu, and Michael A Osbornc. Think global and act local: Bayesian optimisation over high-dimensional categorical and mixed search spaces. International Conference on Machine Learning (ICML) 38, 2021.
Xingchen Wan, Cong Lu, Jack Parker-Holder, Philip J Ball, Vu Nguyen, Binxin Ru, and Michael Osborne. Bayesian generational population-base training. In International Conference on Automated Machine Learning, pp. 14 1. PMLR, 2022.
Linnan Wang, Rodrigo Fonseca, and Yuandong Tian. Learning Search Space Partition for Black-box Op- timization using Monte Carlo Tree Search. In Advances in Neural Information Processing Systems, vol- ume 33, pp. 19511-19522 Curran Associates, Inc., 2020.
Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale Bayesian optimization in high-dimensional spaces. In International Conference on Artificial Intelligence and Statistics, pp. 745- 754. PMLR, 2018.
Ziyu Wang, Babak Shakibi, Lin Jin, and Nando de Freitas. Bayesian Multi-Scale Optimistic Optimization. In Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volune 33, pp. 1014. JMLR.org, 2014.
Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando De Feitas. Bayesian optimization a billion dimensions via random embeddings. Journal of Artificial Intelligence Reseavch, 55:361 -377, 2016.
Ya-xiang Yuan. A review of trust region algorithms for optimization. In Ieiam, volume 99, pp. 271-271121 -222 2000.
Juliusz Krzysztof Ziomek and Haitham Bou Ammar. Are random decompositions all we need in high dimensional bayesian optimisation? In International Conference on Machine Learning, pp. 43347 43368.
PMLR, 2023.
A
Appendix
A.1 TurBO
In TuRBO, in each iteration, the trust region (TR) is constructed as a hyper-rectangle centered at the optimum data point found so far. Each side length of the TR is initialized with a base side length L, and then scaled with the GP lengthscales in each dimension, while maintaining the overall hyper-volume. The size of the TR is critical, as it needs to be large enough to contain potential solutions, while being small enough ensure the accuracy of the GP surrogate model. Therefore, TuRBO adopts an adaptation mechanism to expand shrink the TR depending on whether the algorithm succeeds fails to find better solution. When TuRBO succeeds in finding better solutions for 7suce consecutive times, the TR incrcases its current size, whercas it decreases its size after Thail consecutive failures. Furthermore, TuRBO also defines minimum and maximum thresholds, denoted as Lmin and Lmax respectively, for the TR base side length. The upper bound Lmax is to prevent the TR from becoming too large, while the lower bound Lminn serves as restart criterion, such that when the side length L Lmius TuRBO discards the current TR and restarts new TR from scratch. These hyperparameters (e.g., Tsuce: Thail Tuminy Trmax) as some fixed values in TuRBO.
21

---

Published in Transactions on Machine Learning Research (01/2024)
A.2 BAxUS
BAXUS starts with low value of the larget dimension dv, and then conducts BO process to search for the optimum in this target space within specific evaluation budget before increasing the target cimensionn Additionally, BAxUS also employs TuRBO as its BO optimizer to perform optimization for high-dimensional problems. Therefore, in each iteration, BAXUS also constructs a TR and applies the TR adaptation mech- anism, similar TuRBO, when optimizing within the target space. BAxUS keeps most of the settings to be the same with TuRBO (e.g,, hyper-rectangle TR, shrinkage factor, success tolerance), but redefines the failure tolerance Ttai to make the search quicker in low-dimensional target spaces.
A.3 Additional Information for the CMA Update Formula
Here, we provide additional information about how to set the hyperparmeters for the CMA formula in Eq. (4) based Hansen & Ostermeier (2001). Given the problem dimension and the population size
A = 4 + [3 In dJ, let us define some additional terms as,
w/ = In X+1 Ind, for I,...,,, 2
Pef = ( w
E wf?
(9)
P 2i= w()²
2i= u+ wl 2
The learning rates coefficient are follows,
Cm
2
9 (d+ - 1.3)2 + Heff
(10)
Meft 2 + 1/ Meff
min 1 (1-a c1.
(d 2)² + Helt
The weight coefficients w; is set as,
1
Σ 1= w w!
if w 0,
W
min
1 +
-
C1
1+
24aa
1 1- cf
wi if wt < 0.
(11)
Cμ.
+
Heff
TCCx
Eianti wy
Regarding the the covariance matrix update (sccond line in Eq. (4)), in practice, the evolution path p() Ei (m m(i-1))/o(() is computed via exponential smoothing. Initialized with (( = 0, the exact formula of the evolution path is as follows,
(t) m (c-1)
(()) = (1 (t-c)) + /ccl2 cc olt-1)
(t-1)
(12)
where
+ Helt/ dd Co = d 4 2
(13)
Regarding the 3 coelficient in the step size update (third line in Eq. (4)), the exact formula is as follows,
300) = exp
Co di
p
-1
-1)
(14)
Vd
22

---

Published in Transactions on Machine Learning Research (01/2024)
where
2+ Heft
Ca + 5 2peff
d, = + 2 max
0, d a+11 -1
Poff
+ Co:
(15)
p (1 co)p! + v c/l2 Cc))ann
m
m (t-1) o(t-1) with p 0.
A.4 Pseudocode of the CMA-based BO Algorithms
We present the pseudocode for the local optimization steps of CMA-BO, CMA-TuRBO and CMA-BAxUS. These are the detailed implementation of line 12 in Algorithm depending the BO optimizer bo bo_opt. Note that, discussed in the base algorithm in Section 4.2, when performing the local optimization step, first need to sample a pool of ne data points following the previous scarch distribution, then perform BO to select the data points. In practice, we set CMA-BD and CMA-TuRBO or de = dy, the current target dimensionalltty, for CMA-BAXUS. This value of n. the same as in TuRBO when selecting sampling data points for the TS acquisition function.
min 5000) where d, = d, the dimensionality of the problem, for
Algorithm Local Optimization for CMA-BO.
1: Input: Objective function f(), search distribution ((m, Σ), local region S, dataset Q number of sampling points Nee
Output: A new observed data {xnextt Mnext
Train GP from Ω
Sample ne data points A = = from N(m, Σ) and constrained within S Propose next observed data Tnexx arg minl o TS (x)
Evaluate the observed data Ynext (2next +
Return (mnextt Ynext
Algorithm Local Optimization for CMA -MA-TURBO.
1: Input: Objective function f()), search distribution M(m, ΣCHA-TUBED), local region SCHA-TURBO, dataset Z number sampling points ne Output: A new observed data (mnexts )next
Train GP from f2
Sample ne data points A= (=;))= from N(m, ECHA-TUREO) and constrained within SCHA-TURBO
5: Propose a next observed data Vnext arg min a (x))
6: Evaluate the observed data Ynext. = I (xnext) +
Return {xnext )next)
Algorithm 4 Local Optimization for CMA -BAXU..
1: Input: Objective function f(),, search distribution (m DOMA-BARUS). local region -y,cuA-Baxus, dataset
Ω number of sampling points ne, embedding matrix Q V + X'
2: Output: A new observed data {ænext, Unext, fnext) } in both X' and V
Train GP from {0; y:) Ω
[2]
Sample ne data points A = {wj}īm from A(my, Eema-BAuu and coustrained within Sy BAx Propose next observed data Unext = arg min a TS (v) and Tnext
Uuext
Evaluate the observed data ynext [2xnext +
Return { (xexxt Unext, fnext) }
23

---

Published in Transactions on Machine Learning Research (01/2024)
A.5 Experimental Setup
We use Matérn 5/2 ARD kernels for the GPs in all methods. The input domains of all problems are scaled to have equal domain lengths in all dimensions as Loshchilov & Hutter (2016). The output observations are normalized following a Normal distribution y ~ -V(0, 1).
For the hyperparameters of the CMA strategy in all CMA-based BO and ES methods, we set them using the suggested values in Hansen (2016). Specifically, the population size À is set to be 4 + [3 In dJ. The initial mean vector m (²)) is selected by minimizing 20 initial data points following Latin hypercube sampling (Jones, 2001). The covariance matrix (' is initialized with an identity matrix Ia. and the initial step size ((0) is set to 0.3(u 1) where u, I denote the upper and lower bounds of the search domain X, i.c., X = [2, u]d.
To ensure fair comparison between the CMA-based BO methods and the corresponding BO optimizers, we set the hyperparameters of the CMA-based BO methods to be the same as those of the corresponding B0 optimizers. Specifically, for BO and CMA-BD, the hyperparameter settings of the GP and the TS acquisition function of these met hods are the same. For TuRBO and CMA-TURBO, we follow the same setting suggested by TuRBO (Eriksson et al., 2019) to set the initial TR base side length Lo, the maximum and minimum TR side lengths Lmax and Lmin, and the success and failure threshold Tanoe and Thail For BAXUS and CMA- CMA-BAXUS, we also follow the same setting suggested by BAXUS (Papenmeier et al., 2022) to set the initial TR base side length Lo, the maximum and minimum TR side lengths mmax and Limin, the success and failure thresholds Tsuce and Thail and the bin size b. All the developed CMA-based BO methods (CMA-BO, CMA-TuRBO, CMA-BAxUS) are implemented using GPyTorch (Gardner al., 2018) with TuRBO and BAXUS. All the Python- based methods are run with the same Python package versions.
A.6 Baselines
To evaluate the baseline methods described in Section 5, we use the implementation and hyperparameter settings provided in the authors' public source code and their respective papers. Note that for DTS- and BADS, the authors' implementation source code is in Matlab, so to ensure consistency in the objective function evaluation process with other baselines, we call Python from Matlab to evaluate the objective function values. All the methods are initialized with 20 initial data points and are run for 10 repeats with different random seeds. All experimental results are averaged over these 10 independent runs. We then report the mean and the standard error of the simple regret or the best optimal value found. Details of the implementation for each baseline in the paper are as follows.
BO. This is the standard BO method with the TS acquisition function. The GP is constructed with the Matérn 5/2 ARD kernel and is fitted using the Maximum Likelihood method. The domains of the input variables in all problems are scaled to have equal domain lengths in all dimensions (Loshchilov & Hutter, 2016). The output observations }y:} are normalized following a Normal distribution /(00 1).
TuRBO (Eriksson et al., 2019). We set all the hyperparameters of TuRBO suggested in their paper. This includes the upper and lower bound for TR side length Lunax 1.6, Liii = 2-7, batch size b = 1 and the TR adaptation threshold Tancn 3, Thail mmax (4/b,d/b)] where d is the dimension of the problem. We use their implementation that is made available at attps://gitann woweuber--eseeere//tteee research/TuRB0.
BAXUS (Papenmeier al., 2022). We set all the hyperparameters of BAXUS as suggested in their paper. This includes the upper and lower bound for TR side length T.max 1.6, Timin 2-7 the TR adaptation threshold Tsuce = 3 and bin size b b=3, 3, budget to input dim mp is set to the maximum budget. We use their implementation that is made available at /etps://githhe com/LeoIV/BAXUS.
LA-MCTS (Wang et al., 2020). We set all the hyperparamcters of LA-MCTS as suggested in their paper. This includes the exploration factor in UCB C, 1, the kernel type of SVM is RBF and the splitting threshold 0 20. For Levy function, use different settings, which is recommended in the author's
24

---

Published in Transactions on Machine Learning Research (01/2024)
implementation code', i.., C 10, polynomial kernel, 0 = 8. We use their implementation that is made available at hetps://githuue cc:///accbboorrrearr//aacCCCC
MCTSVS (Song et al., 2022). We set all the hyperparameters of MCTS-VS as suggested in their paper. This includes the exploration factor in UCB Cp 1, the fill-in strategy of "best-k' with k 20, the feature batch size N, 2, the sample batch size N, 3, the tree re-initialization threshold Niad 5, the node splitting threshold Naplu = 3. We use their implementation that is made available at https: //g.thub.com/lladaa
CMAES (Hansen & Ostermeier, 2001). We use the default settings as suggested in the paper, which is similar to our settings for CMA-BO. This includes the population size À = 4 + |3 + In d where d is the problem dimension, the random initial mean vector m selected from the minimum of the 20 random initial points, the identity initial covariance matrix C = Ia and the initial step-size g = 0.3(u - 1b) where the domain is scaled to uniform bound of u|d. We also activate the restart mechanism of CMA-ES so that the algorithm can restart when converges local minimum. We use their implementation that is made available at https: bttpps//ggttub.... com www.cmammmccceeec. CMAAEM//M-Maaaa
DTS-CMAES (Bajer et al., 2019). We set all the hyperparameters of DTS-CMAES as suggested in their paper. We use the doubly-trained GP configuration with the population size X 8 + |6 + In d] where d is the problem dimension, initial mean vector (") selected from the minimum of the 20 random initial points, the identity initial covariance matrix ( Ia and the initial step-size σ()) 0.3(u lb) where the domain is scaled to uniform bound of [4, u]a and fixed learning rate 3 0.05. We use their implementation that is com/bajeeiieeooooooo... waaejeekkuurrrrre -maee...teeem........
made available at sstps://githhe
BADS (Acerbi & Ma, 2017) We set all the hyperparameters of BADS as suggested in their paper and the Matlab package. We use their implementation that is made available at https: //ge:////tttttt/g www. con/acerbilab/ bads.
A.7 Synthetic and Real-world Benchmark Problems
We conduct experiments on eight synthetic and three real-world benchmark problems to evaluate all methods.
Synthetic Problems. We use Levy-100D, Alpine-100D, Rastrigin- Ellipsoid-100D, Schaffer2-100D. Branin2-500D, and two modified versions, Shifted--Ley--000D and Shifled-Alpine-10000 For Branin2-500D, we use the implementation from Papenmeier et al. (2022); Wang et al. (2016) where the function is created by adding additional 498 dummy dimensions to the original Branin 2D function, resulting in function with the dimension d to be 500. Schaffer2-1000 is implemented similarly with 98 dummy dimensions added to the original Schaffer 2D function. The Levy-1001 Alpine-1001. Rastrigin- 100D and Hiipssoid-100D are common test functions" used in BO research, and we set the dimension d to be 100 for each function. Additionally, for Levy-100D and Alpine-100D, we create two new versions, namely Shifted-Appine 100d and Shifted-Levy, 100D, where we shift the global optimum away from the original global optimum by uniformly random shifting, i.c., we set (s) = (original + 5) with 8 = ji...., 5al E [1, u] 100 and 5; ~ Ull,, u). The scarch domains of these two functions, X' = !1, w]d, are kept the same as in the original functions. The molivation behind including these two shifted synthetic problems for evaluation that, based on our observations, some sparse embedding methods (c.g., BAxUS) have considcrable advantages when the global optimum is at the center of the search domain, thus, we also evaluate all methods on problems where the global optima are not at the search domain's 's center.
Real-world Problems. We use the following real-world problems: Half-cheetah-102D, DASsoDNA-1800 and Rover-100D. For Half-choctah-102D, we use the same implementation as described in Song al. (2022),
functiions/functiinn..y httt:p::/ww..ttt..eeaaarrrraraas
25

---

Published in Transactions on Machine Learning Research (01/2024)
which parameterizes the Half-cheetah-v4 Mujoco environment from the Gym package" into 102D rein- forcement learning (RL) problem. The goal of this problem is to optimize the parameters of linear policy designed to solve the RL task. These Mujoco RL tasks have been used in many works, such as Wang et al. (2020); Nguyen et al. (2020); Song et al. (2022). For LassoDNA-180D, we use the implementation from the Python LassoBench library (Sehić et al.; 2022) as in Papenmeier et al. (2022). The problem LassoDNA-180D solves the Least Absolute Shrinkage and Selection Operator (LASSO) problem using the DNA dataset from microbiology problem. The LassoBench suite has been used in several previous works, such as Papenmeier al. (2022); Ziomek & Ammar (2023). For Rover-100D, (2018). This problem optimizes the locations of 50 points in 2D>plane trajectory of rover, resulting in a 100D benchmark function. The goal is to maximize the reward calculated based on the number of collisions along the rover trajectory. The Rover function has been used in various research works, e.g., Wang el al. (2018); Eriksson et al. (2019); Eriksson & Poloczek (2021); Nguyen et al. (2022).
the implementation provided by Wang al.
A.8 Running time of all methods
The results demonstrate that the running time our proposed CMA-based BO methods is very similar to the running time of the BO optimizers we incorporale. This demonstrate the ellicacy of our proposed mela-algorithm.
We report the average running time per each iteration in Table 1.
Table 1: Average time (in second) for each iteration run in each method.
Average time per iteration (s)
Shifted
Shifted
Alpine 100D
Levy 100D
Ellipsoid
Rastrigin 100D
Schaffer2 100D
Branin2
Alpine
Levy
100D
500D
100D
100D
CMA-BO
3.3
3.28
3.62
3.46
3.21
5.48
3.07
6.8
CMA-TuRBO
3.13
3.17
3.33
3.34
3.24
3.17
2.98
9.19
CMA-BAXUS
8.68
21.27
21.28
21.4
8.27
9.32
5.93
15.19
5.7
5.69
5.63
5.61
5.57
5.6
2.57
1.73
TuRBO
5.32
5.29
5.28
5.53
5.39
5.45
2.46
2.97
BAxUS
9.61
21.79
21.42
20.62
7.56
8.61
7.04
9.36
MCTSVS-BO
0.96
0.36
1.5
1.18
0.48
0.63
0.25
0.43
-CRSVS-TTRBO
0.24
0.27
0.79
0.8
0.27
0.23
0.13
0.21
LAMCTS-TURBO
3.53
3.51
10.03
8.99
3.46
3.45
1.95
2.39
CMA-ES
4.6E-04
3.2E-03
4.4E-04
2.9E-03
4.4E-04
1.7E-03
4.7E-04
1.0E-03
BADS
1.52
1.74
1.78
1.6
1.18
1.7
0.28
0.29
DTS-CMAER
0.06
0.05
0.05
0.05
0.06
0.05
0.03
0.44
A.9
Additional Trajectory Plots of the Local Regions by the CMA-based Meta-algorithi
We Weshow the remaining trajectory plots of the local regions defined by CMA-BO and CMA-TuRBO in the remaining synthetic functions: AD (Fig. 10), Levy-2D (Fig. 11), Shifted-Levy-2D (Fig. 12), Branin-2D (Fig. 13), Ellipsoid-2D (Fig. 14), Schaffer-2D (Fig. 15) and Rastrigin-2D (Fig. 16).
A.10 Additional Results of Schaffer2-100D
As an expansion of Fig. 9, we show in Fig. 17 all the repeats the 2D plot for Schaffer2. 100D functions.
26

---

Published in Transactions on Machine Learning Research (01/2024)
Iteration 0
Iteration 2
Iteration 4
Iteration 7
Iteration 8
Iteration 10
10
10a
10
0
10
-10
10-5
-10
10
10
10
10
10
10
10
10
10
-0
10
Iteration 0
Iteration 2 L=0,2
Iteration 4 L=0.025
Iteration 7 L=0.4
Iteration 8 L=0.4
Iteration 10 L=0.1
L=0.4
10
10a
10-1
0
10-3
-10
10
-10
-10
-10
-10
-0
Global minimum
Ellipse mean vector
Data proposed by BO
Local region
Base local region
Trajectories of the local regions defined by CMA-based BО methods, CMA-BO (upper) and CMA-TURBO (lower), for Alpine 2D function. The function global optimum at [].
Figure 10:
Iteration 0
Iteration 2
Iteration 4
Iteration 5
Iteration 8
Iteration 10
10
101
10-1
0
10 3
-10
10-5
10
10
10
10
-10
10
10
10
10
10
-00
10
Iteration 0
Iteration 2 L=0.05
Iteration 4 L=0.025
Iteration 5 L=0.4
Iteration 8 L=0.05
Iteration 10 L=0.0125
L=0.4
10
101
10-1
0
10
-10
10-5
-10
-10
-10
-10
-10
-10
Global minimum
Ellipse mean vector
Data proposed by BO
Local region
Base local region
Trajectories of the local regions defined by CMA-based BО methods, CMA-BO (upper) and CMA-TuRBO (lower), for Levy 2D function. The function global optimum is [1,
Figure 11:
27

---

Published in Transactions on Machine Learning Research (01/2024)
Iteration 0
Iteration 2
Iteration 4
Iteration 7
Iteration 8
Iteration 10
10
10a
10-1
0
-0 3
-10
10-5
-10
10
10
10
10
10
10
10
10
10
10
10
Iteration 0
Iteration 2 L=0.4
Iteration 4 L=0.05
Iteration 7 L=0,8
Iteration 8 L=0.4
Iteration 10 L=0.1
0.0.8
10
10a
10-1
0
10
-10
10-5
-10
10
-0
-
-00
Global minimum
Ellipse mean vector
Data proposed by
Local region
Base local region
Figure 12:
Trajectories of the local regions defined by CMA-based BО methods, CMA-BO (upper) and
CMA-TuRBO (lower), for Shifted-Levy 2D function. The function global optimum [3, -4].
Iteration 0
Iteration 2
Iteration 5
Iteration 8
Iteration 10
Iteration 12
15
10a
7.5
10
10-3
0
10-5
-5
2.5
10
-
2.5
10
5
2.5
10
-
2.5
10
5
2.5
10
-
2.5
10
Iteration 0
Iteration 2 L=0.4
Iteration 5
Iteration 8
Iteration 10
Iteration 12 L=0.2
L=0.8
L=0.05
L=0,8
L=0.8
15
10
7.5
1
10-3 3
0
10-5
-5
-
2.5
-
2.5
-
2.5
-
2.5
2.5
Global minimum
Ellipse mean vector Local region
Data proposed by
Base local region
Figure 13:
Trajectories of the local regions defined by CMA-based BO methods, CMA-BO (upper) and
CMA-TuRBO (lower), for Branin 2D function. The function has 3 global optima [-m, 12. 12.275), [Tr, 2.275) and 19.42478, 2.475].
28

---

Published in Transactions on Machine Learning Research (01/2024)
Iteration 0
Iteration 2
Iteration 5
Iteration 6
Iteration 8
Iteration 12
10
10a
0
10
10-3
-10
10-5
-10
10
10
10
-10
10
10
10
10
10
-10
10
Iteration 0
Iteration 2 L=0.2
Iteration 5 L=0.0125
Iteration 6 L=0.8
Iteration 8 L=0.4
Iteration 12 L=0.0125
L=0.4
10
10a
0
10-1
10-3
-10
10-5
-10
-10
-10
10
-10
Global minimum
Ellipse mean vector
Data proposed by BO
Local region
Base local region
Trajectories of the local regions defined by CMA-based ВО methods, CMA-BO (upper) and CMA-TURBO (lower), for Ellipsoid 2D function. The function global optimum at [0, 0].
Figure 14:
Iteration 0
Iteration 8
Iteration 20
Iteration 24
Iteration 32
Iteration 40
100
10-1
0
10-2
10-3
-100
-1000
100
1000
100
10000
100
-100
100
1000
100
-1000 100
Iteration 0
Iteration 8 L=0.4
Iteration 20 L=0.1
Iteration 24
Iteration 32
Iteration 40 L=0.025
L=0.4
L=0.4
L=0.4
100
10
0
10-2
10 3
-100
-1000
100
-1000
100
-1000
100 -1000 100 -1000
100 -1000 100
Global minimum
Ellipse mean vector
Data proposed by BO
Local region
Base local region
Trajectories of the local regions defined by CMA-based BО methods, CMA-BO (upper) and CMA-TuRBO (lower), for Schaffer 2D function. The function global optimum at [).
Figure 15:
29

---

Published in Transactions on Machine Learning Research (01/2024)
Iteration 0
Iteration 8
Iteration 14
Iteration 16
Iteration 21
Iteration 30
5
10a
10-1 -1
0
10-3
-5
10-5
-
5
5
5
5
5
-
5
5
5
-
5
Iteration 0
Iteration 8 L=0.1
Iteration 14 L=0,8
Iteration 16 L=0.2
Iteration 21
Iteration 30
L=0.4
L=0.8
L=0.025
5
10
10-1
0
10-3
-5
10-5
-5
-5
5
-5
5 -5
5 -5
-5
5
Global minimum
Ellipse mean vector
Data proposed by BO
Local region
Base local region
Trajectories of the local regions defined by CMA-based BO methods, CMA-BO (upper) and CMA-TuRBO (lower), for Rastrigin 2D function. The function global optimum is at [0, 0].
Figure 16:
30

---

Published in Transactions on Machine Learning Research (01/2024)
5-0, nm37
5-0, 1=11, n-224
5-0, 1-22, n=411
|m33, n
5m0, 1-44,--488 n-785
5-0, 1-55, n
0-
100
5=1, i=0, n-37
1=1, 11. 224
5-1, 1-22 n1411
5- 1. 33 n 598
1, 44, 785
s-1, i-55, n7
100
-100
5=2, i=0,
i=11,
s=2, 1=22, n=411
n= 598
2= n=785
i=55,
100
5-3, 1-0, n-37
5-3, 1-11, 2-2
5-3, 1-22, n-411
5-3, 1-33, nー 598
5-3, 1-44, 7-785
5-3, i-55, n-972
-100
s=4, i=0, n=37
4= i=11, 2=224
$=4, 1=22, n=411
s=4, 4, i=33, n=9
$=4, 1=44, n=785
s=4, i=55, n=972
100
Global minimum
CMA-ES
5=5, 1= 1=1
1= n=411
l=33, 598
5=5, 1=44, n=785
i=55,
100
5-6, 1-0, n-37
6-6, i- 11. 2-224
5-6 1-22, n1411
5- 6, 1-33, nー 598
5-6, i-44, 785
5-6, i-55, n
1100
$=7, I=0, n=37
$=7, |=11,, n=2
s=7, 1=22, n=411
5=7, |=38, n= 598
5=7, |=44, 7=785
5=7, i=55, n-972
100
-100
5m8, |-0, n-377
|-1 1-22, n-411
|m33, 598
n87
1-55,
100
100
9=9, i=0, n=37
9- i-
11. 224
5=9, i=22, n=411
5 9, 3-
n 598
9, i=44, n= 785
9 i=55, 97
0
-100
-100
100
100
100
1.00
100
-100
100
-100
100
-100
100
Figure 17: 2D plots for Schaffer2-1001 projecting on the 2 effective dimensions for 10 repeats. CMA-ES always tends to over explore whilst CMA-BO BD generally finds more meaningful data points.
31


Reviews:
[{"review_id": "OtcoFBo6t7", "reviewer": "Reviewer_6Py4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is reasonably well written but IMO much longer than it needs to be. There is only one key idea being introduced here, and only really two result figures (figs 2 and 3 are redundant and should be combined), so there is no need for 15 pages. All sections would benefit from aggressive editing, but in particular the repetitive descriptions of how CMA is combined with the different underlying BO methods would be better presented as a succinct table. \n\nThese BO methods are in mind a little dissatisfying since they feel a like a collection of heuristics without any theoretical grounding. Indeed, I would have thought CMA ran a risk of being dragged towards wide shallow optima at the expense of deep narrow optima. The empirical results seems strong but it is challenging to assess how well these examples represent the characteristics of real world applications of BO. \n\nMinor comments: \n- what is a \"heterogenous\" objective function? (p1)\n- eq 2 is repetitive with the inline eqns after it\n- TurBO is introduced twice in different sections\n- section 4.1: how is the \"true search distribution\" defined?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider a meta-algorithm for Bayesian optimization (BO) where the underlying BO model is only fit to subregions. They port an approach from evolutionary algorithms, covariance matrix adaptation (CMA), to control the points that are modeled using the underlying BO method. They show promising results on learning synthetic and \"real\" functions.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Edit down to preferably 8, max 10 pages.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "SO3mJpqs3Y", "reviewer": "Reviewer_kmik", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength \n\nThe proposed framework is a rather simple framework and can be applied to many BO baselines. The use of the BO in the surrogate-assisted CMA-ES in this paper is different from other surrogate-assisted CMA-ES, hence I see the novelty. (At the same time, I see the weakness  on this point. See below.) The shown empirical performance is promising. \n\nWeakness\n\n1. The expression of the CMA mechanism is buggy. The authors call $\\sigma$ the overall standard deviation, but it is indeed the overall variance as the authors define $\\Sigma = \\sigma C$ rather than the standard notation of $\\Sigma = \\sigma^2 C$ in the CMA-ES literature. Then, the algorithm description in (4) and below seems wrong. If $\\sigma$ is the variance, (4) is fine but the update of $p$ should not be divided by $\\sigma$. Otherwise it does not make sense. If $\\sigma$ is the standard deviation, (4) is wrong. The second term should be divided by $\\sigma^2$ rather than $\\sigma$ and the third term should not be divided by $\\sigma$. Because of this point, I am not sure how the proposed approach is indeed implemented.\n\n2. The proposed approach breaks one of the main design principle of the CMA-ES, i.e., stationarity. As lambda candidate solutions are selected by BO approaches from the pool of many candidate solutions that are generated from the search distribution, the selected lambda candidate solutions do not follow the Gaussian distribution anymore. Then the update of $m$, $C$, and $\\sigma$ is strongly biased. In the situation of the proposed approach, because the best candidate solutions are selected from BO, if the number of points in the candidate solution pool is large, the selected candidate solutions are located near the single point. Even if the number of the points is not so large, the distribution is expected to be concentrated around a point compared to the original sampling distribution (i.e., Gaussian distribution). Then, $\\sigma$ and $C$ end up converging at nearly the rate of their learning rate. Such an update is known to tend to cause premature convergence to a not-even-locally-optimal point and hence is avoided in the CMA-ES. \n\nTherefore, I feel that the reported efficacy of the proposed approach is at the cost of the final performance on many multimodal problems. Although the test problems have many local minima, it is not clear whether the better performance of the proposed approach is due to the convergence to a better local minima or just due to the faster convergence to the not-even-locally-optimal point. To evaluate the local search ability and the global search ability, the authors should perform more empirical study on unimodal problems (such as ellipsoid) and multimodal problems (such as rastrigin, shaffer, bohachevsky, etc.). One possibility is to use a standard benchmarking framework, such as BBOB,  to compare the proposed approaches with the baselines.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a novel framework for black-box optimization in high dimensional search space. To address the issue of inefficacy of BO, the proposed approach combines BO and CMA in the following way. The search distribution of CMA, i.e., Gaussian distribution, is maintained and many candidate solutions are generated from the search distribution within 3 sigma region. BO then selects the best lambda candidate solutions to evaluate on the objective. The search distribution of CMA is updated using the selected candidate solutions using the same update formula as the original CMA-ES. Overall, the proposed approach can be recognized as a variant of the CMA-ES with a surrogate model, where the surrogate model is used to select candidate solutions to be evaluated. The main difference between existing surrogate-assisted CMA-ES and the proposed framework is that the surrogate model is used to pre-select the good candidate solutions rather than assigning their artificial objective function values.\n\nThe effect of the proposed approach is demonstrated on synthetic problems and three ML-type problems. Different BO approaches with and without the proposed mechanism, as well as the standard CMA-ES and a surrogate-assisted CMA-ES are compared. Compared to the baseline approaches without the proposed mechanisms, the proposed approaches show their superior efficacy with a limited function evaluation budget.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see the weak points pointed out above.", "broader_impact_concerns": "no concern", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "5YIRuuckL6", "reviewer": "Reviewer_pHRT", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- Well written text.\n- Comparison to an extensive suite of baselines.\n- Good coverage of related work with perhaps one exception (Unbounded Bayesian Optimization via regularization).\n\nWeaknesses:\n- Minor worry about the sentence: \"the rationale behind this step is that in the CMA strategy [...] the observed data points are required to be sampled from the previous search distribution.\" This argument somewhat breaks when the points are then filtered by Thompson sampling from the GP. The actual distribution from which the points are sampled is a factored distribution where the first is the CMA's normal distribution and the second is p(x* | D) under the GP prior.\n- Algorithm 1 could be a bit cleaned up, see next section.\n- Experimental results could be a bit cleaned up, see next section.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This manuscript proposes to combine the CMA strategy (from CMA-ES) with that of BO to dynamically adjust trust regions given observed data in order to make high-dimensional BO more tractable. On 5 synthetic objective functions and 3 simulated physics control tasks, they demonstrate that their approach beats many existing baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Algorithm 1: It's not super clear what the different roles of Omega and D_lambda are, perhaps that can be clarified in both the pseudo-code and the text.\n- The 2 experimental results figures are very cluttered and could probably be refactored into a clearer 3-4 plots with more targeted comparisons (e.g. CMA-XYZ vs XYZ then CMA-XYZ vs MCTS-XYZ, etc).", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This manuscript considers the problem of high-dimensional Bayesian optimization, wherein we seek to optimize an expensive objective function that is typically modeled as a "black-box." This is an important and well-studied problem. The authors attack this problem with a strategy that attempts to partition the high-dimensional domain into smaller local regions that can be explored more effectively. This strategy has the advantage of being compatible with many existing approaches and performs admirably in a series of well-designed experiments.

The reviewers universally agree that the work presented in this manuscript is high quality and of interest to the TMLR audience. They are universal in their recommendation of acceptance.

Throughout the review and discussion period, the reviewers provided feedback and suggestions that the authors have faithfully incorporated into the manuscript. I believe these changes have considerably strengthened the work. I commend both authors and reviewers for their engagement throughout this process.

I believe the manuscript is both suitable for and ready for publication.
Decision: Accept (poster)